Pierre Ailliot

Valrie Monbet

Markov-switching autoregressive models for wind time series

Keywords: Stochastic weather generators, Wind time series, Markov-switching autoregressive model, Multiscale model, Overdispersion

published or not. The documents may come

Introduction

This paper develops stochastic models for wind time series over different time scales. A particular impetus for this study was the need to generate realistic wind time series at different meteorological stations located nearby potential wind farms in France, with the aim of assessing various quantities related to the wind power production (see e.g. [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF], [START_REF] Castino | Stochastic modelling of wind velocities time series[END_REF]). More generally, the proposed model have many other risk forecasting applications since they can be used to generate synthetic time series of unlimited length for a location consistent with historical data at the same location. The synthetic time series can then be used as inputs into environmental and ecosystem models (see e.g. [START_REF] Ailliot | Long term object drift forecast in the ocean with tide and wind[END_REF]), among many other applications (see also [START_REF] Monbet | Survey of stochastic models for wind and sea-state time series[END_REF]).

[20] provides a recent review of stochastic models which have proposed in the literature for wind time series. The most classical approach consists in applying the general Box-Jenkins methodology ( [START_REF] Box | Time series analysis, forecasting and control[END_REF]). Wind times series are generally not stationary, with typically important seasonal and daily components but also interannual variability. The first step of the Box-Jenkins methodology consists in achieving stationarity. Interannual components are generally neglected and various methods are available in the literature for modeling daily and seasonal components, for instance differentiation or scaling with mean and variance functions which evolve periodically (see e.g. [START_REF] Brockwell | Introduction to Time Series and Forecasting[END_REF] and [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF]). Another common approach for treating seasonality in meteorological applications consists in blocking the data by short time periods, typically one month or one season depending on the amount of data available, and then assuming that the changes due to the season can be neglected on this time period. Then, after achieving stationarity, a Box-Cox transformation is generally applied to get a time series with marginal distribution close to a Gaussian distribution and an ARMA model is fitted to this transformed time series (see e.g. [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF]).

Box-Jenkins methodology generally leads to a good description of the marginal distribution and second order structure of the original wind time series but fails to reproduce non-linearities which may exist in the dynamics. One well known source of non-linearity in many meteorological time series is induced by the existence of "weather types". They correspond to typical pressure and frontal patterns and induce regime shifts in local weather conditions. For example, for the specific time series considered in this paper, at least two regimes can be easily identified when looking at the data. In the first one, anticyclonic conditions are prevailing leading to steady and low wind speed, whereas in the second one moving low pressure systems are dominating leading to more important time variability in the wind conditions. This induces heteroscedasticity in wind time series, and GARCH models have been proposed in this context (see [START_REF] Toll | Autoregressive conditional heteroscedasticity in daily wind speed measurements[END_REF]). In this work, in order to get a more physically-based model, we introduce explicitly the regime shifts through a hidden variable. As concerns meteorological applications, the idea of introducing a latent variable which represents the weather type goes back to [START_REF] Zucchini | A hidden Markov model for space-time precipitation[END_REF] where Hidden Markov Models (HMM) were proposed for modelling the space-time evolution of daily rainfall. HMM have then been extensively used for modelling rainfall (see [START_REF] Ailliot | Space time modeling of precipitation using a hidden markov model and censored gaussian distributions[END_REF] and references therein). HMM are characterized by various conditional independence assumptions on the joint dynamics of the hidden variable and the observations. These assumptions imply in particular that successive observations are conditionally independent given the hidden variable. In practice, it means that all the dynamics should be explained by the weather type, and such assumptions seems rather unrealistic for wind time series since the correlation between successive observations is generally high.

In this paper, we propose to use Markov-Switching AutoRegressive (MS-AR) models. This family of model, which was initially proposed in [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF] to describe econometric time series, is a generalization of both HMM and autoregressive models. Indeed, they combine different autoregressive models to describe the evolution of the process at different periods of time, the transition between these different autoregressive models being controlled by a hidden Markov chain like in HMM. Other forms of regime dependent autoregressive models have been proposed in the literature for modelling meteorological time series (see e.g. [START_REF] Zwiers | Regime-dependent autoregressive time series modeling of the southern oscillation[END_REF], [START_REF] Parlange | An extended version of the richardson model for simulating daily weather variables[END_REF]).

Hereafter, we focus on a particular time series of wind speed measured on the Island of Ouessant (48 0 27 ′ 36 ′′ N, 5 0 6 ′ 0 ′′ W). This time series, which is shown on Figure 1, consists of 51 years of data, from 1948 to 1998, with a data every 6 hours which corresponds to the mean wind speed on a 20 minutes time period. We have chosen to consider a long time series in this study in order to discuss the modelling of interannual variability, but good results were also obtained with the methodology described in this paper on shorter wind time series at various locations in France with different climatologies. The results are satisfactory except that the model tends to underestimate the observed interannual variability in the wind condition. This is the well known "overdispersion" phenomenon which is a weakness of many stochastic weather generators (see e.g. [START_REF] Katz | Overdispersion phenomenon in stochastic modeling of precipitation[END_REF]). One possible explanation may be the existence of interannual variations which are not taken into account in the model. In Section 5, we thus explore interannual components and show that there is a clear trend in the probability of occurrence of the different weather type. This lead us to propose another MS-AR model which includes a trend in the parameter values. It is shown that this new MS-AR permits a better description of the observed interannual variability although it is still underestimated. Finally we conclude in Section 6.

MS-AR models 2.1 Model description

A MS-AR process is a discrete-time process with two components {S t , Y t } where, for our particular application, {Y t } denotes the wind speed process with values in (0, +∞) and S t ∈ {1, ..., M} represents the latent weather type at time t. A MS-AR process is then characterized by the two conditional independence assumptions below:

• the conditional distribution of S t given the values of {S t ′ } t ′ <t and {Y t ′ } t ′ <t only depends on the value of S t-1 . In other terms, we assume that the weather type {S t } is a first order Markov chain which evolution is independent of the past wind conditions.

• the conditional distribution of Y t given the values of {Y t ′ } t ′ <t and {S t ′ } t ′ ≤t only depends on the values of S t and Y t-1 , . . . , Y t-p . For our particular application, it means that the wind speed process {Y t } is an autoregressive process of order p ≥ 0 which coefficients evolve in time with the weather type sequence.

When p = 0, we retrieve the usual HMM, and the various conditional independence assumptions are summarized by the directed graph below for p = 1:

• • • → S t-1 → S t → S t+1 → • • • ↓ ↓ ↓ • • • → Y t-1 → Y t → Y t+1 → • • •
MS-AR models were initially introduced in [START_REF] Hamilton | A new approach to the economic analysis of nonstationary time series and the business cycle[END_REF] to describe econometric time series, the regimes corresponding to the different states of the economy, and then used for other applications (see e.g. [START_REF] Ephraim | Hidden markov processes[END_REF] and references therein). MS-AR models for wind time series were initially proposed in [START_REF] Ailliot | Modèles autorégréssifs à changements de régimes Makoviens[END_REF] and [START_REF] Monbet | Survey of stochastic models for wind and sea-state time series[END_REF] and then, in different contexts in [3] [22].

In the usual applications of HMM and MS-AR models, the hidden Markov chain {S t } is supposed to be homogeneous, in which case the transition probabilities P (S t = s ′ |S t-1 = s) are constant in time and the evolution of {S t } is parametrized by the transition matrix Q = (q s,s ′ ) s,s ′ ∈{1,...,M } with q s,s ′ = P (S t = s ′ |S t-1 = s). HMM with non-homogeneous hidden Markov chain were also proposed for meteorological applications, for example in [START_REF] Zucchini | Hidden Markov models for time series : an introduction using R[END_REF] to describe non-stationary components in time series of wind direction and in [START_REF] Hughes | A class of stochastic models for relating synoptic atmospheric patterns to local hydrologic phenomenon[END_REF] to relate the large circulation to local rainfall conditions. In the next sections, we also propose various MS-AR models with non-homogeneous hidden Markov chain in order to describe seasonal and interannual variations.

As concerns the autoregressive models, the most standard MS-AR model is obtained using standard AR(p) models with Gaussian innovations. If S t = s t , it is assumed that

Y t = a (st) 0 + a (st) 1 Y t-1 + . . . + a (st) p Y t-p + σ (st) ǫ t (1)
where a

(s) 0 , a (s)
1 , ..., a

p , σ (s) ∈ R p+1 × (0, +∞) denotes the unknown parameters of the AR(p) model which describes the evolution of the observed process in the regime s ∈ {1, ..., M} and {ǫ t } is a sequence of independent and identically distributed Gaussian variable with zero mean and unit variance independent of the Markov chain {S t }. In other words, it is assumed that the conditional distribution P (Y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p , S t = s t ) is a Gaussian distribution with conditional mean and variance given respectively by

E(Y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p , S t = s t ) = a (st) 0 + a (st) 1 y t-1 + . . . + a (st) p y t-p (2)
var(Y t |Y t-1 = y t-1 , . . . , Y t-p = y t-p , S t = s t ) = σ (st) 2 (3)
For the application considered in this paper, {Y t } is a process with positive values, and the model with conditional Gaussian distribution may not be appropriate in such situation since it does not permit to restore the constraint Y t ≥ 0. In [START_REF] Ailliot | Modèles autorégréssifs à changements de régimes Makoviens[END_REF], it was proposed to replace the Gaussian distribution by a Gamma distribution and keep (2) and (3) for the conditional moments (the Gamma distribution is also characterized by its two first moments). This model was fitted to various wind time series, and generally good results

were obtained (some of these results are reported in [START_REF] Monbet | Survey of stochastic models for wind and sea-state time series[END_REF]). However, one drawback of the parametrization based on the Gamma distribution is that the additional constraints

a (s)
k > 0 are needed, for k ∈ {0, . . . , p} and s ∈ {1, . . . , M}, in order to ensure that the conditional mean of the Gamma distribution is positive. The tests that we have done on various time series indicates that fitting the MS-AR model with Gaussian innovations (1) generally leads to selecting AR models of order p = 2 with autoregressive coefficients a (s) 2 < 0 (see for example the numerical results given in the next sections) and in such situation using the model with conditional Gamma distribution may not be appropriate.

Hereafter we consider only models with conditional Gaussian distributions. It permits to save computational time and also avoid numerical problems which may exist when using the Gamma distribution. This and other modelling issues related to the parametrization of the autoregressive models will be further discussed in the next sections.

Statistical inference

The most classical method for fitting a MS-AR model, for given values of M and p, consists of using the Expectation-Maximization (EM) algorithm. It was proven that this algorithm, which was first introduced in [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains[END_REF] for HMM and then generalized to other models with hidden variables in [START_REF] Dempster | Maximum likelihood from incomplete data[END_REF], converges to a maximum of the likelihood function under general conditions (see [START_REF] Wu | On the convergence properties of the em algorithm[END_REF]). The description of the particular form of this algorithm for MS-AR models with Gaussian innovations and homogeneous Markov chain can be found in [START_REF] Krolzig | Markov-switching vector Autoregressions. Modelling, statistical inference and applications to business cycle analysis[END_REF]. This is an iterative procedure, starting from an initial value θ(0) for the parameters. Each iteration consists of 2 steps:

• E step: computation of an auxiliary function R(θ, θ(n)) which is defined as the conditional expectation of the complete likelihood given the observations and the current value of the parameters. For all the models considered in this paper, this step can easily be performed using the classical forward-backward recursions (see e.g. [START_REF] Cappé | Inference in hidden Markov models[END_REF] and [START_REF] Zucchini | A hidden Markov model for space-time precipitation[END_REF]).

• M step: computation of θ(n + 1) = argmax θ R(θ, θ(n)). Depending on the MS-AR model under consideration, there are not always analytical expression for θ(n + 1), in which case a numerical optimization procedure is required. In order to get an efficient EM algorithm, it is important to implement carefully the optimization problem. In particular, it is often possible to break the optimization problem into several lower dimensional optimization problems which are much quicker to solve.

More precisely, for all the models considered in this paper, it is possible to separate the parameters related to the evolution of the hidden Markov chain, θ S , and the parameters related to the evolution of the observed process in each regime s ∈ {1, . . . , M}, denoted θ

(s) Y , such that θ = θ S , θ (1)
Y , . . . , θ (M ) Y

. For example, for the MS-AR model with homogeneous Markov chain and AR(p) models with Gaussian distribution (1), we have θ S = (q s,s ′ ) s,s ′ ∈{1,...,M } , with the usual constraints to get a well defined transition matrix, and θ

(s) Y = a (s) 0 , a (s) 1 , . . . , a (s)
p , σ (s) . Then we have a decomposition of the form

R(θ, θ(n)) = R S (θ S , θ(n)) + R (1) Y (θ (1) Y , θ(n)) + . . . + R (M ) Y (θ (M ) Y , θ(n))
which leads to M + 1 separate optimization problems on reduced dimension spaces.

There may exist analytic expression for some of them, e.g. when the hidden Markov chain is homogeneous or when the autoregressive models are parametrized using [START_REF] Ailliot | Modèles autorégréssifs à changements de régimes Makoviens[END_REF].

Otherwise, a standard quasi-Newton algorithm has been used in this work, with an appropriate treatment of the various constraints on the coefficients.

The EM algorithm has several well-known limitations. First, it may converge to a noninteresting local maximum of the likelihood function. In practice, it means that a careful choice of the starting value has to be made; this is further discussed in the next sections.

Another drawback is its slow rate of convergence near the maxima, where it is known that a usual quasi-Newton algorithm is more efficient. An additional advantage of using a quasi-Newton algorithm is that it provides directly an approximation of the Hessian of the log-likelihood function, and thus gives useful information on the variance of the estimates. On the other hand, quasi-Newton algorithms are generally more sensitive to the choice of the starting value and require some programming efforts since they need the gradient of the function to optimize as input to be efficient. Such algorithms have not been implemented in this work.

The stability of MS-AR models and the asymptotic properties of the Maximum Likelihood Estimates (MLE) in HMM and MS-AR models have been studied extensively in the recent years (see e.g. [START_REF] Cappé | Inference in hidden Markov models[END_REF] and references therein). In particular, general conditions which ensure consistency and asymptotic normality of the MLE for MS-AR model with homogeneous hidden Markov chain and the autoregressive models described above with Gaussian or Gamma conditional distributions can be found in [START_REF] Krishnamurthy | Consistent estimation of linear and non-linear autoregressive models with markov regime[END_REF] and [START_REF] Ailliot | An autoregressive model with time-varying coefficients for wind fields[END_REF], but the existing results do not apply to the non-stationary MS-AR models considered in this paper. We have performed various simulation experiments which indicate however that the MLE have a reasonable behaviour for the non-stationary models considered in this paper.

Another important problem in practice, which has received lots of attention in the last few years, is the problem of model selection which aims at finding the "optimal" value of p and M (see e.g. [START_REF] Cappé | Inference in hidden Markov models[END_REF] for a recent review). Hereafter, we have chosen to use the Bayes Information Criterion (BIC) as a first guide. Although its use is not justified for MS-AR models from a theoretical point of view, we found that it generally permits to select parsimonious models which fit the data well. It is defined as

BIC = -2 log L + k log N
where L is the likelihood of the data, k is the number of parameters and N is the number of observations. It can be easily computed from the likelihood which is a natural output of the forward recursions performed in the E-step of the EM algorithm.

MS-AR model for monthly data

A classical approach for treating seasonality for meteorological time series consists in blocking the data by month and fit a separate model each month, assuming that the different realisations of the same month over different years are independent realizations of a common stochastic process. This approach is used in this section and we discuss the results obtained on the wind time series introduced in Section 1.

Model description

Even when focusing on a monthly time period, diurnal fluctuations generally imply that wind time series are not stationary (see the periodogram on Figure [START_REF] Bouette | Wind in ireland: long memory or seasonal effect?[END_REF]). According to Figure 2, for the time series considered in this work, the wind speed is generally higher during the day than during the night, with maximum mean value at noon and the diurnal variations are more important in summer than in winter due to the higher diurnal variations in the temperature. A classical approach (see e.g. [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF], [START_REF] Parlange | An extended version of the richardson model for simulating daily weather variables[END_REF]) for modelling wind time series with diurnal components consists in scaling the data by subtracting the periodic mean function and then eventually dividing by the periodic standard variation function and then assume that the residual time series is an AR(p) process. This is similar to assume that the wind time series is a non-homogeneous AR(p) process with periodically evolving coefficients. Here we propose to use such non-homogeneous autoregressive models in each regime and replace

(2) by ( 4)

E (Y t |S t = s t , Y t-p = y t-p , ..., Y t-1 = y t-1 ) = a (st) 0 (t) + a (st) 1 y t-1 + ... + a (st) p y t-p (4)
where, for s ∈ {1, ..., M}, a (s)

0 (t) = α (s) 0 + α (s) 1 cos 2πt T d -α (s) 2
with the unknown parameters α

(s) 0 ∈ R, α (s)
1 ≥ 0 and α (s)

2 ∈ [0, 2π[ and T d represents the number of observations per day in such a way that (4) defines a periodic function with period one day.

In [START_REF] Ailliot | Space time modeling of precipitation using a hidden markov model and censored gaussian distributions[END_REF], only the intercepts of the AR models are assumed to vary with time. Loosely speaking, it permits to model that the mean of the wind speed exhibits diurnal variation but not its variance. According to Figure 2, this seems to be realistic for the particular time series considered in this work. This model could be obviously generalized by assuming that the other coefficients of the autoregressive models are periodic functions, but the various attempts that we have done in this direction for our particular data set did not improve the results obtained with the simplest model ( 4). The coefficient α (s) 1 is related to the amplitude of the diurnal variations in regime s, whereas α (s) 2 is associated to the time in the day when the wind speed is maximum. In the limiting case α (s) 1 = 0 we retrieve the homogeneous model [START_REF] Ailliot | Long term object drift forecast in the ocean with tide and wind[END_REF]. The model ( 4) allows these characteristics to be different in the weather types. For example, for our dataset we expect more important diurnal variations when anticyclonic conditions are prevailing than in cyclonic conditions and such behaviour can not be restored by the more conventional approach discussed in [START_REF] Brown | Time series models to simulate and forecast wind speed and wind power[END_REF].

Hereafter, we will also consider another approach which consists in including the diurnal component in the dynamics of the hidden Markov chain and assume that the transition probabilities are periodic functions. Such approach was initially proposed in [START_REF] Zucchini | A hidden Markov model for space-time precipitation[END_REF] for wind direction and then used in [START_REF] Ailliot | Modèles autorégréssifs à changements de régimes Makoviens[END_REF] for wind speed. In this work, we have considered simple parametric functions derived from the multinomial logit model and assumed that

P (S t = s ′ |S t-1 = s) ∝ q s,s ′ exp κ s ′ cos 2π T t -φ s ′ (5)
where Q = (q s,s ′ ) s,s ′ ∈{1,...,M } is a stochastic matrix and, for s ∈ {1, ..., M}, κ s ≥ 0 and φ s ∈ [0, 2π[ are unknown parameters. Again, T = T d represents the number of observations per day in such a way that (5) defines a periodic function with period one day. The limiting case κ s = 0 for s ∈ {1...M} corresponds to the homogeneous case whereas for high values of κ s the conditional distribution ( 5) is concentrated around φ s .

Parameter estimation

The various models obtained by combining the different parametrizations discussed above for the hidden Markov chain, which can be homogeneous or not, and the autoregressive models which can be homogeneous or not have been fitted to the 12 data sets obtained by blocking the data by month. There are some missing data in the original time series.

When only one data point is missing, a single linear interpolation method has been used to fill in the gap using adjacent values. It leads to a new time series with missing values only in 1986 and 1991, and these two years include a long time period with no data. We have thus decided to remove these two years from the original time series in order to facilitate the statistical inference. Finally, for each month 49 realizations of length N = 4 * Nd remain, where Nd represents the number of day in the month under consideration, in order to fit and validate the models. These realizations are supposed to be independent and the likelihood function which we consider is obtained as the product of the likelihood over the 49 realizations. The likelihood function has been maximized using the EM algorithm.

In practice, we consider models with a number of regimes M = 1, . . . , 5 and autoregressive models of order p ∈ {1, 2}. In order to initialize this algorithm with realistic parameter values, and thus avoid convergence to non-interesting maxima and save computational time, we have used the nested nature of the models. For example, the non-homogeneous models were initialized using the parameter values of the corresponding fitted homogeneous models and the models of order p = 2 using the models of order p = 1. When such initialization was not available (for example for the homogeneous models of order p = 1), the EM has been initialized using several starting values chosen randomly in a set on physically realistic parameter values.

Results

In this section, we focus on the months of January and July since the results obtained for these two months are representative of the ones for the other months.

According According to Table 2, the first regime corresponds to periods with steady wind conditions, with a low standard deviation for the innovation σ (s) and also a slower decrease to zero of the autocorrelation functions than in the other regimes, whereas the third regime corresponds to periods with important temporal variability in the wind conditions. The comparison of the means of the stationary distributions in the different regimes also indicates that higher wind speed are generally observed in periods with high variability than in period with low variability. The transition matrix exhibits high values on the diagonal and thus the different regimes are relatively persistent (the mean duration of sojourns varies between 2.69 days in the regime 2 and 5.86 days in regime 3). There are also some very small transition probabilities : for example most of the time the Markov chain will transit from regime 1 to regime 3 through regime 2 and vice-versa.

The stationary distribution of the hidden Markov chain indicates that the three regimes have almost the same probability of occurrence.

A useful tool to confirm visually the interpretation of the states consists in computing the smoothing probabilities defined as the conditional distribution of the hidden state given all the observations (y 1 , ..., y N ) available in a given month

P [S t = s|Y 1 = y 1 , ..., Y N = y N ]
for s ∈ {1, ..., M}. A complementary information is provided by the most likely hidden sequence computed using the Viterbi algorithm which leads to a segmentation of the observed time series according to the different regimes (see [START_REF] Zucchini | Hidden Markov models for time series : an introduction using R[END_REF]). An example of such Table 2: Estimated parameters for the homogeneous model with M = 3 regimes and autoregressive models of order p = 2 together with the stationary distribution π (s) of the Markov chain and the mean µ (s) of the stationary solution of the AR models. Asymptotic standard errors computed from the observed information matrix are given in brackets.

Results for January.

segmentation is shown on Figure 3 : we retrieve periods when one of the regimes is clearly the more likely (see for example the stormy period in the middle of the month) and transition periods with several regimes having similar probabilities.

Figure 4 shows the distribution of the wind direction in the different regimes identified using the Viterbi algorithm. The third regime is mainly associated with wind from the South-West : it may correspond to cyclonic conditions when quickly evolving low-pressure systems are coming from the Atlantic ocean. The two other regimes have similar distributions for the wind direction and can be associated with all wind directions. Looking at the distribution of other meteorological variables, such as the sea-level pressure, may help refining the meteorological interpretation of the different regimes.

To further validate the model, we have checked its ability to simulate realistic wind time series since this is an important aspect for the applications which have motivated this work.

For that, we have generated artificial time series from the model and we have compared various statistics computed from these artificial sequences with those computed from the data. Typical results are shown in Figure 5. In order to assist visual comparison, 95%

prediction intervals for the fitted model have been superimposed where these quantities have been computed using Monte Carlo methods. The limits of the intervals correspond to the 2.5th and 97.5th percentiles from 1000 independently sequences of 49 months of January simulated using the fitted model. Results for January.

Figure 5 shows the results obtained for the marginal distributions, the autocorrelation function and the distribution of the sojourn durations above and below some selected thresholds. They are among the most common criteria used to validate stochastic weather generators. These results were compared with the ones obtained using the Box-Jenkins methodology, and we could identify several advantages of using MS-AR models. First, the MS-AR model is able to reproduce the marginal distribution of the process without applying an initial transformation, such as the Box-Cox transformation, to achieve normality. This is not surprising given the distributional flexibility inherent in hidden Markov modelling. However, the model tends to overestimate the probability of low wind speed and can even simulate negative wind speed. Nevertheless, the results remains satisfactory, especially if the simulated time series are used as input to simulate the behaviour of a system with a low sensitivity to light wind conditions, such as the power value of a wind turbine. If low wind speed are important for a particular application or if the statistical inference leads to creating artificial state just to account for the skewness of the marginal distribution then using the model with conditional Gamma distributions discussed in Section 2 or applying a preliminary Box-Cox transformation could be more appropriate. Figure 5 shows that the fitted MS-AR model permits also to reproduce the autocorrelation function and the distribution functions of the sojourn durations above and below some selected thresholds since the sample distribution function always lie in the 95% prediction interval computed from the model. We could not get such good results using Box-Jenkings methodology as concerns the sojourn durations and there are good theoretical reasons for that. Indeed, Box-Jenkins methodology is based on ARMA models and thus assume that, after eventual increasing transformation, the process is Gaussian.

It entails some symmetry in the dynamics of the time series and that the behaviour for low wind speed should be similar to the one for high wind speed. In particular, the durations of the sojourns below the quantile of order p should have the same distribution that the ones above the quantile of order 1p. Figure 5 indicates that this is not true for the time series considered in this paper and that the durations of the excursions below the 25% quantile tend to be longer that the ones above the 75% quantile : we find again that the time series exibits more variability at high level than at low level. MS-AR models mix different AR models and thus allow, for example, different dynamics at low and high levels. It leads to a good reproduction of the sojourn durations (see Figure 5). Similar results were obtained for other months. For the months with important diurnal variations, we also checked that the fitted MS-AR model with homogeneous hidden which corresponds to the 25% quantile of the marginal distribution. Bottom right: logsurvivor-functions of the time duration of the sojourns above 12ms -1 which corresponds to the 75% quantile of the marginal distribution. Time is expressed in days. The full lines corresponds to the sample functions and the dashed lines to the fitted model with a 95% prediction intervals (dotted line). The distributions for the fitted model were obtained by simulation. Results for January.

Markov chain but non-homogeneous AR models can reproduce the characteristics of these variations. For example, Figure 6 shows that the fitted model can reproduce both the fluctuations of the mean wind speed and the peak at 1day -1 in the periodogram for the month of July. The model has also been successfully validated at other locations with different climatologies. In particular, a MS-AR model was fitted on wind data in the Aegean Sea and a systematic comparison with the Box-Jenkins methodology and a nonparametric resampling method is performed in [START_REF] Monbet | Survey of stochastic models for wind and sea-state time series[END_REF]. It is shown that the MS-AR model and the non-parametric resampling method improve the results obtained with the Box-Jenkins methodology and both lead to similar results. An advantage of using MS-AR models over non-parametric resampling methods is their interpretability which leads to open structure and allows for more physical models as it is illustrated in the next sections.

MS-AR model with seasonal components

In the previous section, a separate MS-AR model was fitted each month. For many applications, it is necessary to have a model which can simulate the wind speed on a yearly basis. A straightforward combination of the monthly models would lead to a yearly model where the parameters vary as a step function with a break at the beginning of each month. In this section, we propose including seasonality in a more appropriate way into the model.

The results obtained when fitting the MS-AR models introduced in the previous section to each of the 12 months indicate that a MS-AR with M = 3 regimes, homogeneous hidden Markov chain but non-homogeneous autoregressive models is the simplest model which gives satisfactory results for all months. As discussed in Section 3, simplest homogeneous MS-AR models also provide a good description of wind conditions in winter, when daily components can be neglected, but for simplicity reasons we have decided to keep the same model for the different months. Then, in order to be able to follow the seasonal evolution of the parameters, the regimes have been numbered increasingly according to their conditional standard deviations σ (s) , the first regime corresponding to wind conditions with low variability whereas the third one to higher variability. In May the first two states were inverted in order to make the time evolution of the coefficients more consistent.

Figure 7 provides a synthetic view of the seasonal evolution of some of the coefficients of the fitted models and summarizes important features of the climatology. First, the time evolution of α (s) 0 and σ (s) indicates respectively that the mean and the temporal variability are generally higher in winter than in summer. Then, the amplitude of the daily component α

Y t = a (st) 0 (t) + a (st) 1 (t)Y t-1 + a (st) 2 (t)Y t-2 + σ (st) (t)ǫ t (6)
with {ǫ t } a sequence of independent and identically distributed Gaussian variable with zero mean and unit variance independent of the Markov chain {S t }. Again the daily component is modelled assuming that a Then, if f (t) denotes the value of one of the parameters of the AR models at time t (i.e.

f (t) = a (s) 1 (t), f (t) = a (s) 2 (t), f (t) = σ (s) (t), f (t) = α (s) 0 (t), f (t) = α (s) 1 (t) or f (t) = α (s) 2 (t)
for some s ∈ {1, ..., M} ), we assume a smooth seasonal evolution of the form

f (t) = f 0 + f 1 cos 2π T y t -f 2 (7)
with f 0 , f 1 ≥ 0 and f 2 ∈ [0, 2π[ unknown parameters. Since the conditional standard deviation σ (s) (t) and the amplitude of the daily component α (s) 0 (t) should be positive in order to ensure identifiability, the constraints f 0 > 0 and f 1 < f 0 were added when Due to the complexity of the model and the length of the time-series under consideration, it is important to initialize the EM algorithm with realistic parameter values. Indeed, each iteration of the EM algorithm requires important CPU time and thus the number of iteration needs to be reasonable. Furthermore, using arbitrary values is very likely to lead the algorithm to converge to a non-interesting maximum of the likelihood function.

f (t) = σ (s) (t) or f (t) = α
In practice, the parameters have first been estimated using the least square methods and the parameter values obtained when fitting separately the models to each month of data and then reestimated using the EM algorithm on the whole time series.

We obtain a non-stationary model, which includes both daily and seasonal components and which can be used to generate long wind time series with no discontinuity problems at the beginning of each month. Again, like in Section 3.3, the realism of the simulated time series has been checked by comparing various statistics computed from the synthetic time series to the ones of the original data. We first performed validation on a monthly basis (we considered both calendar month and also periods from the 15th of one month to the 15th of the following month), and the results were similar to those reported in Section 3.3. This is not surprising since, according to Figure 7, the restriction of the fitted seasonal model to monthly time periods are relatively close to the models fitted on monthly data.

We also performed validation on a yearly basis, and in particular we checked the ability of the model to reproduce the interannual variability of the wind conditions. Figure 8 shows that the fitted model underestimates the observed variability in the yearly mean.

Monthly or seasonal validation leads to similar results. This is a well known feature of many stochastic weather generators which is termed "overdispersion" in the literature (see e.g. [START_REF] Katz | Overdispersion phenomenon in stochastic modeling of precipitation[END_REF] and [START_REF] Katz | Mixture model for overdispersion of precipitation[END_REF]). Two possible sources of overdispersion are identified in [START_REF] Katz | Overdispersion phenomenon in stochastic modeling of precipitation[END_REF].

The first one is an inadequate modelling of high-frequency variations and the second one is the presence of low-frequency variations in the climate, on an interannual time scale, which are not taken into account by the model (see also [START_REF] Bouette | Wind in ireland: long memory or seasonal effect?[END_REF]). The results given in Section 3.3 indicate that the model is able to reproduce the short-term dynamics and in particular the autocorrelation function up to time lags of one month (see Figure 5). As a consequence, in absence of interannual components, the variability of the monthly mean should also be well described by the model since for a second-order stationary process the variance of the sample mean can be deduced from the autocovariance function. In the next section we thus investigate the presence of interannual components in the time series under consideration and show that including it into the model help reproducing the interannual variability. It is well known that many long meteorological time series exhibit non-climatic (or artificial) sudden changes due, for example, to an instrument change or a change in station location or exposure and that this may affect the study of the climatic trends. Since we have few information on the existence of such changes for the time series considered in this paper or access to an homogenized version of this time series, the results have been compared with those obtained using reanalysis data. More precisely, we have used the ERA-40 data set which consists in a global reanalysis with 6-hourly data covering the period from 1958 to 2001. This reanalysis was carried out by the European Centre for Medium Range Weather Forecast (ECMWF) and can be freely downloaded and used for scientific purposes at the URL: http://data.ecmwf.int/data

Model with interannual components

The seasonal model described in Section 4 was then fitted to the time series retrieved from the ERA-40 data set for the same location than in situ-data. We obtained generally similar estimation for the parameters of the model, except for the conditional standard deviations σ (s) which are systematically lower for ERA-40 time series: reanalysis data tend to be smoother than in-situ data. Figure 9 also shows the 7-year running mean of the conditional expectation associated to the smoothing probabilities computed using ERA-40 data. We also observe a clear trend and comparing with in-situ data indicates a good overall agreement. For comparison purpose, Figure 9 also shows the Atlantic multidecadal oscillation (AMO) index which is significantly positively correlated with the running means of the smoothing expectation shown on the same figure: periods with higher values of the AMO index seems to coincide with less frequent steady wind conditions. This may be an indication that the observed trend may partly be explained by climatic variations. However, there are also differences between the results obtained using ERA-40 and in-situ data which may correspond to non-climatic breaks in one of these two time series (see also [START_REF] Pirazzoli | Changes in storminess and surges in western france during the last century[END_REF]).

In order to study the impact of the interannual variations on the overdispersion, we have chosen to focus on the time period from 1973 to 1998 when the running mean for in-situ data shown on Figure 9 exhibits a clear increasing trend, and replaced (5) by where, for s ∈ {1...M}, λ s is an unknown parameter which describe possible trends in the probability of occurrence of the regime s. Here, we assume that the long-term climatic variations only impact the probability of occurrence of the different weather types but not the dynamics inside the weather types. Non-homogeneous hidden Markov models, based on similar conditional independence assumptions, have already been proposed in the literature for statistical downscaling (see e.g. [START_REF] Hughes | A class of stochastic models for relating synoptic atmospheric patterns to local hydrologic phenomenon[END_REF] and [START_REF] Vrac | Statistical downscaling of precipitation through nonhomogeneous stochastic weather typing[END_REF]), in which case the hidden weather type is used to link the large-scale circulation to local weather condition.

P (S t = s ′ |S t-1 = s) ∝ q s,s ′ exp κ s ′ cos 2π T y (t -φ s ′ ) + λ s ′ t (8)
Again, the model has been fitted using the EM algorithm. Figure 10 indicates that the fitted model is able to reproduce the observed trend in the probability of occurrence of the different weather types. The estimation for the parameter λ 1 is negative, and it coincides with the fact that regime 1 is less and less likely whereas positive values for λ 2 and λ 3

indicate that the two regimes with more variability become more and more likely.

Figure 8 shows that the model with interannual components better reproduces the variability of the observed mean values compared to the model without interannual component, but still underestimates the observed variability. Using a more sophisticated model for the interannual components could again improve these results. For example, we could replace the linear trend in (8) by a polynomial function, include covariates such as the AMO index, or consider models where the interannual components also modify the dynamics inside the regimes. An alternative in order to improve the description of the interannual variability could consist in using more sophisticated models for the seasonal component (see e.g. [START_REF] Sansom | A hidden seasonal switching model for high-resolution breakpoint rainfall data[END_REF]).

Conclusions

This paper investigates the use of MS-AR models to describe wind time series and it is shown that these models have several virtues. First, thanks to their distributional versatility, they are able to describe the marginal distribution of the time series and thus pre-processing the data, like applying the Box-Cox transformation, is not needed. Then, these models have the ability to model diverse time scales which are present in wind time series and improve the description of important properties of the dynamics such as the durations of calm or stormy conditions. This is an important aspect for many applications of these models. Other related modelling issues which remain to be investigated include the inclusion of the wind direction in the model and the generalization to a space-time model to describe the wind evolution at several locations.

Figure 1 :

 1 Figure 1: Wind speed in ms -1 (y-axis) at Ouessant between 1948 and 1998.

Figure 2 :

 2 Figure 2: Boxplot of the wind speed distribution as a function of the time in the day (x-axis) for January (left) and July (right).

  has M = 2 regimes, but the difference with the model with M = 3 regimes is low. A more precise investigation of these two models shows that the model with M = 3 regimes permits to better reproduce some important properties of the data such as the durations of the storms than the model with M = 2 regimes. The models with M ≥ 4 regimes had states with very low probability of occurrence or the fitted states included two very similar states. This led us to restrict attention to M = 3 and select the model with homogeneous hidden Markov chain and autoregressive models of order p = 2 which has 18 parameters.

Figure 3 :Figure 4 :

 34 Figure 3: Top panel : smoothing probabilities P [S t = s|Y 1 = y 1 , ..., Y N = y N ] for s = 1 (dotted line), s = 2 (dashed line) and s = 3 (full line) for one month of January. Bottom panel : wind speed (y 1 , ..., y N ) for the same month of January. The line style corresponds to the most likely state identified using the Viterbi algorithm with the same convention as on the top panel.

Figure 5 :

 5 Figure 5: Top left: QQ-plot of the marginal distribution for the fitted model (y-axis) against the sample distribution (x-axis). Top right : autocorrelation functions. Bottom left: log-survivor-functions of the time duration of the sojourns below the threshold 6ms -1

Figure 6 :

 6 Figure 6: Left panel : diurnal variations of the mean wind speed. Right panel : periodogram (on x-axis in day -1 ) . The full line corresponds to the sample functions and the dashed line to the fitted model with a 95% prediction intervals dotted line). The distributions for the fitted model were obtained by simulation. Results for July.

  in spring and summer in regime 1 and 2 and at the end of summer in regime 3. The comparison of the values of α (s) 1 and α (s) 0 in the different regimes s ∈ {1, ..., M} shows that the contribution of the daily component to the mean wind speed is more important in the regimes with low temporal variability. Finally, the diagonal coefficients of the transition matrices indicate that the third regime is more persistent in winter than in summer, and thus cyclonic conditions may generally last longer in winter than in summer.

Figure 7

 7 Figure 7 also suggests to let the coefficients of the model evolve smoothly in time instead of using step functions. In this work, we use simple parametric forms to describe the seasonal evolution of the different coefficients. More precisely, we used again (5) for the transition matrix but with T = T y the number of observations in one year in order to obtain a periodic function with period one year. Then we use AR models of order p = 2 with time varying coefficients to model the conditional evolution of the wind speed

Figure 7 : 1 (

 71 Figure 7: Seasonal variations of σ (s) (top left), q s,s (top right), α (s) 0 (bottom left) and α(s) 1 (bottom right). The full line corresponds to regime s = 1, the dashed line to s = 2 and the dotted line to s = 3. The thin line corresponds to the values obtained when fitting separately the models to monthly data whereas the thick line corresponds to the values obtained after fitting the seasonal model on yearly data. The x-axis represents the time in month.

  for s ∈ {1, ..., M}. The model has M(M + 1) + 3M(p + 4) parameters.

Figure 8 :

 8 Figure 8: Distribution function of the annual mean (left) and annual maxima (right). The full line corresponds to the sample function and the dashed line to the fitted models with a 95% prediction intervals (dotted line). The thin lines correspond to the seasonnal model without trend (see Section 4) and the thick lines to the model with trend (see Section 5). The distribution for the fitted model was obtained by simulation. Results for the time period 1973-1998.

Figure 9

 9 Figure 9 shows the 7-year running mean of the conditional expectation E[S t |Y 1 = y 1 , ..., Y T = y T ] associated to the smoothing probabilities for the seasonal model introduced in Section 4. This time series exhibits a clear trend, with low values in the years 1950-1955 and 1970-1975, high values for the years 1960-1965 and a tendency to increase from 1970 to 2000. According to the interpretation of the various states given in Section 3.3, higher expectations may correspond to periods with higher temporal variability in the wind conditions and thus to periods with more frequent cyclonic conditions.

Figure 9 : 7 -

 97 Figure 9: 7-year running mean of the smoothing expectation E[S t |Y 1 = y 1 , ..., Y T = y T ] for the seasonal model together with the AMO index (dotted line). The full line corresponds to the model fitted on in-situ data whereas the dashed line corresponds to the model fitted on ERA40 data at the same location. The three time series have been scaled by removing the mean and dividing by the standard deviation in order to facilitate the comparison.

Figure 10 :

 10 Figure 10: Annual mean of the smoothing expectations E[S t |Y 1 = y 1 , ..., Y T = y T ] for the seasonal model with interannual component (full line) and annual mean of the expectation of the Markov chain with transition probabilities (8) (dashed line) with a 95% prediction intervals (red dotted line). The expectation for the fitted model was obtained by simulation. Results for the time period 1973-1998

  Finally, their interpretability leads to open structure which allows for more physical models. In this work, this is used to include various time scale, from daily to interannual components, in a realistic manner into the model. The model have important limitations. First, it fails in reproducing the lowest part of the marginal distribution and can even simulate negative wind speed. Using alternative autoregressive models or applying an initial transformation on the data may lead to substantial improvements. Further, it fails in reproducing the observed interanual variability and the results obtained in Section 5 suggest looking for a more sophisticated trend model. A trend model which can incorporate scenarios for the future climate would be of particular interest since it would permit to produce simulations consistent with climate change.

Table 1 :

 1 to Table1, BIC clearly favours MS-AR models with autoregressive models of order p = 2 and a number of regimes M between 2 and 4. In January, BIC selects a model with homogeneous hidden Markov chain and homogeneous autoregressive models, that is a homogeneous model with no diurnal component. This seems consistent with Figure2which suggests that the diurnal components are less important in January. In July, when diurnal components are more significant, BIC favours models with homogeneous hidden Markov chain but non-homogeneous autoregressive models. It indicates that it is more appropriate to model diurnal components inside the dynamics of the weather types than in the dynamics of the weather type. This seems also more natural from a physical point of view since the hidden variable is interpreted as a surrogate of the large scale atmospheric situation which may not be affected by diurnal components but only by longer term components (seeSections 4 and 5). BIC values for the various MS-AR models fitted for the months of January and July. The first column indicates if the hidden Markov chain is homogeneous (H) or non-Let us now first focus on January. According to the BIC values given in Table1, the best model

	January

homogeneous (NH), the second column indicates if the AR models are homogeneous (H) or non-homogeneous (NH). The last column gives the number of parameters Npar.

Acknowledgments

The authors are grateful to the three anonymous referees for their constructive comments and suggestions which led to significant improvements in the paper.


M T Boudjiet

V Cuisset

J Bertrand
email: johan.bertrand@andra.fr

I Dufour

PRELIMINARY RESULTS OF THE FEASIBILITY OF HYDROGEN DETECTION BY THE USE OF UNCOATED SILICON MICROCANTILEVER-BASED SENSORS

Hydrogen is a key parameter to monitor radioactive disposal facility such as the envisioned French geological repository for nuclear wastes. The use of microcantilevers as chemical sensors usually involves a sensitive layer whose purpose is to selectively sorb the analyte of interest. The sorbed substance can then be detected by monitoring either the resonant frequency shift (dynamic mode) or the quasi-static deflection (static mode). The objective of this paper is to demonstrate the feasibility of eliminating the need for the sensitive layer in the dynamic mode, thereby increasing the long-term reliability. The microcantilever resonant frequency allows probing the mechanical properties (mass density and viscosity) of the surrounding fluid and, thus, to determine the concentration of a species in a binary gaseous. Promising preliminary work has allowed detecting concentration of 200ppm of hydrogen in air with non-optimized geometry of silicon microcantilever with integrated actuation and read-out.

INTRODUCTION 1.1 Hydrogen release in radioactive waste disposal facility

Many countries have already selected deep geological disposal as the reference solution for the management of nuclear high level and intermediate level long-lived (HL and IL-LL) radioactive waste. After having concluded a feasibility study on deep geological disposal for HL and IL-LL radioactive waste in 2005, the french National Radioactive Waste Management Agency (Andra) was charged by the Planning Act n°2006-739 to study the design and the creation of an industrial center for geological disposal called Cigéo which must be reversible for at least a century-long period. Within the framework of this geological repository project, the observation and surveillance must contribute to acquire the knowledge required to run the disposal and its reversible management.

Hydrogen release is expected in the radioactive waste disposal facility. It originates from (i) radioactive waste release (ii) and anoxic corrosion of metallic materials. In fact, some radioactive wastes (containing α, β, γ radioactivity) resulting from the reprocessing of irradiated fuels are embedded in bitumen matrix. The self-irradiation of the bituminized waste leads mainly to the production of radiolytic hydrogen (75-95 vol. % of produced radiolytic gas). In the deep geological disposal environments steel components would corrode to more oxidized corrosion products and hydrogen gas. The kinetic of anoxic corrosion of Fe components is not well known, it is actually considered an average rate about 10 μm/year. This assumption still awaits experimental verification in the Underground Research Laboratory (sited at Bure in the Meuse district) which aims at studying the feasibility of the reversible geological disposal of high-level and long-lived intermediate-level radioactive waste in the Callovo-Oxfordian clay formation (east of France). Despite the fact that hydrogen releases are expected to be small (in the order of 430 mmol/hour for each intermediate level nuclear waste), when ventilation stops with cell closure, concentrations would slowly and regularly increase. The first calculation gives 4 % hydrogen content in the atmosphere of the radioactive waste cell in less than one year by the end of the oxic corrosion period.

Monitoring of repository structures contributes to security, safety and reversible management of the repository [START_REF] Buschaert | Developing the tools for geologic repository monitoring -Andra's monitoring R&D program[END_REF]. It is important as it relates to the guidance of the disposal process and to the corresponding decision-making process. In the upstream repository-design phases, the hydrogen monitoring system was planned to work under normal operating conditions and to withstand radiation exposure in case of an accidental event. There exist different kinds of hydrogen sensors, including those that are commercially available and others still under development [START_REF] Hübert | [END_REF]. In the present paper, the preliminary results concerning the development of microcantilever-based sensors without sensitive coating for hydrogen detection are presented.

Microcantilever-based sensors without sensitive coating

In recent years, interest in microcantilever-based chemical sensing systems has risen due to their projected high sensitivity [START_REF] Goeders | Microcantilevers: sensing chemical interactions via mechanical motion[END_REF][START_REF] Boisen | Cantilever-like micromechanical sensors[END_REF][START_REF] Li | Integrated microcantilevers for high-resolution sensing and probing[END_REF][START_REF] Sepaniak | Microcantilever transducers: a new approach in sensor technology[END_REF][START_REF] Zhu | Microcantilever sensors in biological and chemical detections[END_REF]. The large ratio of surface area to volume makes the microcantilever extremely sensitive to surface processes. For chemical detection, the microcantilever is usually coated with a chemically sensitive layer whose purpose is to selectively sorb the analyte or molecule of interest. The sorbed substance can then be detected by monitoring either the resonant frequency shift (dynamic mode) or the quasi-static deflection (static mode). A simplified way to explain the basic principle of such sensors is to say that in the case of dynamic mode operation, the change in mass associated with the sorption of analyte onto/into the sensitive coating causes a shift in resonant frequency, which may be correlated to the ambient concentration of the target substance. For staticmode operation, the sorption of analyte causes a cantilever deflection that is induced by surface stress or by the tendency of the sensitive coating to expand or contract upon analyte sorption or desorption (modification of strain and stress in the coating). Typical examples of such microsensors for hydrogen detection are: (i) palladium-coated cantilevers operating in static mode [START_REF] Baselt | Design and performance of a microcantilever-based hydrogen sensor[END_REF] and exploiting the expansion of palladium associated to the formation of hydride compounds during hydrogenation, (ii) and resonator based-microcantilever operating in dynamic mode, using carbone nanotubes (CNTs) as sensitive layer [START_REF] Ono | Mass sensing of adsorbed molecules in subpicogram sample with ultrathin silicon resonator[END_REF] and exploiting hydrogen storage capacity of this material.

The major drawbacks of the use of a sensitive coating come from the fact that they are subjected to environmental effects including temperature, humidity and aging effects. Then, the limited long-term stability of the coatings and the resulting aging affect the reliability of the sensor. To eliminate problems associated with the sensitive coatings, an alternative method to detect and quantify chemical species without the use of a sensitive coating has been proposed [START_REF] Tétin | Modeling and performance of uncoated microcantilever-based chemical sensors[END_REF]. In this paper this principle will be tested for hydrogen detection.

MODELING

The geometry of the cantilever (Fig. 1) is defined by the width b, thickness h, and length L. Coordinate x is measured along the beam length. The properties of the cantilever material are the Young's modulus E and the mass density . The fluid (gas or liquid) properties are defined as  f , the fluid's mass density, and , the fluid's dynamic viscosity.

Resonant frequency shift due to modification of fluid properties

When a microcantilever vibrates in a viscous fluid (gas or liquid), the fluid offers resistance to the motion. The force per unit length, F fluid , which is the consequence of all normal and tangential stresses (hydrodynamic pressure and viscous shear) exerted by the fluid on all the surfaces of the cantilever, can be written in the frequency domain as [START_REF] Blevins | Formulas for natural frequency and mode shape[END_REF]:

    2 1 2 ( , ) , , ( , ) fluid
F x j g x g x w x              ( 1
)
where x is the longitudinal coordinate,  the radial frequency of vibration, w the microcantilever deflection, j is the unit imaginary number, and g 1 and g 2 are functions depending on the fluid properties and the microcantilever cross-sectional geometry. They may be interpreted, respectively, as the distributed damping coefficient of the fluid and the distributed effective fluid mass (per unit length of beam). In general, these quantities are both frequency-dependent and x-dependent; however, in what follows g 1 and g 2 will be assumed uniform in x as this assumption has been shown to be sufficiently accurate in many cases of practical interest [START_REF] Chon | Experimental validation of theoretical models for the frequency response of atomic force microscope cantilever beams immersed in fluids[END_REF].

It is common to normalize F fluid into a dimensionless form called the hydrodynamic function,  =  r ()-j i (), where [START_REF] Sader | Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope[END_REF]:

2 2 4 ( ) ( ) r f g b      (2) 1 2 ( ) 4 ( ) i f g b       (3)
The hydrodynamic function, , represents the total hydrodynamic force per unit length normalized by the force per unit length required to excite (at the same amplitude and frequency) a circular cylindrical volume of diameter b and density  f . The subscripts r and i in Eqs. 2-3 denote the real and imaginary parts of the hydrodynamic function, respectively.

The fluid effects (viscous term g 1 , and inertial term g 2 ) influence the dynamic response of the beam; in particular, the resonant frequency, f r , and the quality factor associated with viscous losses, Q visc , may be expressed as [START_REF] Chon | Experimental validation of theoretical models for the frequency response of atomic force microscope cantilever beams immersed in fluids[END_REF][START_REF] Sader | Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope[END_REF]:

0, 2 2 1 1 1 2 1 r v a c f f Q Lg m    (4) 2 0, 1 2 1 / / visc vac Lg m Q f Lg m    (5)
where f 0,vac is the undamped natural frequency of the microcantilever in vacuum,

    2 2 0, 1.875 2 (12 )  vac f h L E   , Q
is the total quality factor of the cantilever/fluid system (incorporating all losses), m is the total microcantilever mass and L is the microcantilever length.

The hydrodynamic function of a transversely vibrating beam of ribbon-like cross-section (i.e., width b much larger than the thickness h) was previously determined in [START_REF] Tuck | Calculation of unsteady flows due to small motions of cylinders in a viscous fluid[END_REF] by using the method of moments to solve for the velocity field in the fluid and in [START_REF] Sader | Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope[END_REF] by performing a fit of the numerical results of [START_REF] Tuck | Calculation of unsteady flows due to small motions of cylinders in a viscous fluid[END_REF].

A more accurate approximation of the total hydrodynamic force acting on a transversely vibrating microcantilever in liquids should include both the effects of the pressure and shear stress exerted by the fluid on all faces of the beam of finite cross-sectional dimensions. Due to the symmetry of the problem, the hydrodynamic forces acting in the x or y direction on a beam with a rectangular crosssection vibrating transversely will cancel each other out. Only the hydrodynamic forces acting in the z direction will then affect the response of the beam. These forces are the pressure forces acting on the large faces (of dimension b) and the shear forces acting on the small faces (of dimension h). It can be assumed that the microcantilever is long enough so that the hydrodynamic force acting on the surface at the unsupported end of the beam is negligible. Recently, results for the hydrodynamic function in the case of transverse bending that accounted for both the beam's aspect ratio (h/b) and the Reynolds number (Re) of the fluid flow have been published [START_REF] Brumley | Oscillation of cylinders of rectangular cross section immersed in fluid[END_REF]. A Taylor series of the analytical function or a fitting of the numerical results can be made, and then the hydrodynamic function  can be written as

      0 1 2 1 1 Re Re r a b h a b h a b h     (6)
      0 1 2 1 1 Re Re i b b h b b h b b h     (7)
with a 0 , a 1 , a 2 , b 0 , b 1 and b 2 being functions of aspect ratio b/h, and the Reynolds number, Re, defined as:

2 Re 4 f b     (8)
For the ribbon-like case, a 0 , a 1 , a 2 , b 0 , b 1 and b 2 have constant numerical values. In [START_REF] Maali | Hydrodynamics of oscillating atomic force microscopy cantilevers in viscous fluids[END_REF] a fit of the results of [START_REF] Sader | Frequency response of cantilever beams immersed in viscous fluids with applications to the atomic force microscope[END_REF] has been made and the obtained numerical values are:

      (9)
The above equations illustrate the fact that the resonant frequency depends on the fluid properties via the Reynolds number, Re, which appears in the hydrodynamic function  (Eqs. 6-7) and consequently in the terms g 1 and g 2 .

Fluid properties of binary mixture

First, in order to apply the theoretical equations, it is necessary to know the mass density,  f , and viscosity, , of a binary mixture.

Using the definition of the mass density, in the case of a binary mixture, it can be expressed as a function of the mass densities of each component:

  1 1 1 2 1 f x x       (10)
with  1 and  2 the mass densities of gas 1 and gas 2, respectively, and x 1 the volume fraction of gas 1.

Contrary to the case of mass density, the viscosity of a binary mixture is not easy to obtain using the viscosity definition. An approximate expression for the viscosity of gas mixtures was derived by Sutherland in 1895. It is based on a simple mean-free-path model and, for a binary gas mixture, is given by [START_REF] Sutherland | The viscosity of mixed gases[END_REF]:

    1 2 1 1 1 1 12 1 1 21 1 1 1 x x x x x x             (11)
with  1 and  2 being the viscosities of gas 1 and gas 2, respectively.

Concerning the expressions for 12  and 21  , many attempts to predict these parameters have been made. Wilke derived the following expressions [START_REF] Wilke | A viscosity equation for gas mixtures[END_REF]:

2 1 2 1 2 1 4 1 1 2 12 2 2 1 8 1 1 M M M M                                        (12) 2 1 2 1 2 1 4 2 2 1 21 1 1 2 8 1 1 M M M M                                        (13)
where M 1 and M 2 are the molar masses of the gases.

The average deviation between calculated and experimental data for thirteen mixtures involving different gases (CO 2 , O 2 , CO, H 2 , CH 4 , N 2 , etc.) is 1.9% [START_REF] Wilke | A viscosity equation for gas mixtures[END_REF]. Brokaw gave other more complicated expressions resulting in almost the same numerical results [START_REF] Brokaw | Predicting transport properties of dilute gases[END_REF]: for 280 mixtures involving 25 gas pairs (including H 2 /N 2 ), the average error when comparing with experimental data is 0.7% [START_REF] Brokaw | Predicting transport properties of dilute gases[END_REF]. Both of these methods use only the viscosity and molecular mass of the two pure components.

Numerical frequency shift for hydrogen detection

Considering that hydrogen and nitrogen are perfect gases, their mass densities can be calculated for a temperature of 20°C and a pressure of 10 5 Pa. For these temperature and pressure conditions, the mass density of hydrogen and nitrogen are, respectively:

3 2 H 0 082055kg/m   ( ) . ; 3 2 N 1 1488kg/m ( ) .   (14)
The viscosities of the gases also depend on the temperature [START_REF] Gosse | Propriétés de transport des gaz à pression modérée[END_REF]:

6 2 (2,2)* 2.6693 10 MT      (15)
with M the molar mass (g.mol -1 ), T the temperature (K),  a length for the expression of the interaction potential (Å) and (

T T T e e        (16)
with the normalized temperature * T defined by:

* kT T   ( 17
)
where k is the Boltzmann constant and  the energy in the action potential.

Using equations 15-17 and the numerical values of the different parameters for each gas, the viscosities at 20°C for hydrogen and nitrogen are, respectively: Using equations 10, it is possible to calculate the mass density of the binary mixture of 1% of hydrogen in nitrogen and its relative variation compared to pure nitrogen:

3 2 2 1 H 99 N 1 1381kg/m f   ( % / % ) . ; 2 2 2 1 H 99 N 0 928 100 N ( % / % ) . % ( % ) f f d    ( 19
)
Using equation 11, it is possible to calculate the viscosity of the binary mixture of 1% of hydrogen in nitrogen and its relative variation compared to pure nitrogen:

d    ( 20
)
Considering small variations of fluid density and fluid viscosity, it is possible to obtain the analytical expression of the resonant frequency shift due to these fluid property variations. Using equations 4-9, the resonant frequency shift can be approximated with a first-order approximation:

1 0 8 2 2Re f f r r f d df a b d a f h                   (21)
Using this equation, the numerical results obtained in equations 19-20 and the fact that the Reynolds number of microcantilever vibrating in gas is usually high, it comes that the resonant frequency shift is essentially due to the fluid mass density variation. Moreover, from equation 21, one may conclude that the relative frequency shift due to the mass density variation of the gas increases with the ratio of the microcantilever width to the microcantilever thickness (b/h).

To estimate the magnitude of the relative resonant frequency shift due to the increase of 1% of hydrogen in nitrogen, let's consider a silicon cantilever ( = 2330kg/m 3 ) with a ratio of width to thickness (b/h) equal to 25. Using this particular case, a relative resonant frequency shift of 50ppm is obtained. This variation is small but such measurement is proven to be possible. Then test of hydrogen detection has been conducted.

EXPERIMENTS

Examples for the detection of hydrogen in nitrogen using silicon cantilevers not dedicated to this application have been made in order to validate the principle of hydrogen detection with uncoated microcantilever.

Microcantilever design and gas line setup

The microcantilevers have been fabricated by ESIEE group (Paris, France, www.esiee.fr) using silicon technologies [START_REF] Fadel | Chemical sensing: millimeter size resonant microcantilever performance[END_REF].

The excitation of the microcantilever is made by electromagnetic actuation. The excitation arises due to the interaction between a magnetic field created by a magnet placed close to the chip and an electrical current passing through a conducting strip placed on the microcantilever. The Lorentz electromagnetic force induces the structural motion. In order to detect the microcantilever oscillation, semiconductor strain gauges (boron-doped piezoresistors) are fabricated during the process. They are arranged in a half Wheatstone bridge configuration: a first gauge is located where the strains are maximum (the clamped-end of the beam) and the other is on the rigid substrate (Fig. 2). The way to measure the resonant frequency shift is to actuate the microcantilever at a fixed frequency close to the resonant frequency, and to read out the piezoresistive signal phase. When the resonant frequency shifts, the phase of the microcantilever signal is modified. Near the resonant frequency, the phase is almost a linear function of the frequency; thus, the measurement of the phase shift at a fixed frequency allows one to determine the resonant frequency shift.

The microcantilevers are placed in a gas chamber (total volume: 500 µl) under a controlled gas flow (100cm 3 /min). Gas streams containing binary mixtures of the desired species (hydrogen) and nitrogen are fed to the chamber using bottles of gas and a set of mass-flow controllers.

Hydrogen detection

The silicon cantilever used for this test has a length of 2 mm, a width of 400 µm and a thickness of 15 µm. The resonance frequency of the uncoated sensor is approximately 5 kHz.

Figure 3 shows examples of resonant frequency shift measurements for various hydrogen concentrations in nitrogen (2%, 1%, 0.6% and 0.2%). As predicted by Eq. 21, due to the fact that hydrogen is lighter than nitrogen, the resonant frequency increases with the hydrogen concentration. The resonant frequency shift is less than the expected one. This is probably due to the fact that some gas adsorption occurs on the surface of the cantilever.

Nevertheless, other measurements made with different cantilever geometries have proven that, as stated by the theoretical equations, the relative resonant frequency shift is almost proportional to the ratio (width/thickness) and to the variation of the mass density of the hydrogen/nitrogen mixture [START_REF] Cuisset | Utilisation de micropoutres sans couche sensible en vibration pour la détection d'espèces gazeuses. Rélaisation d'un réseau gaz pilotable[END_REF].

In Figure 3, the duration of each concentration step is 3 minutes, this proves that the response time of the sensor without sensitive coating is very small. The 0-90% sensor response time to 1% of H 2 in N 2 is 12 s which is much faster than measurements made with the same gas chamber with polymer coated cantilevers for VOC detection (more than 1 min) [START_REF] Lochon | Silicon made resonant microcantilever: dependence of the chemical sensing performances on the sensitive coating thickness[END_REF] and palladium coated cantilever in static mode for hydrogen detection (90 s) [START_REF] Baselt | Design and performance of a microcantilever-based hydrogen sensor[END_REF].

In order to verify that the measurement is due to a modification of the mass density of the gas, similar test have been conducted with carbon dioxide in nitrogen. The results corresponding to different concentrations (2%, 1.5%, 1% and 0.5%) are presented in Figure 4. In this case, the resonant frequency decreases with carbon dioxide concentration, due to the fact that this gas is heavier than nitrogen. This experiment allows validating that the measured shifts of resonant frequency are induced by the modification of the mass density of the gas. The sensor noise has been quantified by calculating the standard deviation of measurement on one concentration step: it is equal to 0.9 mHz. The estimated limit of detection or the resolution corresponds to a frequency shift three times larger than the noise, which corresponds to 0.02% of H 2 in N 2 (200 ppm), since a resonant frequency shift of 20 mHz has been measured for 0.2% of H 2 in N 2 . In order to validate experimentally this estimated value of the limit of detection with this cantilever, measurements with smaller hydrogen concentrations have been conducted (Fig. 5). It can be seen that, without any signal processing, the detection of 200 ppm of hydrogen is possible and probably the detection of smaller concentration can be achieved. The experimental conditions of the measurement presented in Figure 5 were quite different from the one of Figure 3: the gas flow has been increased from 100mL/min to 250mL/min in order to be able to generate the small concentration of 0.02% of H 2 , and the temperature was approximately 26°C instead of 23°C. These experimental conditions can explain the difference of the resonant frequency shift corresponding to 0.2% of H 2 (0.020Hz and 0.028Hz in Figure 3 and Figure 5, respectively). The sensitivity of such sensors depends on the temperature because both the gas density and the cantilever Young's modulus are temperature dependent. Then, in order to have repeatability in the measurement the temperature has to be the same (this was the case for the two resonant frequency shifts corresponding to 0.2% of H 2 in Figure 5).

CONCLUSION AND FUTURE WORKS

Compared to classical resonant microcantilever-based chemical sensors with sensitive coatings, the uncoated microsensors exhibit shorter response times because there is no analyte sorption into the coating affecting the sensor response. Moreover, the absence of the sensitive coating leads to a more reliable and reversible behavior because there is no significant absorption and desorption phenomena.

The above results indicate that uncoated microsensors may serve as viable devices for the detection of specific concentrations of one gas in a binary mixture. The sensitivity and resolution of such sensors will be larger for those cases in which the difference between the mass densities of the two gases is higher (which is the case for hydrogen in nitrogen or air). The major drawbacks of such sensors is that there is no intrinsic selectivity; thus, these sensors can only be used for specific applications, such as when monitoring environments in which it is known that only one gas concentration can vary, while the relative concentrations of the other gases remain fixed. Another intrinsic drawback concerns the sensitivity which is smaller than coated microcantilevers, although, as presented in this paper, the sensitivity is sufficient to detect 200ppm of hydrogen in nitrogen.

These preliminary results are promising and the future work will consist to optimize the geometry of the vibrating microstructure (increase of the width to thickness ratio, paddle-shaped geometry, etc.) in order to increase the sensitivity, decrease the noise and then improve the limit of detection. Another study will concern the influence of parasitic variations (temperature, third gas, etc.) and the compensation of these variations.

Figure 1 .

 1 Figure 1. Schematic representation of the microcantilever geometry and of the transverse bending vibration

Figure 2 .

 2 Figure 2. Schematic representation of the microcantilever electromagnetic actuation and piezoresistive read-out

Figure 3 .

 3 Figure 3. Example of detection of hydrogen in nitrogen with an uncoated 2000 x 400 x 15 µm 3 silicon cantilever

Figure 4 .

 4 Figure 4. Example of detection of carbon dioxide in nitrogen with an uncoated 2000 x 400 x 15 µm 3 silicon cantilever

Figure 5 .

 5 Figure 5. Example of detection of small concentration of hydrogen in nitrogen with an uncoated 2000 x 400 x 15 µm 3 silicon cantilever