Franck Pigeonneau

Coupled modelling of redox reactions and glass melt fining processes

, but differs in the numerical method. After a description of the model, we present the evolution of a bubble population with time and also apply the numerical method to an experimental tool where bubbling is used to equilibrate the partial pressures between the bubbles and between the melt. We find, with the model, an equilibrium time longer than seen experimentally. The possible origins of the disagreement are investigated and discussed.

Introduction

Molten soda-lime-silica glass is obtained from, among other things, carbonaceous raw materials such as CaCO 3 and N a 2 CO 3 . These release CO 2 gas which is then found in the molten glass partly in bubble form and partly dissolved in the melt. Nemec and Klouzek [START_REF] Němec | Modelling of glass refining kinetics. Part 1: single bubbles[END_REF] indicate that the bubble concentration can approach 10 9 bubbles per cubic meter. Since the final product must have the fewest bubbles possible, the goal of the fining process is to eliminate a large number of the bubbles present in the glass.

Bubble evacuation is a difficult task in a viscous fluid such as molten glass due to the low rise velocity especially for small bubbles with diameters less than 0.5 mm. This velocity increases with temperature due to decreasing viscosity and with increasing bubble size. Moreover, though the addition of "fining" agents, the bubble size can be increased due to the migration of "fining" gases dissolved in the molten glass towards the bubbles. Finding good fining agents is the art of a glass maker. Optimization of these agents can be achieved through the use of experimental devices and by modelling the fining process.

As previously mentioned, most fining processes are governed by the oxidation-reduction reactions of fining agents. In many investigations [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF] these two phenomena are treated separately with, in the first stage, redox behavior in the glass being determined and, secondly, the evolution of the bubbles being studied. But due to the large number of bubbles, the consumption of dissolved gases can significantly change the redox state. Only a few publications by Nemec and his co-workers [START_REF] Kloužek | Modelling of glass refining kinetics. Part 2: Bubble distribution models and methods of measurement of refining properties[END_REF][START_REF] Stahlavsky | The bubble effect on the redox state of glass[END_REF] try to take into account the coupling of redox and fining. This article carries out the coupled modelling of these two processes in a system assumed to be spatially homogeneous.

After the presentation of the basic equations of the model in section 1, the results are given in section 2 followed by a discussion and conclusion.

Modelling of redox and fining processes

In a glass melt at high temperatures, chemical reactions play an important role in the fining process [START_REF] Doremus | Glass Science[END_REF]. Many gas species, such as water, oxygen, carbon dioxide and sulphur, react with glass and as indicated in the introduction, the fining of gas bubbles is supported by the addition of fining agents which react as the temperature increases. The chemical agents classically used are: sulphate, arsenic(V) oxide, antimony(V) oxide.

Redox model

Knowledge of the concentrations of dissolved gases is necessary in order to describe the evolution of bubbles. These concentrations depend on reduction-oxidation reactions. We assume that locally these reactions are in a state of chemical equilibrium, an assumptionthat has been justified by Russel [START_REF] Rssel | Redox reactions during cooling of glass melts -A theoretical consideration[END_REF]. Consider a polyvalent chemical element, M . The reduction reaction of the polyvalent ion M is given by [START_REF] Schaeffer | Oxidation state of equilibrated and non-equilibrated glass melts[END_REF]

M (n+k)+ + k 2 O 2-⇋ M n+ + k 4 O 2 . (1)
M n+ and M (n+k)+ represent, respectively the reduced and oxidized states of M . The equilibrium constant of this reaction is given by the following equation [START_REF] Beerkens | Mechanism of oxygen diffusion in glassmelts containing variable-valence ions[END_REF]

K M = C M n+ C k/4 O2 C M (n+k)+ . (2)
C M n+ , C O2 and C M (n+k)+ are molar concentrations. The modelling of reduction-oxidation reactions is presented here under a general form where we assume that R reactions take place in the glass melt. The number of ionic species is N I and the number of fining gas species is N f g . The generic reaction r is thus written

NI i=1 ν ′ ri A i + N f g j=1 β ′ rj G j ⇋ NI i=1 ν ′′ ri A i + N f g j=1 β ′′ rj G j . (3)
A i is the ion i where i varies from 1 to N I and G j corresponds to the gas j where j varies from 1 to N f g . The stoichiometric coefficients for ionic and gas species taken into account in reaction r are defined by the following relations [START_REF] Prigogine | Thermodynamique. Des moteurs thermiques aux structures dissipatives[END_REF]

ν ri = ν ′′ ri -ν ′ ri , for i = 1 to N I . ( 4a
)
β rj = β ′′ rj -β ′ rj , for j = 1 to N f g . (4b)
The equilibrium constant of reaction r is given by [START_REF] Prigogine | Thermodynamique. Des moteurs thermiques aux structures dissipatives[END_REF]

K r = NI i=1 C νri Ai N f g j=1 C βrj Gj = exp -∆G r RT , (5)
where ∆G r is the Gibbs free energy of reaction r defined by

∆G r = ∆H r -T ∆S r . (6)
∆H r and ∆S r are, respectively, the enthalpy and entropy of reaction r, R is the universal gas constant and T the temperature. The determination of ionic and gaseous concentrations is carried out by setting the chemical affinity equal to zero for each of the R reactions due to the assumption of chemical equilibrium. Before presenting how this is achieved, the fining model will be described.

Fining model

Gas bubbles contained in the glass melt must be removed. As noted by Beerkens [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF], primary fining is the removal of bubbles through the action of fining agents releasing gases in order to enhance bubble growth. Since the rise velocity of a bubble is proportional to its radius squared, the time required for a bubble to reach the glass melt surface decreases with bubble growth. Primary fining is relevant at high temperature while secondary fining is the re-absorption of bubbles during cooling of the melt. As temperature decreases, the chemical solubility of species like O 2 , SO 2 increases. This effect leads to a reduction in bubble size.

Both bubble growth and shrinkage are controlled by mass transfer between bubble and melt. This process also takes into account multispecies diffusion. A detailed model can be found in reference [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF]. Here, we begin by outlining the mass transfer laws of gases diffusing into bubbles moving in a melt. Afterwards, the equations describing the dynamic variation of bubble position, radius and partial pressures will be specified.

Mass transfer

Assuming that bubbles are composed of N g gases (where N g = N f g + N nf g and N nf g corresponds to the number of non-fining gases such N 2 or H 2 O), the number of moles of gas species G j can be determined by

dn Gj dt = 4πa 2 k Gj (C Gj -C s Gj ), (7)
where n Gj is the number of moles of gas species G j , a the bubble radius, k Gj the mass transfer coefficient, C Gj the bulk molar concentration, C s Gj the molar concentration at the bubble surface. The mass transfer coefficient is derived from the Sherwood number.

Sh Gj = 2ak Gj D Gj . ( 8
)
D Gj is the diffusion coefficient of species G j . Mass transfer is also influenced by the velocity of the bubble relative to the molten glass. The importance of this motion is quantified by the Péclet number

P e Gj = 2au b D Gj . ( 9
)
u b is the bubble velocity relative to the molten glass. If we assume that the bubble interface is rigid, the Sherwood number must be proportional to the Peclet number to the power of one-third.

Here we use the law given in the textbook of Clift et al [START_REF] Clift | Bubbles, Drops, and Particles[END_REF]:

Sh Gj = 1 + (1 + P e Gj ) 1/3 . (10)
This equation agrees with the numerical solution within 2% over a large range of Péclet numbers. This is established for Stokes flow in an infinite domain however equation ( 10) can only be used for Re ≤ 1 (see for details [START_REF] Clift | Bubbles, Drops, and Particles[END_REF]) where Re is the Reynolds number given by

Re = 2au b ρ µ . ( 11
)
ρ is glass density and µ dynamical viscosity. The surface concentration, C s Gj , is determined using Henry's law:

C s Gj = L Gj P α G j Gj , (12)
where L Gj is the solubility coefficient, P Gj is the partial pressure of the species G j in the bubble and α Gj is a stoichiometric coefficient equal to unity for most species but is equal to 1/2 for water [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF].

Bubble growth and shrinkage

To describe the evolution of radius and partial pressure in a bubble, we use the ideal gas law:

P V = RT Ng j=1 n Gj , ( 13
)
where P is the total pressure in the bubble and V bubble volume. For each species G j , we have

P Gj V = n Gj RT. ( 14
)
The total pressure and partial pressures are linked by Dalton's law:

P = Ng j=1 P Gj , (15)
and the total pressure is defined by the following expression

P = P 0 + ρg(h -z) + 2σ a , (16)
where P 0 is the atmospheric pressure above the melt; h glass height; z the vertical position and σ the surface tension.

The equations describing the change of radius and partial pressure over time are obtained by the derivation, with respect to time, of the relations ( 13) and ( 14) and by inserting the expressions ( 7) and (16) (see [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF]). The rise velocity, growth/skrinkage rate and change in partial pressure are given by

dz dt = 2ρga 2 9µ , (17a)
da dt = 2ρ 2 g 2 a 3 27µ(P -2σ/3a) + RT P -2σ/3a Ng j=1 k Gj C Gj -L Gj P α G j Gj + aP 3(P -2σ/3a) 1 T dT dt - 1 P dP 0 dt , (17b)
dP Gj dt = 3RT a k Gj C Gj -L Gj P α G j Gj + P Gj T dT dt - 3P Gj a da dt , for j = 1 to N g -1. (17c)
In these equations, we have taken into account the possibility of varying atmospheric pressure and temperature.

Coupled models

In the two preceding subsections, the reduction-oxidation and fining processes have been presented independently. Due to the large quantity of bubbles, the consumption of dissolved gases changes the redox state of glass. Here, we present a coupled model similar to that previously presented by Klouzek and Nemec [START_REF] Kloužek | Modelling of glass refining kinetics. Part 2: Bubble distribution models and methods of measurement of refining properties[END_REF] and Matyas and Nemec [START_REF] Matyas | Behaviour of bubble files in glass melting space[END_REF]. The domain is assumed to be spatially homogeneous which means that spatial derivatives of all quantities are equal to zero. In other words, only time is taken into account as a parameter.

In order to describe both the reduction-oxidation state and the fining process, we must first determine the behavior of ionic and gas species in molten glass with respect to time. The temporal evolution of ionic and dissolved gases can be written as (see [START_REF] Kloužek | Modelling of glass refining kinetics. Part 2: Bubble distribution models and methods of measurement of refining properties[END_REF][START_REF] Matyas | Behaviour of bubble files in glass melting space[END_REF][START_REF] Prigogine | Thermodynamique. Des moteurs thermiques aux structures dissipatives[END_REF])

dC Ai dt = R r=1 ν ri dζ r dt , for i = 1 to N I . ( 18a
)
dC Gj dt = R r=1 β rj dζ r dt + S b,Gj , for j = 1 to N f g (fining gases), ( 18b
)
dC Gj dt = S b,Gj , for j = N f g + 1 to N g (non-fining gases). ( 18c
)
dζ r /dt is the rate of reaction r -given in mol/(m 3 •s) -and S b,Gj is the sink term due to the consumption of dissolved gas G j . If the gas G j goes from bubble to bulk, S b,Gj becomes a source term. S b,Gj is determined using the evolution of the number of moles n Gj given by equation [START_REF] Kloužek | Modelling of glass refining kinetics. Part 2: Bubble distribution models and methods of measurement of refining properties[END_REF]. If the number of bubbles with radius a is noted N b , S b,Gj is given by [START_REF] Kloužek | Modelling of glass refining kinetics. Part 2: Bubble distribution models and methods of measurement of refining properties[END_REF][START_REF] Matyas | Behaviour of bubble files in glass melting space[END_REF]

S b,Gj = -4πa 2 k Gj C Gj -L Gj P α G j Gj N b . ( 19
)
To determine the evolution of C Ai and C Gj , we must know the rate of reaction, dζ r /dt. In order to find these values, we use the equilibrium constants for each reaction r defined by the equation [START_REF] Doremus | Glass Science[END_REF]. By taking the logarithmic derivate of (5), we have :

NI i=1 ν ri C Ai dC Ai dt + N f g j=1 β rj C Gj dC Gj dt = d ln K r dT dT dt , for r = 1 to R. (20)
Finally, using (18a) and (18b), we obtain a linear system of R equations:

R k=1 M rk dζ k dt = d ln K r dT dT dt - N f g j=1 β rj C Gj S b,Gj , for r = 1 to R. ( 21
)
The coefficients of the symmetric matrix M, M rk , are given by

M rk = NI i=1 ν ri ν ki C Ai + N f g j=1 β rj β kj C Gj . (22)
Thus, the rates of reaction are not equal to zero if the temperature varies or if there is a consumption of the dissolved gases is present. At this stage, it is possible to numerically couple the redox behavior and fining process. The method is based on the numerical solution of the two systems (17) and ( 18) with (21) and the integration of the differential systems is carried out using the fourth order Runge-Kutta method.

Results

Fining evolution with coupled model

We consider a glass where the following two reduction-oxidation reactions take place

F e 3+ + O 2- 2 ⇋ F e 2+ + O 2 4 , (23a)
SO 2- 4 ⇋ SO 2 + O 2 2 + O 2-. (23b)
The enthalpy ∆H r and entropy ∆S r for these reactions are given in Table 3. These values come from reference [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF]. The assumed concentrations of ionic and gaseous species at T = 1200 • C are given in Table 4.

The coupled model has been tested on bubble behavior in a soda-lime-silica glass melt. The height of the glass bath is taken to be 1 meter and the atmospheric pressure equal to 1 bar. The bubbles are composed of five gas species, two fining gases: O 2 and SO 2 and three non-fining gases: CO 2 , N 2 and H 2 O. Initially, we assume that the bubble contains only CO 2 and that the radius is equal to 100 µm with the solubilities and diffusion coefficients for these elements taken from the data of Beerkens [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF]. The concentrations of CO 2 , N 2 and H 2 O in the glass are given in Table 5. We carried out two computations at T = 1400 • C: the first one with N b = 0 and the second with N b = 10 8 bubbles/m 3 .

In Figure 1, we present the evolution of bubble radius as a function of the time for the two cases. The calculation is stopped when the bubble reaches the free surface. As can be seen, due to consumption of dissolved gases, the growth rate of the bubble decreases with an increasing concentration of bubbles. Indeed, the quantity of dissolved gas decreases with time, therefore, the gradient between the bulk and the bubble also decreases. We see also that the time needed by the bubble to reach the upper boundary is greater when N b = 10 8 bubbles/m 3 . For N b 0, this time is equal to 2.69 h and for N b = 10 8 bubbles/m 3 , it is equal to 3.21 h (i.e. a 20% of increase).

The influence of bubble concentration can also be observed on the evolution of molar fractions in the bubble presented in Figure 2 for the five gas species. Table 6 and Table 7 give the final molar fraction in each of the two cases. The molar fraction of O 2 in the bubble at the end of evolution is more important when N b = 10 8 bubbles/m 3 .

Equilibrium of gases in glass by bubbling

In this last subsection, we used the coupled method in order to find the time needed in order to equilibrate the glass with bubbles at certain initial concentrations. This is important if we wish to experimentally determine the equilibrium constants of reduction-oxidation reactions. The question is : how much time is necessary to reach a state of equilibrium? It is possible to use the coupled model presented in this article to realize this "numerical" experiment. We inject bubbles at a fixed gas flow rate into a crucible filled with one kilogram of glass. In order to reproduce bubbling flow, when a given bubble reaches the surface of the melt it is re-injected through the crucible bottom with its initial concentrations. The size and the concentration of bubbles are determined using bubbling laws previously detailed by Chavanne et al [START_REF] Chavanne | Bubbling laws in a glass furnace[END_REF]. Computational tests have been performed at four temperatures: 1200, 1300, 1400 and 1500 • C. The radius and concentration of the bubbles are given in Table 8. The gas flow rate is fixed at 2 l/min under standard conditions (T = 25 • C and P = 101325 Pa) and we use the physical properties of a flat glass composition with the solubilities and diffusion coefficients being taken from Beerkens [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF]. We assume air bubbles to mean that the molar fraction of O 2 is fixed at 0.21 and N 2 at 0.79. It should be specified that only the gases O 2 , SO 2 and N 2 are taken into account in this computation.

Figure 3 presents the evolution of O 2 activity in the glass as a function of time for the four temperatures. The evolution of P O2 has three stages: in the first, the activity does not change; after that, we observe increasing P O2 followed by a decrease and in the last stage, the P O2 tends to the partial pressure fixed by the bubble. The equilibrium time reduces as the temperature increase. All these facts are physically intuitive, nevertheless, the equilibrium times found in this calculation are too large. Indeed, at T = 1500 • C, we find that the time to obtain equilibrium is approximatively equal to 16.7 days! A value which reaches 83 days at T = 1400 • C. In the "real" laboratory experiments, the order of magnitude is a few hours.

In order to improve these results, we have modified the Sherwood law. In fact, due to the large bubbles, the Reynolds number can be greater than one. With the sizes given in Table 8, Re is in the range of 9 to 120. We therfore use a different Sherwood law given by Clift et al [START_REF] Clift | Bubbles, Drops, and Particles[END_REF] :

Sh = 0.991P e 1/3 1 + Re 4 0.27 . ( 24
)
This equation is an extension of Levich's law [START_REF] Levich | Physicochemical hydrodynamics[END_REF] for Reynolds numbers greater than one. However, even with this modification, the equilibrium times are not significantly changed. At 1600 • C, for example, equilibrium time is only reduced to 16 days, which is the same order of magnitude. As presented by Beerkens and Waal [START_REF] Beerkens | Mechanism of oxygen diffusion in glassmelts containing variable-valence ions[END_REF] and by Beerkens [START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF], we have used an effective diffusion coefficient of O 2 in order to take into account the reaction of oxygen in the melt. This effective diffusion is evaluated using the same method proposed by Beerkens and Waal [START_REF] Beerkens | Mechanism of oxygen diffusion in glassmelts containing variable-valence ions[END_REF]. We find that the effective diffusivity is 6 times smaller than diffusion coefficient of O 2 on its own. However, the Sherwood number is proportional to the Péclet number to the power of one-third which means that decreasing the diffusivity by a factor 6 leeds to only a factor 1.8 on the Sherwood number. Consequently, the disagreement between numerical results and experimental data presented here cannot explained by this effect. We need to take into account the effect of sulphate reaction and in that way couple the O 2 and SO 2 diffusivities.

Moreover, to describe bubble ascension, we have used Stokes's law. However, the size of the bubble is not negligible compared to the crucible and therefore hydrodynamic interactions might impose a drag force on the bubbles (see reference [START_REF] Feuillebois | Some theorical results for the motion of solid spherical particles in a viscous fluid[END_REF]). Bubble velocity can thus, be decreased and this will affect also mass transfer. This improvement is outside of the scope of this article and will be addressed in the future.

The disagreement between model and experiment could also be due to the fact that, in the present method, the domain is spatially homogeneous and the mass transfer might be in fact controlled to local phenomena. This kind of description must take into account local advection and diffusion processes.

Conclusion

The objective of this investigation was the modelling of fining and redox processes. Due to the large quantity of bubbles present in glass during the first stage of melting, the consumption of dissolved gases can modify the redox state. In order to take this into consideration, we have presented a coupled model of redox reactions and fining processes. The method is based on the sequential solution of sets of differential equations, describing redox reaction and fining evolution.

The present model has been applied to observe the influence of the concentration of small bubbles present in the glass melt on growth rate. We have shown that a large density of bubbles changes the bubble growth rate and that the rise time increases when bubble population is large. Moreover, the average composition of the bubble also changes with bubble concentration.

Finally, we have tried to apply our model to describe the gas concentration in a glass melt due to bubbling. The model predicts the transient evolution of redox state which tends to reach an equilibrium between the two phases; however, the equilibrium time is very large compared to experimental reality. The calculations have been also performed taking into account the effect of redox reactions on the diffusion process of O 2 by changing the diffusion coefficient. We have also modified the mass transfer coefficient due to the vigorous convection of bubbling flow. Both of these modifications do not, however, significantly improve the results.

In order to improve the coupled model in the future, it is essential to take convection into account. This insertion of our model into a CFD software might be one way to achieve such a goal. Reac.

Nomenclature

∆H r (J.mol -1 ) ∆S r (J.mol -1 .K

8 TFigure 1 :Figure 2 :

 812 Figure 1: Increase of bubble radius as a function of time with N b = 0 and N b = 10 8 bubbles/m 3 .

Figure 3 :

 3 Figure 3: Evolution of O 2 activity of float glass melt as a function of time for T = 1200, 1300, 1400 and 1500 • C during bubbling with air.

Table 3 :

 3 ∆H r et ∆S r of iron and sulfate reactions from[START_REF] Beerkens | Analysis of advanced and fast fining processes for glass melts[END_REF].Species C Ai or C Gj (mol/m 3 )

	-1 )

Table 4 :

 4 Ionic and gas concentrations at T = 1200 • C.

	Species C Gj (mol/m 3 )
	CO 2	0.682
	N 2	2.1 • 10 -2
	H 2 O	38

Table 5 :

 5 Non-fining gas concentrations in the glass.

	Gas	SO 2	O 2	CO 2	N 2	H 2 O
	Molar fraction 0.28 0.10 0.47 0.08 0.06

Table 6 :

 6 Molar fraction of gas species in the bubble at the end of evolution with N b = 0.

	Gas	SO 2	O 2	CO 2	N 2	H 2 O
	Molar fraction 0.31 0.15 0.45 0.02 0.06

Table 7 :

 7 Molar fraction of gas species in the bubble at the end of evolution with N b = 10 8 . N b (m -3 ) 3.50 • 10 3 7.07 • 10 3 1.22 • 10 4 1.82 • 10 4

	T ( • C)	1200	1300	1400	1500
	a (mm)	30.5	24.3	20.3	18.0

Table 8 :

 8 Calculated a and N b at T = 1200, 1300, 1400 and 1500 • C for a gas flow rate equal to 2 l/min under standard conditions.


Ph.D Kévin Rousseeuw
email: kevin.rousseeuw@gmail.com

Member, IEEE Émilie Student

Alain Poisson-Caillault
email: emilie.caillault@lisic.univ-littoral.fr

Denis Lefebvre
email: alain.lefebvre@ifremer.fr

Hamad

PREPRINTS JSTARS 2014

Keywords: Hybrid Hidden Markov Model, marine water monitoring, Phytoplankton blooms, spectral clustering. I

Phytoplankton is an important indicator of water quality assessment. To understand phytoplankton dynamics, many fixed buoys and ferry boxes were implemented, resulting in the generation of substantial data signals. Collected data are used as inputs of an effective monitoring system. The system, based on unsupervised Hidden Markov Model (HMM), is designed not only to detect phytoplancton blooms but also to understand their dynamics. HMM parameters are usually estimated by an iterative Expectation-Maximisation approach. We propose to estimate HMM parameters by using spectral clustering algorithm. The monitoring system is assessed on database signals from MAREL-Carnot station (Boulogne-sur-Mer, France). Experiment results show that the proposed system is efficient to detect environmental states such as phytoplankton productive and non productive periods without a priori knowledge. Furthermore, discovered states are consistent with biological interpretation.

forecasting issue using a hybrid HMM. The specific objective of this work is to design a system able to model phytoplankton dynamics from large database and no prior knowledge. For this purpose, a fully unsupervised HMM is built using spectral clustering algorithm to generate HMM symbols and states.

The paper is organized as follows. Section II describes the monitoring system and the proposed hybrid HMM with three parts. Part II-A discusses about usual unsupervised techniques to build HMM, and spectral clustering approach is argued to estimate HMM static parameters (states and symbols). Part II-B details HMM symbol generation by a self tuning fast K-means proposed algorithm. Part II-C defines HMM state generation by spectral clustering algorithm. Section III describes protocol of experimentations and collected data used from the IFREMER MAREL-Carnot station that registers water characteristics at HF resolution. First a fixed 2-state HMM is built in order to assess symbol and state generations thanks to an artificial labeling. Thus, our algorithms are compared with other machine learning techniques. Then, in section IV experiment results of a fully unsupervised N-state HMM are presented, and are related to examine biological expectations. II. MONITORING SYSTEM BASED ON HYBRID HMM Fig. 1 presents the proposed monitoring system architecture. Data collected at high frequency resolution from 2005 to 2008 are first pretreated. Then, the clustering step is applied in order to find environmental states. The final step relies on temporal information between these states to develop a phytoplakton dynamics model. The built model is used to predict a new or forthcoming phytoplankton bloom, or specific states (classification/alert box in Fig. 1).

A. Hidden Markov Model

According to normal course of phytoplankton succession highlighted by Margalef [START_REF] Margalef | Life-forms of phytoplankton as survival alternatives in an unstable environment[END_REF] and Reynolds et al. [START_REF] Reynolds | Towards a functional classification of the freshwater phytoplankton[END_REF] works, we assume that phytoplankton biomass is constrained by a high level of dependence among successive observations. Besides, it may be viewed as the result of a probabilistic walk along the environmental states. So let us see how to design one ergodic Hidden Markov Model to characterize the dynamics of phytoplankton blooms from physico-chemical and biological parameters in an unsupervised context.

A HMM noted λ = λ(S, V, π, A, B) is defined with 2 static sets (S, V) and 3 computed sets of probabilities (π, A, B) that we recalled as follows [START_REF] Rabiner | A tutorial on Hidden Markov Models and selected applications in speech recognition[END_REF]:

1) S = {s 1 , s 2 , . . . , s N } is the set of states with N the number of distinct states. For instance: non-productive period, pre-bloom, bloom, post-bloom and other rare events, like dam opening, factory or agricultural discharge. The number of states is generally set by expert people in relation with the applications, or automatically determined by penalized maximum likelihood criterion. 2) V = {v 1 , v 2 , . . . v M } is the set of symbols and M is the number of distinct symbols. In the simpliest cases, observation symbols correspond to the system outputs. Environmental state is not characterized by a unique representative per state, and two system outputs can belong to several states. Thus, it is necessary to build a codebook of symbols by vector quantization [START_REF] Koo | An improved vq codebook design algorithm for hmm[END_REF]- [START_REF] Debyeche | Improved vector quantization approach for discrete hmm speech recognition system[END_REF]. So data space will be represented by this codebook V of M symbols. 3) π = {π i } of size N , defines the initial probability distribution, π i = P (s(t = 1) = s i )}. There is no information about the state that will be predominant during data acquisition. A priori initial states are equiprobable. 4) A = {a ij } of size N ×N , defines the transition matrix with a ij = P (s(t) = s i |s(t-1) = s j ) the conditional probability.

Therefore, the number of times that we move from a state s i to state s j is estimated, then A is normalized in row. 5) B = {b ik } of size N × M , defines the emission probability with b ik = P (v(t) = v k |s(t) = s i ). From a finite observation sequence without any labeled states, HMM symbols, transition and emission matrices [START_REF] Liao | Identification of ground targets from sequential high-range-resolution radar signatures[END_REF] are adapted iteratively. Expectation-Maximisation approach (EM) is used with entropy criterion with Minimum Description Length When new data o t is collected, it is associated with its nearest symbol. Viterbi algorithm [START_REF] Viterbi | Error bounds for convolutional codes and an asymptotically optimum decoding algorithm[END_REF], [START_REF] Forney | The viterbi algorithm[END_REF] is then applied to estimate its environmental state.

B. Symbol generation

MAREL-Carnot database consists of 26,280 × 19 parameters per year as from November 2004. To discover underlying states in this large database, instance selection is required. K-means algorithm is a well-adapted vector quantization method, and is popular for data clustering too [START_REF] Jain | Data clustering: 50 years beyond k-means[END_REF]. The main idea is to build vector prototypes from a set of observations denoted O = {o 1 , o 2 , . . . , o Np } of N p data points preserving HF information. The fast K-means algorithm [START_REF] Shindler | Fast and accurate k-means for large datasets[END_REF] is modified to obtain a self tuning K-means based on Hartigan-Wong algorithm [START_REF] Hatrigan | A k-means clustering algorithm[END_REF], and on Elbow criterion: the number of clusters K is incremented until a fixed percentage of explained variance or a Kmax number of retained prototypes (i.e. symbols) is met. The principle of this proposed algorithm, named Self Tuning Fast K-Means (STFKM), is described in Fig 3 . The center initialization for large databases (number of points N p > 20, 000) is very important here to speed up the process convergence. Kmax is the maximum number of reduced points specified by the user or in the default case, the number of measures in the time series. varExplained is the explained variance desired by the user, by default this number is set to 95 percent.

k = k + 1; 10:
Step 1: Initialization of k centers 11:

..Cut Data in n subsamples of 20,000 points 12:

..Compute K-means (K=k) on each subsample 13:

..Select the k clusters centers from the best partition according M S(within)/M S(between) 14:

Step 2: Decide the class memberships of the N p points by assigning them to the nearest center.

15:

Step 3: Re-estimate the k cluster centers, by assuming the new memberships found 16:

Step 4: If none of N p points changed membership in the last iteration, Otherwise goto 2. 17:

Step 5: vE = M S(between)/M S(total) return k obtained centers 20: end procedure Fig. 3. Outline of Self Tuning Fast K-Means algorithm noted STFKM. The variance, MS(.) for Mean Square defined the variance between or within groups, or the total.

C. State generation by spectral clustering

After STFKM procedure on the pretreated data, M symbols summarize the entire database. From these M symbols, N states are detected by unsupervised clustering. Each MAREL-Carnot physico-chemical parameter follows a stochatic, nonlinear and non-stationary process (except sea-level), see section III. They have not-gaussian distributions, and environmental state characterisation is unknown. So, SC method is the best way to avoid some assumptions about data shape. SC is capable of classifying data which are connected, but which are not necessarily compact, or clustered within convex boundaries. Indeed, the key idea of SC is to transform the input data space into a new feature space where K-means clustering could be applied. The most typical method by Ng, Jordan et al. [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] is recalled in Fig. 4. Lap = D -1/2 WD -1/2 ; Laplacian matrix 6:

Select the K-largest eigenvectors x of Lap; Form the Y matrix from the row-normalisation of X thus

y ij = w ij /( j x 2 ij ) 1/2 9: l =K-means(Y,K) each row of Y is a point 10:
Assign original Point o i to the cluster l i 11:

return label vector : l 12: end procedure The number of clusters K in input of the SC algorithm and the way to build the Gram Affinity matrix W have both significant effects on the classification result. Gaussian kernel function is the most widely used function for constructing W = {w ij } defined as:

w ij = exp(- o i -o j 2 2σ 2 ) (1)
The scaling parameter σ helps to sparse the matrix and tends to obtain an ideal case with a robust eigen-decomposition (i.e. in the ideal case, the first K eigenvalues are equal to one). However, a bad choice of σ brings an incorrect classification. Zelnik and Perona (ZP) [START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF], or Kong et al. [START_REF] Kong | Robust and smart spectral culstering from normalized cut[END_REF] proposed a local scale parameter σ i for each data o i based on its neighborhood, instead of selecting a uniform scaling parameter σ. The ZP affinity matrix W is chosen with a z-neighborhood (o nz the z th neigborhood of the point o i ):

w ij = exp(- o i -o j 2 2σ i σ j ) with σ i = o i -o nz (2)
Many authors proposed to overcome the choice of the number of clusters K by analysing either eigenvalues magnitude (equal or nearest to one) or eigengap, or eigenvectors (see [START_REF] Kong | Robust and smart spectral culstering from normalized cut[END_REF], [START_REF] Sanguinetti | Automatic determination of the number of clusters using spectral algorithms[END_REF], [START_REF] Xiang | Spectral clustering with eigenvector selection[END_REF]). To select the number of states N for HMM topology, the eigengap method is used: it is the simpliest one to implement, and it has the lowest complexity.

From the spectral clustering of the M symbols V = {v 1 , v 2 , . . . v M } issued from the STFKM step, we assign the observation data O = {o 1 , o 2 , . . . , o Np } thanks to the label s i = l k of their cluster center v k .

III. DATA AND FIXED 2-STATE HMM VALIDATION

One HMM is built according to the scheme described in the previous section from MAREL-Carnot multivariate marine signals in order to model the phytoplankton dynamics in the French Channel coast around Boulogne-sur-Mer. Data and their curves are available on the website (http://www.ifremer.fr/difMarelCarnot/) with authorisation request. These data are first presented, then the experiment validation protocol follows.

Without ground truth on the environmental states and in order to assess our system, we decide to create an automatic data labeling based on the monitoring sampling strategy for the EU Water Directive Framework (EU-WFD). Thus, data from March to October are labeled s 1 , corresponding to the productive period (in terms of biomass production capacity), and the others s 2 for the non-productive period. Furthermore, this labeling will allow to compare our sytem with other machine learning algorithms. I. Note that Fluorescence signal is not taken into account, but is used to validate clustering results since we have no ground truth: it is a presence indicator of phytoplankton cells.

A. Data presentation

When an observation data contains at least one missing value (among in R 10 ), this point does not participate to the generation of symbols and states. Centering and standard deviation scalings are achieved on each parameter to avoid parameter range influence.

Data from 2005 to 2008 are considered to build HMM parameters. To test the time modeling, data from 2009 will be tested with the same pretreatment protocol.

B. Vector quantization validation

The number of symbols V required to characterize a state is first analyzed. Selection of the M symbols from the data is performed by 100 random drawings. A 1-Nearest Neighbor algorithm (1-NN) is used to evaluate the approximation of a mixture of components per state. 1-NN is first applied on each parameter, and then on the multidimensional matrix (from 2005 to 2008). K-means algorithm is also applied to build K = M representatives per state.

Two scores are considered: Rate Recognition (RR) and the monthly Overlap defined by the following equation: |.| is the cardinal operator and |s 1 (i)| defines the number of points labeled s 1 during the i th month. Phytoplankton productive and non-productive periods are expected to have no overlap according to EU-WFD.

Overlap = i [|s 1 (i)| + |s 2 (i)| -max(|s 1 (i)|, |s 2 (i)|)] |s 1 ∪ s 2 | (3)
For the monodimensional analysis, water temperature ETCO1 is the most discriminative parameter, with a recognition rate from 75.1% (± 3.5) for one representative per state, to 77.8% (± 0.4) for 1,000 representatives. For the multidimensional analysis, Table II summarizes the mean and standard deviation of the two scores, RR and Overlap, for different M-values. Approximating data distribution with one unique random representative gives poor recognition rate (68.1%) and often an important Overlap (18.4%) of the two desired environmental states. To decrease this Overlap around 10%, more than 100 random representatives are required.

K-means is a geometric approach adapted for linearly separable data sets. This algorithm requires to know the desired number of symbols (centers) M . Here with 10 symbols per state, the Overlap is around 10%. To reduce this Overlap around 5%, 100 representatives per state are required.

The proposed STFKM automatically searches the number of symbols that describes the data structure, and it is fast running. The impact of the STFKM selection is tested with a learning machine: a Support Vector Machine (SVM). A SVM model (radial basis kernel) was trained on 2005-2008 data and was tested on 2009 data with 10 cross-validation. Three experiments for the SVM training were led: on all training data, on 1,000 randomly selected representatives per state (M =2,000) of this data, and on symbols issued from STFKM vector quantization.

Table III summarizes SVM capacities of training and generalisation (test) for these 3 studies with RR and Overlap scores. In generalisation, STFKM algorithm allows to keep a similar recognition rate (92.6%) and Overlap score (7.4%) to no sampling-SVM. STFKM-SVM gives better training capacity than a random selection of M symbols. In this supervised context, we can conclude that the obtained vector quantization by STFKM is a relevant data reduction.

The stability of STFKM algorithm is then assessed according to the Rand Index (RI) of 10 achieved symbol generations. RI score [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF] is a measure of similarity between two data clusterings P 1 and P 2 of a given set of n elements E = {e 1 , . . . , e n }. Note that the number of clusters in each partition can be different. RI is computed according to the following equation:

RI(P 1 , P 2 ) = (a + b) n 2 (4)
where a (resp. b) is defined by the number of pairs of elements in E that are in the same set in P 1 (resp. in different sets) and in the same set in P 2 (resp. in different sets) .

In a fully unsupervised context, STFKM approach allows to respect the high frequency information without losing data structure. Table IV shows that RI scores of the 10 obtained partitions and the number of retained symbols are quite similar. The RI score near to one shows that STFKM algorithm gives a robust vector quantization.

C. State generation validation

Then, the spectral clustering algorithm is compared to usual unsupervised approaches : EM, Hierarchical Clustering (HC) for a set number of states N = 2 with the same STFKM symbols. For experiments, a 7-neighborhood is considered for the scaling parameter of the similarity matrixin SC. We used EM and HC algorithms implemented in R-Gui (library Mclust and stats: http://www.r-project.org/). Only the best or runnable options are retained: Expectation-Maximisation (EII model), HC with complete linkage. EM is also performed on each of the 10 parameters, only results of the most discriminative parameter (salinity CSAL1) are presented. 1-NN algorithm is used to label data from the built model. The 2005-2008 period is named building database, and the year 2009 test database.

Table V presents RR and Overlap scores to analyse jointly. In spite of its very low Overlap score around 0.5%, hierarchical clustering (cutting obtained tree to 2 clusters) does not separate productive and non-productive periods. State 1 represents 84,118 of 84,614 points, 99% of the building database and state 2 concerns few points in August, September and November for the building database and February, April and June for the test database without any biological or sensor interpretation. EM approach offers the best RR results for the two databases. STFKM approach gives lower RR than EM one, but its Overlap for the largest database (building database) is reduced. STFKM-SC approach is a balanced one for EU-WFD labeling, and will be relevant for a number of states greater than 2. Indeed, we expect that our hybrid HMM system can detect more than 2 states like phytoplankton spring bloom or automnal bloom, rare events.

D. Time modeling validation, fixed 2-state HMM

We evaluate, through experiments, the reliability of our hybrid model: the entire procedure for building one HMM from clustering is repeated 10 times. The 10 partitions have a mean RI score of 0.95, so the whole STFKM-SC step (symbol and state generation) is robust. We keep the symbol and state partition with the smallest normalised multi-cut, MNCut [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] to build HMM. For experiments, the number of states is set: N = 2 and other parameters are automatically tuned: explained variance is fixed to 95 percent. According to the MNCut criterion, HMM built has M = 2794 symbols.

According to the EU-WFD labeling, 79.3 percent of building database is well recognized with an Overlap of 11.7%, and 82.1 percent of 2009 test database is well recognized with an Overlap of 6.7%. Fig. 6 illustrates the distribution of labeled states by HMM prediction for the building database and the 2009 test database. In 2005-2008, state s 1 in red color ties in with the period from March to December with a dominant April-November period whereas state s 2 in green color is dominant in the November-April period. Over the period 2009 state s 2 ties in with the period from April to October, whereas state s 2 in green color is dominant in the December-April period. Many data in March and August-November have no estimated state (noted NA in black color in Fig. 6) due to one least missing value in R 10 ; that means the system forecasts confusions in transition periods. But states s 1 and s 2 match well with the two main environmental EU-WFD states: s 1 and s 2 characterize the phytoplankton dynamics, the productive and the non-productive periods.

To validate the relevancy of the built model from 2005-2008 and its symbol quantization in a fully unsupervised context, MAREL-Carnot signals Obs on the year 2009 are reconstructed and compared with the original data Obs. For an observation obs t , the system estimates its state s i with the higher likelihood. Then the most present symbols v k in this state are retained, see Fig. 7. A similarity sim score is defined by the following equation 5: From 2005 to 2008, each year participates to HMM building, and reconstructed signals are compared to the original ones so as to assess the modeling power. Table VI shows the similarity scores, Eq.( 5), between original and reconstructed signals for each year over the period 2005-2008 and their related mean. The last column corresponds to the year 2009. This similarity score is bounded by 0 and 1: 1-value implies that reconstructed signals are exactly the same as the original ones. The similarity scores and their mean are greater than 0.83 from 2005 to 2008. So, we conclude that reconstruted signals are very close to the original data, and that the vector quantization algorithm for HMM states is efficient. Therefore, the proposed system has an interesting generalisation power, for the year 2009, which does not participate in HMM building. Indeed, the time series are built with a high similarity score: 0.79 according to the choice of the most probable symbol of the state. Fig. 8 illustrates original and reconstructed signal of one parameter, dissolved oxygen concentration C O21 for 2009.

sim(Obs, Obs) = 1 |Obs| |0bs| t=1 1 obs t -obs t + 1 (5)

IV. N-STATE HMM FOR PHYTOPLANKTON DYNAMICS

Considering that HMM is now validated on a fixed 2-state biological dynamics, the next step is to increase the number of states to refine and to try to better understand the bloom determinism and its dynamics. To interpret this model and to relate a ecological meaning, we analyze the state sequencing and the characterization of each state by a correlation analysis (Table VII).

State s 6 (yellow) clearly highlights high salinity values, ranging from 33.9 to 36.2 with a mean value of 35.4 (Fig. 5). These values are more representative of offshore waters. And then we can conclude in this coastal zone that mainly state s 6 corresponds to salinity anomalies (sensor failures). Nevertheless, s 6 may sometimes be explained when west winds persist and consequently bring more offshore waters to the coast. State s 3 (blue) is representative of the winter non-productive period, with high nutrient concentrations and low temperature (Table VII). The initiation of the main phytoplankton bloom, and then the growing phytoplankton stage (between February and May: interannual variability of the bloom) are characterized by s 1 (red). High oxygen concentrations, explained by a high phytoplankton production (photosynthesis) are observed during this stage (Table VII). During state s 1 , phytoplankton mainly uses the winter nutrients stock, and consequently this state corresponds to the new production period [START_REF] Gentilhomme | Seasonal cycle of nitrogen and phytoplankton biomass in a well-mixed coastal system (eastern english channel)[END_REF]. States s 2 and s 4 follow state s 1 , and are identified as the regenerated production period when phytoplankton production is based on regenerated nutrients (transformation of the organic matter from the previous bloom -state s 1 -into new available nutrients).

States s 5 (pink) and s 7 (grey) correspond to rare or short events, respectively, with high turbidity during storm events, and high phosphate and silicate concentrations (C PO1 and C SI1 in Table VII). More investigations are required to better understand the main processes involved during these periods.

Fig. 11 illustrates the predicted states for the year 2009 with their sequencing. The 7-state HMM succeeds in predicting phytoplankton biomass dynamics. The state sequencing matches our assumption with a pre-bloom winter period (state s 3 mainly) followed by the main phytoplankton bloom based on external nutrient inputs (state s 1 ), and the regenerated bloom (states s 2 and s 4 ).

V. CONCLUSIONS AND FUTURE WORKS

Two N-state HMM was built in order to forecast phytoplankton blooms near the French Channel coast from MAREL-Carnot signals (IFREMER, Boulogne-sur-Mer) without any biological knowledge. HMM building requires at least two parameters: a number of states, a number of symbols that characterize states. These parameters are commonly estimated iteratively by Expectation-Maximisation. We propose a one-pass process to estimate HMM symbols and states in a fully unsupervised context. A proposed Self Tuning Fast K-Means STFKM algorithm extracts symbols from observation data. From this vector quantization, a spectral clustering approach, with no tuning too, generates HMM states that allows to treat non-convex data. A signal reconstruction approach is proposed to assess HMM prediction.

Result analyses from the MAREL-Carnot buoy data first demonstrate interests and the stability of each used algorithm (state and symbol generation) throughout the monitoring chain. The high resolution information is preserved. Built 2-state HMM permits to detect the main productive and non-productive periods, as used for the purposes of the EU Water Framework Directive to assess good environmental status. A 7-state HMM was proposed to refine knowledge about phytoplankton bloom dynamics in a temperate ecosystem, temporarily dominated by a harmful algae (Phaeocsytis globosa). The obtained state sequencing coincides with dynamics described using measurements from low resolution system near the MAREL-Carnot (Rephy/SRN data [START_REF] Lefebvre | Long-term hydrological and phytoplankton monitoring (1992-2007) of three potentially eutrophic systems in the eastern English Channel and the Southern Bight of the North Sea[END_REF]). The proposed HMM system succeeds in characterizing phytoplankton dynamics from new incoming data (in near real-time approach). Using the main statistical characteristics of the parameters underlying the definition of a given state, the system will allow to further increase knowledge about the main controlling or forcing parameters (i.e., nutrient pressure, light availability, turbidity), the environmental status (e.g., phytoplankton biomass), and the direct and/or indirect effects of such blooms (e.g., oxygen concentration). The main limiting step in the monitoring chain is the removing samples with missing values. Indeed, the latter affects the state estimation and characterization (symbol process). Some phytoplankton blooms were not taken into account for HMM building.

Several environmental monitoring and research programmes could benefit from the proposed method to avoid the critical expert labeling step when modelling. It could help to process large multivariate time series as generated by high resolution (in time and/or space) platforms, more and more frequently implemented for the integrated observation of pelagic ecosystems and biogeochemical cycles in the oceans. Moreover, the possibility of identifying environmental states (characterized by a combination of several parameters) is a clear opportunity to better understand what a good environmental status is, as defined and used for the needs of the WFD, the MSFD or other regional sea convention (as OSPAR).

Fig. 1 .

 1 Fig. 1. General scheme of the monitoring system.

Fig. 2 .

 2 Fig. 2. Hybrid HMM building scheme. Dotted line separates HMM modelling on left side from state prediction of a new observation obst on right side.

1 : 3 :

 13 procedure SPECTRALCLUSTERING(O, K) 2: Variable: W, D, Lap, X, Y, l Compute a Gram affinity matrix W M ×M from O

7 :

 7 Form the matrix X = [x 1 x 2 . . . x k ] ∈ R M by stacking the eigenvectors in column 8:

Fig. 4 .

 4 Fig. 4. Spectral clustering algorithm

  MAREL-Carnot station registers 19 signals: 16 water characteristics every 20 minutes, and 3 nutrient levels every 12 hours. These signals are detailed in TableI. Collected data signals come from different sensors. They require pretreatments to respect sensor range, and sensor time alignment. In case of sensor failures, its measurements are not retained (all pH values and measured sea-level values are removed). The sensor ranges are adapted to Boulogne-sur-Mer ecosystem. Time alignment is obtained by a moving average technique on a small sampling rate (20 minutes). Nutrient parameters are duplicated to obtain the same time resolution. From this step, signal database for 2005-2008 contains 105,192 points in R19 . Only half of the data (48,157 points) have no missing values (due to sensor default). To reduce missing value, a moving average over one week (temporal scale of a bloom) is applied, leading to a completed data of N p =84,614 points.

Fig. 5 .

 5 Fig. 5. Five completed MAREL-Carnot signals from 2005 to 2009 by mobile average: C SAL1, E TU1, C NI1, C PO1, C SI1. Vertical line separates data into 2 subsets. Subset O corresponds to 2005-2008 period which is used for HMM building. Subset Obs relates to the year 2009 which is used for generalisation and accordingly does not participate in HMM building.

Fig. 5

 5 Fig. 5 illustrates 5 signals after this data completion: Salinity, Turbidity and nutrients (Nitrate, Phosphate and Silicate); residual missing measurements of nitrate concentrations, for instance, can be visualized. After a correlation analysis, N f =10 physico-chemical parameters, not correlated, are retained and detailed in NC column of the TableI. Note that Fluorescence signal is not taken into account, but is used to validate clustering results since we have no ground truth: it is a presence indicator of phytoplankton cells.When an observation data contains at least one missing value (among in R 10 ), this point does not participate to the generation of symbols and states. Centering and standard deviation scalings are achieved on each parameter to avoid parameter range influence.Data from 2005 to 2008 are considered to build HMM parameters. To test the time modeling, data from 2009 will be tested with the same pretreatment protocol.

  (a) Period 2005-2008, HMM building database. (b) 2009, test database.

Fig. 6 .

 6 Fig. 6. State distribution per month for a typical seasonal cycle: state s 1 is represented in red color and s 2 in green color. The black color, named NA, concerns measures whose state is not estimated (one least missing values).

Fig. 7 .

 7 Fig. 7. Time series reconstruction from HMM classification. Signal Obs is reconstructed from original observation Obs thanks to the M symbols.

  A 7-state HMM with M =2,884 symbols is built from 2005-2008 database, N = 7 according to the eigengap technique. Time modeling validation is achieved according to the same protocol as the one of the 2-state HMM. Reconstructed 2009 signals have a similarity score above 0.8 with the original data. Fig. 9 is the color-state projection of the estimated states over the period 2005-2008 on Fluorescence signal; black color denotes unlabelled observations due to missing value. The state distribution and sequencing are illustrated in Fig. 10 with the same color standard.

  (a) C O21 original signal. (b) C O21 estimated signal.

Fig. 8 .

 8 Fig. 8. 2009 -Dissolved oxygen concentration signals: original signal (a) and estimated signal (b) by HMM. In gray color, time measurements are not estimated by HMM due to at least one missing parameter at this time.

Fig. 9 .

 9 Fig. 9. Clustering results: Color-state sequencing projection on the fluorescence signal for 2005-2008. Black color, named NA concerns measures whose state is not estimated (one least missing value).

  (a) State distribution per month for a typical seasonal cycle over the period 2005-2008. (b) Measure sequencing in each state over the period 2005-2008.

Fig. 10 .

 10 Fig. 10. Clustering results: NA concerns measures whose state is not estimated (one least missing value).

  (a) Color-state projection on Fluorescence signal. (b) State sequencing.

Fig. 11 .

 11 Fig. 11. Predicted state results in 2009 by HMM: Predict state results in 2009 by HMM: (a) Color-state projection on Fluorescence signal and (b) state sequencing. NA (black color) corresponds to measures with no estimated state.

TABLE I LIST

 I OF MAREL-CARNOT SIGNALS: ACRONYM, NAME AND MEASUREMENT FREQUENCY.

	Acronym	Name	Frequency	RD NC
	E TA	Air temperature	20 minutes
	C O21	Corrected dissolved oxygen	20 minutes
	CSAL1	Salinity	20 minutes
	CSAT1	Oxygen saturation percentage	20 minutes
	E CO1	Conductivity	20 minutes
	E LU1	P.A.R. Photosynthetically	20 minutes
		Available Radiation
	E O21	Non-corrected dissolved oxygen 20 minutes
	E PH1	pH	20 minutes
	E TU1	Turbidity	20 minutes
	E VDM	Direction wind	20 minutes
	E VVR	Gust wind speed	20 minutes
	E VVM	Average wind speed	20 minutes
	ECHL1	Fluorescence	20 minutes
	EMAHH	Sea-level (measured)	20 minutes
	ETCO1	Water temperature	20 minutes
	XMAHH Sea-level (calculated)	20 minutes
	C PO1	Phosphate concentration	12 hours
	C NI1	Nitrate concentration	12 hours
	C SI1	Silicate concentration	12 hours

TABLE II 1

 II -NN RR AND OVERLAP SCORES (MEAN AND STANDARD DEVIATION IN PERCENT) FOR M REPRESENTATIVES PER STATE ACCORDING EU-WFD

		LABELING : 2005-2008 DATABASE
	M values	Random selection			K-means
		RR	Overlap		RR	Overlap
	1	68.1 (8.9) 18.4 (6.8)	82.7 (0.1)	11.3 (0.1)
	10	79.4 (4.0) 18.5 (3.7)	89.8 (0.8)	10.0 (0.7)
	100	87.6 (0.9) 12.2 (0.9)	95.1 (0.2)	4.9 (0.2)
	1,000	94.7 (0.2)	5.2 (0.2)	98.4 (0.1)	1.6 (0.1)
			TABLE III
	STFKM INFLUENCE MEASURED WITH A SVM APPROACH: RR AND OVERLAP SCORES (IN PERCENT)
	Sampling		Training		Test
			RR	Overlap	RR	Overlap
	No sampling-SVM	97.9	2.1		92.9	7.0
	Random 1,000-SVM	93.3	6.7		92.2	6.7
	STFKM-SVM	95.4	4.6		92.6	7.4

TABLE VII CORRELATIONS

 VII BETWEEN PARAMETERS AND STATES (IN BOLD: THE HIGHEST CORRELATION COEFFICIENTS).

	State	s 1	s 2	s 3	s 4	s 5	s 6	s 7
	Color	red	green	blue	cyan	pink	yellow	grey
	C NI1	-0.14	-0.24	0.54	-0.24	0.12	0.02	-0.03
	C O21	0.64	-0.16	-0.03	-0.38	0.05	-0.17	0.01
	C PO1	-0.08	-0.03	-0.04	-0.07	-0.02	-0.03	0.57
	C SI1	-0.11	-0.16	0.20	-0.25	0.15	0.04	0.47
	CSAL1	0.12	0.13	-0.36	0.10	-0.25	0.42	0.00
	E LU1	-0.08	0.73	-0.21	-0.24	-0.06	-0.06	0.00
	E TU1	-0.11	-0.14	-0.01	-0.22	0.76	-0.03	-0.05
	E VVR	-0.25	-0.03	0.16	-0.05	0.31	-0.06	-0.03
	ETCO1	-0.50	0.32	-0.37	0.56	-0.14	0.07	0.05
	XMAHH -0.00	0.03	0.02	-0.03	0.01	0.01	0.00

This work was supported in part by IFREMER, in part by the Artois-Picardie Water Agency with a Ph.D. Grant, and in part by the DYMAPHY Interreg IV A 2 Mers Seas and Zee en program. The IFREMERs Department REM and its unit Research and Technological Development have operational control for servicing the MAREL system. K. Rousseeuw is both with the LISIC Lab (LISIC : dOpale address: ULCO/LISIC, BP 719, FR-62228 Calais cedex ) and with French Research Institute for Exploitation of the Sea (IFREMER, Nord, BP 699, FR-62321 Boulogne-sur-Mer, France


In this paper a new model is proposed for Gate Oxide Short defects based on a non-split MOS transistor. Because the MOS is not split, this model allows to simulate minimum transistors in realistic digital circuits. The construction of the model is presented in details using a comprehensive and didactic approach. It is demonstrated that the electrical behavior of the proposed model matches in a satisfactory way the defective transistor behavior.

Introduction

The advent of integrated circuit technology has introduced electronics in many aspect of present-day life. As the use of electronic components increases, the expectation of lower cost, better accuracy, and higher reliability increases. Lower cost and better accuracy is achieved by putting more transistors per unit of silicon, using design automation, increasing device operation speed, and reducing its power consumption. However, these design steps cannot guarantee reliability. In fact, as the circuit density increases, the probability of a manufacturing defect increases. The higher expectation of reliability can only be met by more thorough and comprehensive testing.

Due to the complexity of IC technological process, many physical defects occur during the manufacturing of any system. The typical defects encountered in today technologies and modeled in yield simulators are the socalled spot defects that may cause shorts and/or breaks at one or more of the different conductive levels of the devices. Test generation for any type of defect is obviously not feasible due to the huge amount of CPU time and memory size required. Instead, test generation relies on fault models that are supposed to both represent the defect behavior and allow easy generation of test vectors.

Classical fault models (stuck-at, stuck-open, stuckon, …) have been proved to be efficient for the analysis of many of these faults. However, it is well-known that these fault models cover only partially the spectrum of real failures in today's integrated circuits. The increasing demand of low ppm defect rates requires the derivation of ever more accurate fault models. In particular, a special attention must be paid to defects that exhibit complex behavior not accurately represented by classical fault models and defects with a high probability of occurrence. Gate Oxide Short (GOS) defects belong to both categories since (i) they change some of the electrical features of the transistor, and (ii) they are predominant defects in today technologies in which devices are scaled down and oxide thickness reduced. A number of research works have been conducted in the past years dealing with the electrical characterization [START_REF] Hawkins | Electrical Characteristics and Testing Considerations for Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Soden | Test considerations for Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Hawkins | Reliability and Electrical Properties of Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Syed | Defect Analysis, Test Generation and Fault Simulation for Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Rodriguez-Montanes | Current vs. Logic Testing of Gate Oxide Short, Floating Gate and Bridging Failures[END_REF][START_REF] Segura | Approach to the Analysis of Gate Oxide Shorts in CMOS Digital Circuits[END_REF][START_REF] Segura | A Detailed Analysis of GOS Defects in MOS Transistors: Testing Implications at Circuit Level[END_REF][START_REF] Renovell | A Complete Analysis of the Voltage Behaviour of MOS Transistor with Gate Oxide Short[END_REF][START_REF] Renovell | Analysing the Characteristics of MOS Transistors in the Presence of Gate Oxide Short[END_REF][START_REF] Renovell | Boolean and Current Detection of MOS Transistor with Gate Oxide Short[END_REF] and modeling [START_REF] Syrzycki | Modeling of Spot Defects in MOS Transistors[END_REF][START_REF] Syrzycki | Modeling of Gate Oxide Shorts in MOS Transistors[END_REF][START_REF] Segura | Analysis and Modeling of MOS Devices with Gate Oxide Short Failures[END_REF][START_REF] Champac | Fault Modeling of Gate Oxide Short, Floating Gate & Bridging Failures in CMOS Circuits[END_REF][START_REF] Hao | On the Modeling and Testing of Gate Oxide Shorts in CMOS Logic Gates[END_REF][START_REF] Hao | Analysis of Gate Oxide Shorts in CMOS Circuits[END_REF][START_REF] Segura | A Detailed Analysis and Electrical Modeling of Gate Oxide Shorts in MOS Transistors[END_REF] of this kind of failure.

Both for the characterization and modeling process, one should distinguish between two different types of gate oxide short, i.e. GOS connecting directly the gate to the drain or source or GOS connecting the gate to the channel. For the first type of defect, a realistic model based on the addition of a short resistance between the gate the drain (or source) has been proposed in [START_REF] Syrzycki | Modeling of Gate Oxide Shorts in MOS Transistors[END_REF]. For the second type of defect, several models have been proposed based on the split of the faulty device in several smaller devices [START_REF] Syrzycki | Modeling of Spot Defects in MOS Transistors[END_REF][START_REF] Syrzycki | Modeling of Gate Oxide Shorts in MOS Transistors[END_REF][START_REF] Segura | Analysis and Modeling of MOS Devices with Gate Oxide Short Failures[END_REF][START_REF] Champac | Fault Modeling of Gate Oxide Short, Floating Gate & Bridging Failures in CMOS Circuits[END_REF]. Although these models perfectly match the electrical characteristics of the defect behavior, they suffer from a strong limitation: they cannot be used to study GOS defects in realistic digital CMOS circuits. Indeed, these models are not able to cope with minimum transistors affected by a GOS failure. However, most of the transistors involved in a digital cell library are today designed at minimum-size. It is therefore the objective of this paper to present a new model permitting to handle minimum transistors affected by gate-to-channel shorts.

The paper is organized as follows. GOS failures and their associated behavior are described in section 2. Defect modeling is addressed in section 3 and limitations of the previously proposed models are discussed. In section 4, we introduce the new non-linear non-split model and show how this model permits to accurately represent the defect behavior while preserving the length of the original transistor. Finally, concluding remarks are given in section 5.

Modeling Gate Oxide Short Defects in CMOS Minimum Transistors

Defect behavior

One of the main specificity of a GOS defect is that the pinhole that shorts the gate to one point of the underlying silicon creates a new device in which an important gate current can flow. It is therefore possible to study the IG vs VG characteristics of a defective transistor. As commented by the authors of many previous papers [START_REF] Hawkins | Electrical Characteristics and Testing Considerations for Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Soden | Test considerations for Gate Oxide Shorts in CMOS ICs[END_REF][START_REF] Hawkins | Reliability and Electrical Properties of Gate Oxide Shorts in CMOS ICs[END_REF], two distinct electrical behaviors are observed depending on the defect type. The first type of defects connecting the gate to drain or source diffusion zone is referred as an ohmic defect and exhibits a linear behavior in the IG vs VGS characteristics. In contrast, the second type of defects connecting the gate to transistor channel exhibits a non-linear behavior in the IG vs VGS characteristics. Figure 2 shows typical IG vs VGS characteristics for both the ohmic and non-ohmic types of defects. Regarding the ID vs VDS characteristics, all experiments have demonstrated that the drain current of a defective transistor only slightly resemble the typical MOS transistor drain current. Figure 3.a and 3.b show the typical ID vs VDS characteristics of non-defective and defective transistors. Basically, the defect manifests itself through two main phenomena:

• reduction of the maximum drain current at high VDS,

• appearance of a negative drain current when VDS is small in comparison with VGS. For a given VGS, we will denote IDmax the maximum drain current and IDmin the minimum drain current at VDS=0.

f VDS f ID f (b) 0,8 1,2 1,6 2 2,4 IDmax f IDmin
NMOS with a gate-to-channel GOS Fig. 3: ID vs VDS characteristics

Electrical MOS models of GOS

Linear Non-Split MOS Model

A simple linear model has been proposed in [START_REF] Syrzycki | Modeling of Gate Oxide Shorts in MOS Transistors[END_REF] based on the addition of a short resistance between the gate and drain or source electrodes as illustrated in figure 4. This model is perfectly adequate to represent GOS defects connecting the gate to drain or source diffusion zone in a n-channel transistor. This model has been extended in [START_REF] Hao | On the Modeling and Testing of Gate Oxide Shorts in CMOS Logic Gates[END_REF][START_REF] Hao | Analysis of Gate Oxide Shorts in CMOS Circuits[END_REF][START_REF] Segura | A Detailed Analysis and Electrical Modeling of Gate Oxide Shorts in MOS Transistors[END_REF] to cope with both n-and pchannel transistors. One main advantage of these models is that the length of the original transistor T o is preserved, implying that the model can handle minimum transistors. However, these models cannot apply for gate oxide shorts between the gate and the channel of a transistor as they present a different set of electrical properties. In particular, they are not able to represent the non-linear behavior depicted in the IG vs. VGS characteristics of figure 2.

Non-Linear Split MOS Model

In order to study GOS defects connecting the gate to the channel, an electrical model of the fault-free transistor based on a lumped-element model has been proposed in [START_REF] Syrzycki | Modeling of Spot Defects in MOS Transistors[END_REF]. In this model represented in figure 5, the nondefective channel is split and become a two-dimensional array of MOS transistors. Based on this lumped-element model, a GOS defect can be introduced in the transistor by connecting a short resistance R GOS between the common gate G and one of the internal nodes of the network, denoted (i f ,j f ). Figure 7 gives an example of a central defect introduced in the transistor. The defect resistance and location can then be varied through the value of R GOS and position (i f ,j f ) in the network. This model has been demonstrated to perfectly describe the electrical behavior of gate-to-channel shorts [START_REF] Renovell | A Complete Analysis of the Voltage Behaviour of MOS Transistor with Gate Oxide Short[END_REF][START_REF] Renovell | Analysing the Characteristics of MOS Transistors in the Presence of Gate Oxide Short[END_REF]. For illustration, figure 8.a shows the simulated ID vs VDS characteristics of a defective transistor with a 1Ω central defect. For a given VGS, for instance VGS=2V, we observe a reduction of the maximum drain current IDmax f compared to the non-defective one IDmax o together with the appearance of a negative drain current for small VDS. This model also permits to simulate the IG vs VGS characteristics, as depicted in figure 8.b. Note that these characteristics are in agreement with the experimental behavior observed for gate-to-channel shorts (see figure 2).

line 1 (i,j) S D G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G line i line m col. n col. j col. 1 (1,1)

Limitations

As described in the previous section, the split MOS model for GOS perfectly describes the behavior of gateto-channel defects and for this reason has been widely studied by different authors [START_REF] Renovell | A Complete Analysis of the Voltage Behaviour of MOS Transistor with Gate Oxide Short[END_REF][START_REF] Renovell | Analysing the Characteristics of MOS Transistors in the Presence of Gate Oxide Short[END_REF][START_REF] Renovell | Boolean and Current Detection of MOS Transistor with Gate Oxide Short[END_REF][START_REF] Syrzycki | Modeling of Spot Defects in MOS Transistors[END_REF][START_REF] Syrzycki | Modeling of Gate Oxide Shorts in MOS Transistors[END_REF][START_REF] Segura | Analysis and Modeling of MOS Devices with Gate Oxide Short Failures[END_REF][START_REF] Champac | Fault Modeling of Gate Oxide Short, Floating Gate & Bridging Failures in CMOS Circuits[END_REF]. However this model suffers from severe limitations as discussed below.

As everybody knows, SPICE simulations are performed using card models describing the transistor parameters in the targeted technological process. Using these card models, transistors of any length can be simulated as far as they are larger than the minimum (L min ) allowed by the technological process. Any attempt to simulate transistors smaller than the minimum would make the simulation results (if any) unreliable.

Consequently, when using the split MOS model for simulating GOS defects, the elementary MOS transistor obtained after the splitting of the original transistor cannot have a length smaller than L min . This obviously implies that the original transistor cannot be at minimum length. Even the smallest network composed of two transistors in series as suggested in [START_REF] Segura | Analysis and Modeling of MOS Devices with Gate Oxide Short Failures[END_REF][START_REF] Champac | Fault Modeling of Gate Oxide Short, Floating Gate & Bridging Failures in CMOS Circuits[END_REF] corresponds to an original transistor with length equal or higher that 2L min . It is therefore impossible to study minimum transistors affected by GOS failures using a split MOS model. This is an extremely strong limitation taking into account that in a classical digital cell library, all transistors are designed at minimum length. In this context, although the split MOS transistor permits to accurately model the defect behavior, it clearly appears inadequate to simulate GOS defects in realistic digital circuits.

Non-Linear Non-Split MOS Model

This section describes the new model we developed for GOS defects presenting a non-linear behavior, i.e. gate-to-channel shorts. Our first objective is not to split the original transistor as proposed in [START_REF] Hao | On the Modeling and Testing of Gate Oxide Shorts in CMOS Logic Gates[END_REF] for the linear model in order to preserve the length of the original transistor. Our second objective is to introduce non-linearity in the model in order to properly represent the behavior of this type of defect. Consequently, our nonlinear model is constructed by adding non-linear components to the original transistor.

For the demonstration, we reuse the original transistor of figure 5 with L o =2.1µm and W o =3.5µm in a 0.25µm technology. Because the split model has been proved very accurate [START_REF] Renovell | A Complete Analysis of the Voltage Behaviour of MOS Transistor with Gate Oxide Short[END_REF][START_REF] Renovell | Analysing the Characteristics of MOS Transistors in the Presence of Gate Oxide Short[END_REF], we consider that the characteristics of the faulty transistors are given by figure 8. Note that these characteristics obtained with the split model are just used here as a reference.

Our non-linear non-split MOS model is now constructed in three steps with the aim to mimic as precisely as possible the characteristics of figure 8.

Step 1

As commented in section 2.2, a first specificity of a transistor affected by a GOS failure resides in the reduction of the drain current. For a given VGS, we observe that the maximum drain current in the defective transistor is smaller than the maximum drain current in the non-defective one: IDmax f < IDmax o (see figure 6 and8). Consequently, we start our model construction with a transistor with the same length L o than the original transistor but a smaller width W m (W modified) to reduce the drain current. Note that the length of L o of the original transistor is preserved as required by our objectives; only the width W m is modified and calculated to match the faulty current IDmax f . Figure 9 illustrates this first step. Step 2

The second specificity of a transistor affected by a GOS failure is the appearance of a negative drain current for low VDS. Consequently, the second step of our model construction consists in connecting an additional component allowing a current to flow from the gate to the drain. This additional component should act only at low VDS but not modify the I-V characteristics at high VDS. We propose to insert an additional transistor T a connected as illustrated in figure 10. For the sake of clarity, the behavior of the above circuit is analyzed for VGS=2V. The results will be then extended to any VGS. When VDS is higher than VGS=2V, the additional transistor T a is OFF and has no influence on the I-V characteristics. When VD is smaller than VGS=2V, the additional transistor T a has its drain in node G and its source in node D as illustrated in figure 10. Under these conditions, the transistor is ON and operates as a rectifier with Vgd a =0. The corresponding Ids a vs Vds a characteristics are given in figure 11 Step 3

Finally, the last specificity of a transistor affected by a gate-to-channel short is the presence of a non-linear gate current. Consequently, the last step in our model construction consists in connecting another additional transistor T b allowing a current to flow from the gate to the source, as illustrated in figure 14.

Conclusions

In this paper, we proposed a new non-linear non-split MOS model for gate-to-channel GOS defects. The fundamental idea is not to split the original transistor in order to handle minimal-length transistors and therefore enable the simulation of GOS defects in realistic digital circuits. Based on the analysis of the electrical defect behavior, we developed a comprehensive method for the construction of the model. The model is constructed by first modifying the width of the original transistor, and then inserting two additional transistors between the gate and drain/source of this transistor:

• the width modification (W m ) of the original transistor allows to adjust the level of the drain current at high VDS (IDmax f ), • the conductance W a /L a of the first additional transistor allows to adjust the level of the negative drain current at low VDS (IDmin f ), • the conductance W b /L b of the second additional transistor allows to adjust the level of the gate current (IG f ). It has been shown that the behavior of the proposed model matches in a satisfactory way the behavior of a defective transistor model.

Finally, future work will concentrate on the tuning of the model parameters (W m , W a /L a , W b /L b ) to reproduce GOS defects of various resistance, location and size.

Fig. 2 :

 2 Fig.2: IG vs VGS characteristics of a MOS with a GOS

Fig. 4 :

 4 Fig.4: Linear non-split MOS model

Fig. 6 :

 6 Fig.6: Comparison of the ID vs VDS characteristics of the original transistor and the split MOS model

Fig. 7 :

 7 Fig.7: Split model of a NMOS with a GOS

8 :

 8 Simulated I-V charateristics of a defective NMOS using the split model

  original transistor (b) modified transistor Fig.9: Non-linear non-split model -Step 1

Fig. 10 :

 10 Fig.10: Non-linear non-split model -Step 2

  .a. It is to note that Vds a = VGS -Vds m so in figure 11.b and 11.c the characteristics are flipped and shifted to have compatible x-axis with figure 8.a. Ids a vs Vds m =(VGS-Vds a ) Fig.11: I-V characteristics of the additional transistor T a In node D, we can write the current equation as follows: ID f + Ids a = Ids m ⇒ ID f = Ids m + (-Ids a ) That means that the current ID f flowing in the complete device can be obtained by adding the two characteristics of figure 12.a and 12.b. The resulting ID f current is represented in figure 12.c.It is now just necessary to adjust the value of the IDmin f current of figure12.c to exactly fit the current IDmin f of figure8.a. The value of this current actually just depends on the additional transistor size T a since the current Ids m is equal to 0 at point VDS=0. Consequently, the current IDmin f in our model is simply set by tuning the W a /L a parameter of the additional transistor T a . Note that in figure12.c, the W a /L a parameter has been tuned to match the IDmin f value of figure8.a with W a =3.6µm and L a =2.1µm.

Fig. 13 :

 13 Fig.13: ID f vs VDS characteristics of the non-linear non-split model

Fig. 14 :

 14 Fig.14: Non-linear non-split model -Step 3 This second additional transistor actually acts symmetrically to the first additional transistor, creating an IG current in the model. As performed in the previous step with the IDmin f current, we adjust the W b /L b conductance of the second additional transistor T b to mimic the IG f vs VGS f characteristics of figure 8.b. For example, when VDS=0, the current IG is maximum and figure 15 gives the equivalent circuit. The current IGmax f is equal to Ids a +Ids b . The current Ids a has already been fixed and so IGmax f is adjusted through W b /L b . To match the IGmax f value of figure 8.b, we set the additional transistor T b to with W b =3.6µm and L b =2.1µm.

Fig. 15 :Fig. 16 :

 1516 Fig.15: Equivalent circuit for the model at VDS=0Finally, figure16gives the complete set of the IG vs VGS characteristics simulated with our model. These characteristics correctly reproduce the IG-VGS behavior of a gate-to-channel GOS.

. GOS Failures 2.1. Defect description

	M. Renovell, J.M. Gallière, F. Azaïs and Y. Bertrand 2A GOS is a transistor defect that cause a relatively
	low impedance path between the CMOS gate and the Laboratoire d'Informatique Robotique Microélectronique de Montpellier LIRMM-UMII underlying silicon. Depending on the defect location, the Université de Montpellier II: Sciences et Techniques du Languedoc GOS can be seen as a short circuit between the gate and UMR C5506 CNRS -161, rue Ada 34392 Montpellier Cedex 5 France transistor channel or as a short circuit between the gate Tel: (33)467418523 -Fax: (33)467418500 -Email: renovell@lirmm.fr and drain or source diffusion zones. Figure 1 illustrates
	these two types of defect. For both types, an undesired
	path of current through the oxide of the gate appears thus
	creating a violation of the gate isolation principle. It is
	generally admitted that GOS may have different origins
	such as lithographic defect on mask, field failure due to
	ESD…
		Gate			Gate
	Source	poly	Drain	Source	poly	Drain
	N+		N+	N+		N+
	P			P
	(a) gate-to-channel	(b) gate-to-source (drain)
		Fig.1: Pinhole in a MOS gate