David Ginsbourger
email: ginsbourger@emse.fr

Delphine Dupuy

Anca Badea

Laurent Carraro

Olivier Roustant

A note on the choice and the estimation of Kriging models for the analysis of deterministic computer experiments

Keywords: Metamodeling, Kriging, Maximum Likelihood, Determinisitic Drift, Additive Models

Our goal in the present work is to give an insight on some important questions to be asked when choosing a Kriging model for the analysis of numerical experiments. We are especially concerned about the cases where the size of the design of experiments is small relatively to the algebraic dimension of the inputs. We first fix the notations and recall some basic properties of Kriging. Then we expose two experimental studies on subjects that are often skipped in the field of computer simulation analysis: the lack of reliability of likelihood maximization with few data, and the consequences of a trend misspecification. We finally propose an example from a porous media application, with the introduction of an original Kriging method in which a non-linear additive model is used as external trend.

Linear predictors for spatial interpolation of numerical simulators

We study a deterministic numerical simulator as a function z : D ⊂ R d → R, where x ∈ D is the vector of inputs variables. We denote the set of the design points (or "design") by X = {x 1 , ..., x n } and by Z = {z(x 1 ), ..., z(x n )} the set of simulator responses associated with X. Kriging is a class of methods coming from the field of geostatistics [START_REF] Matheron | Principles of geostatistics[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF], known as linear optimal prediction in classical statistics. It provides at each point x ∈ D a prediction Ẑ(x) linearly depending on Z, where the weights depend on the design and on the Kriging model but not on the observations. The way the weights are defined varies as a function of the type of Kriging -Simple (SK), Ordinary (OK), Universal (UK), etc-and many parameters such as the trend functions, the covariance kernel and their own parameters: threshold (or "sill"), scales, nugget, etc... denoted by the r-dimensional vector ψ. In the following, we will concentrate on the parameters of sill and scale (r = 2), denoted respectively either by ψ 1 , ψ 2 or by σ 2 , p ∈ [0, +∞[. Most classic Kriging types (including SK, OK, UK, and more) can be interpreted as random process interpolation relying on the assumption that:

∀x ∈ D, z(x) = t(x) + ε(x) (1)
where t is a numerical deterministic function and ε(x) is one path of a centered stationary Gaussian Process (GP) with known stationary covariance kernel k : h ∈ R d -→ k(h) ∈ R. t is generally known up to a set of parameters or a semi-parametric structure to be estimated within Kriging. Several founder works [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] on the application of Kriging to computer simulations start off with an extremely simplified version of (eq.1). They assume that the trend is an unknown constant (Ordinary Kriging, i.e. t(x) = µ ∈ R) and that k is a generalized exponential kernel [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF], letting the stochastic part of (eq.1) account for the variability of z. Then the covariance parameters ψ are estimated by maximizing the Gaussian likelihood of the observations Z.

On the other hand, recent approaches [START_REF] Jourdan | Approches statistiques des expériences simulées[END_REF][START_REF] Martin | Use of kriging models to approximate deterministic computer models[END_REF] try to take advantage of more complex trends, from linear and polynomial functions to Fourier series. In other respects, [START_REF] Martin | A monte carlo simulation of the kriging model[END_REF] as well as [START_REF] O'hagan | Bayesian analysis of computer code outputs: a tutorial[END_REF] present an application of bayesian analysis to Kriging interpolation of computer codes.

The motivation of this article is to raise some basic questions that should become crucial when applying Kriging techniques with few observations regarding the dimension of inputs, which is quite often the case in numerical simulation. The two coming sections, based on toy experiments, put a focus on the estimation of the covariance parameters ψ and on the choice of the trend t. The two following sections are dedicated at presenting an original combination of additive models and Simple Kriging, with a heuristic fitting methodology. The efficiency of this technique is illustrated on a 3-dimensional example from a porous media simulation test case.

Fitting covariance parameters by MLE with a small sample

The Maximum Likelihood (ML) estimation method is widely used in Kriging to choose covariance parameters on the basis of observations. Following the assumptions from (eq. 1), ML estimation relies on the maximization of the density of the observed values Z, seen as a function of the vector ψ:

L(ψ; Z) := f (Z|ψ) = (2π) -n 2 det(K ψ ) -1 2 e -1 2 (Z-t) T K -1 ψ (Z-t) (2)
where K ψ is the covariance matrix of Z(X) = {Z(x 1 ), ..., Z(x n )} provided that ψ is the true vector of covariance parameters, and t is the vector of values of t at X. The obtained result ψ = argmax ψ {L(ψ; Z)} is closely depending on Z, i.e. on the observed realization of Z(X). The behaviour of ψ relatively to ψ when the sample of observations fluctuates is a of importance. We recall that Z is assumed to be one realization of a multivariate Gaussian random vector with given trend, covariance structure, and covariance parameters ψ. Then L(.; Z) becomes a random function (fig. 1), and ψ = argmax ψ {L(ψ; Z)} becomes a random vector as well. The distribution of ψ has been studied in detail within the theory of likelihood [START_REF] Sweeting | Uniform asymptotic normality of the maximum likelihood estimator[END_REF][START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF]. A first order Taylor expansion leads to an asymptotic result [START_REF] Abt | Fisher information and maximumlikelihood estimation of covariance parameters in gaussian stochastic processes[END_REF] based on Fisher's Information Matrix I(ψ) (denoted by FIM in the sequel):

       ψ L -→ N (ψ, I(ψ) -1 ) I(ψ) = E ∂ ln(L(.;Z)) ∂ψ i (ψ) ∂ ln(L(.;Z)) ∂ψ j (ψ) i,j∈[1,r] (3)
In many computer experiments, one first picks a covariance kernel from a parametric family: Gaussian, Exponential, Matèrn, etc...(the Gaussian covariance kernel is often chosen for its simplicity and regularity properties) and the associated covariance parameters are then automatically fitted by ML. However, the efficiency and robustness of this estimation method when few data are available are rarely discussed. Our concern is to check in what measure the first order asymptotic results hold with small samples. To do so, we computed the theoretical FIM of Z:

∀i, j ∈ [1, r] (I(ψ)) ij = 1 2 tr K -1 ψ ∂K (.) ∂ψ i (ψ)K -1 ψ ∂K (.) ∂ψ j (ψ) (4)
To obtain comparable results for different values of ψ, we introduce a relative inverse FIM: (J (ψ)) ij = (I -1 (ψ)) ij /(ψ i ψ j ). J is in fact the asymptotical covariance matrix of ψ ψ , where the division is made component by component. We conduct experiments with vectors taken from simulated monodimensional Gaussian Processes to compute empirical means and variances of the ML estimators. For each simulation, we compute covariance parameters estimated by ML and the Integrated Squared Error (ISE) between simulated (z(x)) and interpolated ( Z(x)) data:

ISE = 1 vol(D) D |z(x) -Z(x)| 2 dx (5)
where vol(D) is Lebesgue's measure of the set D. ISE is approximated by averaging the squared errors on a fine grid (i.e. 200 points). We finally collect the averages and variance matrices of the relative values of the estimated covariance parameters, the averages and variances of ISE (ISE is random since it depends on the realization z), and the covariances between ISE and

ψ rel i = ψ i -ψi ψ i .
The latter two indicators are not presented in the tables. We focus here on GPs with covariance kernels c g (h) = σ 2 e -h 2 p 2 (Gaussian) and c e (h) = σ 2 e -|h| p (Exponential). The covariance parameters reduce to ψ = (σ 2 , p) ∈]0, +∞[×]0, +∞[ and the design X is chosen among uniform subdivisions of [-1, 1]:

X n = {-1, -1+ 2 n-1 , . . . , -1+ 2(n-2) n-1 , 1} (n ∈ N\{0, 1}
). We restrict our experiments to the designs X 5 and X 10 with both c g and c e , and covariance parameters ψ 1 = σ 2 ∈ {5, 10}, and ψ 2 = p ∈ {0.3, 0.4, 0.5, 0.6}.

= ψ i -ψi ψ i , i = 1
, 2 and for X = X 5 . The second column shows that the relative ML estimates are almost unbiased even with 5 observations. On the contrary, a comparison between the third and fourth columns illustrates that the ψ rel i are clearly more dispersed than given by the asymptotical approximation based on the FIM. In the case of a Gaussian covariance, we observe a negative relative bias1 (-3.4% to -25.6%) in the estimation of ψ 1 = σ 2 . This bias is decreasing with the number of design points #X (see table 2 where the negative relative bias varies between -4.2% and -7.5%), which seems in accordance with the asymptotic unbiasedness of MLE. On the other hand, the relative bias of ψ2 has a small order of magnitude when #X = 5 and slightly oscillates around 0 when #X = 10.

The empirical covariance matrices of the ML estimates offer some surprising results. In particular, the relative variances of ψ1 present huge fluctuations: they vary sometimes of an order of more than 10 between two samples of 1000 realizations issued from the same GP; for instance by resimulating a GP with . Since it is in contradiction with normality and the order of magnitude given by (eq.4), we shall analyze this phenomenon in detail. First, we observe that the extreme values of V ar[ ψ1 ] are caused by some outliers, highly perturbating the non-robust estimate of variance. Second, the histogram in (fig. 2) illustrates that the distribution of the ψ1 's is rather lognormal than normal. Finally, the comparison with the relative FIM shows that the empirical variance of ψ1 is clearly larger than predicted by the second order Fisher approximation, in particular with the smallest designs.

Concerning the relative variances of ψ2 , the results are much more regular: they decrease monotonically with ψ 2 and with #X, both for the empirical and theoretical quantities. Once again, the empirical variances tend to match the theoretical variances as #X grows, even if the first ones are still typically two times larger than the second ones for a sample of size 10. In other respects, both tables illustrate some fundamental properties of the mean squared error. Obviously decreasing with #X, the ISE is also decreasing with the range ψ 2 and linearly increasing with the variance ψ 1 . Finally, we quantify the linear dependence between the underestimation of both covariance parameters by MLE and the ISE (not in the tables). It is worth noticing that ψ 1 and ψ 2 play drastically different roles here: it seems that a bad estimation of ψ 1 is weakly correlated with the ISE. This result seems natural when considering that the OK predictor is not depending on the process variance, see [START_REF] Cressie | Statistics for spatial data[END_REF]. Conversely, the correlation between the ISE and the relative MLE error on ψ 2 is significantly positive: it varies between 40.1% and 55.7% when #X = 5 and between 15% and 62.5% when #X = 10. This coincides with our previous qualitatives observations of larger ISE when the range is much underestimated. and the asymptotic law (black lines) for the scale parameter for increasing size of the design. The boxplots for the experimental laws have been done using 1000 simulations, with Gaussian covariance function of parameters ψ = (5, 0.5). For the asymptotic law the median is represented in continuous line and the first and third quartiles in dashed lines. Right: Histogram of the logarithm of relative errors obtained when estimating ψ 1 by ML using 1000 GP realizations. The shape of this histogram suggests that the distribution of relative errors is far closer to a lognormal law than to a Gaussian.

A similar study with exponential covariance function gives very different results both for the bias and the variances of ML estimates (the corresponding tables are not presented here). Indeed, we observe very regular variances of ML estimates while the bias reaches impressive orders of magnitude. However, the behaviour of the ISE and the correlations between ISE and relative MLE errors follow the same sketch as in the Gaussian case.

To sum up this section about ML estimation:

• Fisher's first order asymptotical results must be applied with much care concerning the sample size. More precisely, it has been observed here for n ≤ 5 that the distribution of the estimated range parameter is asymmetrical with a higher variance than the inverse of Fisher's information, but quickly stabilizes to a Gaussian when n increases (from 5 to 13).

• On the other hand, the distribution of the estimated variance parameter has a very large right tail but its shape is far from being gaussian when n is very small (n ≤ 5). Furthermore, these results still hold when n increases (from 5 to 13) and it seems that the Gaussian approximation becomes reasonable only for larger values of n.

Estimating covariance parameters by ML with few data appears to produce very dispersed results. Hence, it seems unreasonable to neglect the uncertainty associated with this phase of estimation when performing Kriging.

Bayesian techniques are a way to address this issue [START_REF] O'hagan | Bayesian analysis of computer code outputs: a tutorial[END_REF][START_REF] Goria | Evaluation d'un projet minier: approche bayésienne et options réelles[END_REF]. In other respects, [START_REF] Abt | Estimating the prediction mean squared error in gaussian stochastic processes with exponential covariance structure[END_REF] investigates an extended Kriging variance taking the estimation of parameters into account; a this stage the latter relies on the first order approximation. To finish with, frequentist approaches based on the maximization of penalized likelihood functions seem very promizing since they provide estimators with the same asymptotic properties as ML in addition to a more robust behaviour with few observations [START_REF] Li | Analysis of computer experiments using penalized likelihood in gaussian kriging models[END_REF].

3 Kriging with trends: a bless or a curse?

Now we wish to examine another difficulty encountered when Kriging based on few data: the selection and the estimation of deterministic trends. In computer experiments, the most commonly used Kriging model seems to be Ordinary Kriging. However, OK reaches one of its limits when the stationarity assumption does not hold any longer, i.e. when non constant trends t(x) are impossible to ignore. In this case, we are back to the general decomposition of (eq.1), where z is assumed to be the sum of a deterministic trend t and one realization of a centered GP ε. At this stage, we may consider several subcases. If t is known and the parameters of ε have to be estimated, a straightforward solution is to perform Simple Kriging of the residuals {z(x)t(x)} x∈D . If t is unknown, it is common to distinguish between a linear and a more general non-linear framework. The case in which t depends linearly on its parameters and ε has a kwown covariance structure has been intensively studied: it is well known as Universal Kriging [START_REF] Martin | Use of kriging models to approximate deterministic computer models[END_REF]. When the covariance parameters ψ are known and the trend is a linear combination of some chosen basis functions

f j (j ∈ [1, b], b ∈ N\{0}), the only unknowns are the parameters of the trend (∀j ∈ [1, b], β j ∈ R); indeed, if t(x) = b j=1 β j f j (x)
, the β j 's can directly be estimated by Generalized Least Squares (GLS):

β(ψ 2 ) = (F T K -1 ψ F) -1 F T K -1 ψ Z = (F T R -1 ψ 2 F) -1 F T R -1 ψ 2 Z (6)
where F denotes the evaluation of f (x) = [f 1 (x), . . . , f b (x)] at the n design points and R ψ 2 = (1/ψ 1 )K ψ (proportionality since the observations are noisefree) is the correlation matrix of Z(X).

In practice, however, one has seldom the value of the covariance parameters at disposal previous to performing UK. So one has to estimate a model with linear trend and unknown covariance parameters ψ (in the following we will also refer to this case as "UK", like many practitioners do). Hence ψ and β have to be estimated within Kriging. At a first sight, this is likely to create a circularity problem: one needs a known trend to work on the residuals and thus estimate ψ. On the other hand, estimating t without taking the residuals into account may lead to unadapted trends (the estimation of the trend parameters would rely on Ordinary Least Squares instead of GLS).

Fortunately, ML estimation gives a way to escape this vicious circle. Assuming, like in section 2 that the covariance parameters to be estimated are ψ = (σ 2 , p), and using MLE (and the same formula (6) for β), one can get a straightforward formula for σ 2 , explicitly depending on ψ 2 :

σ 2 (ψ 2 ) = (1/n)(Z -F β(ψ 2 )) T R -1 ψ 2 (Z -F β(ψ 2 )) (7)
By injecting ( 6) and [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] in the expression of the likelihood, one can obtain a concentrated likelihood function L(ψ 2 , σ 2 (ψ 2 ), β(ψ 2 )) which clearly depends only on ψ 2 and which has to be maximized to get ψ2 . The Kriging predictor with plugged-in covariance parameters is then given by:

Ẑ ψ2 (x) = f T (x) β( ψ2 ) + r T (x)R -1 ψ2 Z -F β( ψ2 ) (8)
where r(x) is a vector of correlation values between Z at an unknown point x and at the points of the design X. Most of the time (apart in Bayesian Kriging) the variability due to the estimation ψ 2 is not propagated, and one uses the regular UK prediction variance.

UK appears as a very convenient means to incorporate known deterministic trends within Kriging. By the way, we will see in the next section that overcoming the circularity problem is not easy in a more general non-linear framework. Now we would like to go one step deeper in practical considerations and raise a naive but complex question which has to be handled in real-world applications, and particularly in high-dimensional problems: how can one come back to the nature of the trend from raw data? As soon as neither prior information nor obvious graphical clue is available, one has indeed to select a trend on the basis of (X, Z). What means does he have to do so, and what risk does he run in case of a bad choice? In order to show that these questions are crucial, let us first perform some toy experiments. The set-up is the following. A realization of a one-dimensional GP with known covariance function and parameters is simulated on a regular grid (401 points on [-1, 1]) and an affine trend is added; From this set we choose different subsets of points and perform three types of Kriging : OK, UK with linear trend and UK with quadratic trend. The Kriging predictor fails to capture the shape of the realization at the center of the domain. The third design is made of three points clustered at the center of the domain; OK automatically comes back to the mean value outside of the design and dramatically miss the actual trend.

We choose at first a subset of 5 regularly distributed points. Due to the fact that the points are regularly spaced on the grid, all the three kriging give similar good results, even if in two of the three cases the trend is misspecified (fig. 3 left for the case of Ordinary Kriging). This may lead to the conclusion that specifying the trend is not very important and we could obtain good results using OK. But if we perform the same Krigings on different designs, where there are few points concentrated either on the boundaries or in the center of the domain, then the results are very bad (due to the ratio between the parameter ψ 2 and the subdivision length) when the trend is misspecified, see (fig. 3, middle and right). The covariance parameters used for the simulated process in (fig. 3) and (fig. 4) are ψ = (5, 0.2). The results are even ). The GP is here represented between -2 and 2, such that all three cases can be referred to as extrapolations. UK with affine trend (center) gives accurate results on the whole domain. On the contrary, both OK (left) and UK with quadratic trend (right) give good results between -1 1 but dramatically fail in extrapolation.

worse if we use the Kriging predictor given by OK or by UK with quadratic trend in extrapolation (fig. 4, left and right). In the one dimensional case the choice of the trend doesn't seem to be essential while interpolating data which are not very distant one from another with respect to the frequency of variation of the process. On the contrary, when the design is not regular and we are in extrapolation, the performances of Kriging are very sensitive to the adequacy between the real trend of the process and the Kriging trend.

Hence it seems enough to properly fill the space to avoid the risks caused by the choice of trend functions. But what is possible in one or two dimensions becomes unrealistic when the dimension increases: a design with only one point at each vertex of a cubic domain [0, 1] d has 2 d points, i.e. 1024 points in 10 dimensions and more than a billion points in 30 dimensions. As we usually dispose of 10 × d observations per dimension, which is already an optimistic case, choosing a trend based on data only appears as a very difficult task. Let us see nevertheless what would be possible in order to choose a trend starting from a data set (X, Z): the classical frame of linear regression offers a panel of diagnostic tools dedicated to validating both assumptions on the trends and on the model of residuals. For instance, commonly used indicators include R 2 (and R 2 adjusted), the F-ratio and the p-values for each estimated regression coefficients, and numerous criteria to check the adequacy of the residuals to the underlying model. In most cases the Gaussian likelihood of the residuals is considered among the relevant criteria of model selection (some model testing techniques or even based upon it).

Now it seems necessary to recall that the latter measures are exclusively done at the design of experiments, also called "training sample" or "learning sets" in the literature of statistical learning, see [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Selecting only on the basis of a R 2 fit would lead for instance to the systematic choice of models interpolating (X, Z). However such models are not meant to be good in prediction outside the design of experiments. This warning leads to the double message:

• Model complexity must be taken into account in selection procedures

• Testing the model at some test points not used in the model fitting could be worth: this is for instance what cross-validation does.

The following experiment is performed in an intent to illustrate the first point.

The second point will be illustrated in the next section. Here we investigate on a simple case how trend selection may be misleading when likelihood is the only criterion, without any consideration of model complexity. To do so, we compute, for each trend form of the Kriging model, the optimal parameter p by ML, we compare the corresponding values of the likelihoods and we select the kriging model having the highest value of likelihood. In table 3, we compare three Kriging models (OK, UK with linear trend, UK with quadratic trend) for three different functions: one realization of a one-dimensional GP with 11 points and with Gaussian covariance function (ψ = (5, 0.4)), the same realization plus a linear trend 0.5 + 5x, and the same realization plus a quadratic trend 0.5 + 5x + 5x 2 .

Here it is essential to notice that the likelihood values are necessarily larger when adding more degrees of freedom to a statistical model. This constitutes a misleading incentive to always choose the model with the largest number of parameters within a given family. This happens for instance between Kriging models with first order and second order polynomial trends. As can be observed in table 3, L always increases (i.e. the values of -2 ln(L(p)) will decrease) with the complexity of nested model. What we should really compare are maximum likelihood values between models with the same number of degrees of freedom. On the last line of table 3, in the cases of the GP without trend and of the GP with linear trend , the estimated values β are very close to but different from zero. Thus the model obtained by automatically selecting the Kriging with highest likelihood will perform badly in extrapolation because of the higher order terms of the polynomial. The same phenomenon applies in an even more pronounced way with a linear trend in the case of a centered GP (first column, second row).

As a conclusion to this section, we have pointed out that OK and UK may seem to deliver similar results when the design is dense [START_REF] Michael | Interpolation of Spatial Data, Some Theory for Kriging[END_REF], but modeling the trend matters in extrapolation situations [START_REF] Journel | When do we need a trend model in kriging?[END_REF]. Since working in highdimensional spaces means that we will practically always be in extrapolation, we need exploratory and visualization tools dedicated at finding trends in multivariate data. Recent methods of data mining and functional analysis may help [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. We propose now to use additive models within spatial interpolation.

Using non-linear additive models as external drift

Linear models are often used by practitioners of quantitative disciplines since they are simple to interprete and to assess. Additive models (AM ) are an extension of linear models. A precise description of these models can be found for instance in the book [START_REF] Hastie | Generalized Additive Models[END_REF]. The advantage of AM is to conserve the feature of non-interacting predictors, but they allow much more flexible inference for each univariate problem, using kernel smoothers for instance [START_REF] Wahba | Spline Models for Observational Data[END_REF]. The generic expression for an additive model is the following:

   Z i = z(x) + ε i z(x) = α + d j=1 f j (x j ) The ε i are n.i.i.d. ( 9
)
and the f j s are arbitrary univariate functions, one for each predictor but possibly not the same kind of function for each dimension. Hence additive models deal with additive functions observed in a Gaussian noise. x is in fact assumed here to be a control variable, and not a random variable as in [START_REF] Hastie | Generalized Additive Models[END_REF]. This model may be used to approximate deterministic computer experiments, provided that the response surface can reasonably be decomposed in an additive way. Once the nature of the f j s is chosen, they can be estimated using a powerful iterative procedure called backfitting algorithm, see [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. Backfitting means that f 1 is estimated on the basis of all data (X, Z), then f 2 is fitted to the residuals Zf 1 (X), and so on. Under mild assumptions, the backfitting algorithm converges and finds the unique solution of the additive decomposition of (eq.10). In this section, we propose a combination of additive model and Kriging that offers the great flexibility of AMs and yet interpolates the data. It seems very natural to combine both models by using the following decomposition:

   z(x) = t(x) + ε SK (x) t(x) = α + d j=1 f j (x j ) ε SK (x)
is a GP realization like in (eq.1) (10) This identity may at first seem similar to the equation of Universal Kriging. However, in this case the non-linear nature of the trend prevents one from solving the estimation globally. Indeed, a likelihood maximization would lead to an optimization problem in infinite dimension:

max ψ,(f j ) j∈[1,d] L(ψ, t; Z) (11)
To our knowledge such a problem is analytically intractable. On the other hand, the backfitting algorithm is not suited anymore if we take the Kriging part into account. Indeed, kriging the residuals after fitting a smoother in one dimension would lead to an interpolation and thus end the iterative procedure without fitting the additive parts in the other dimensions.

Kriging with external trend [START_REF] Cressie | Statistics for spatial data[END_REF] seems to constitute a good alternative for solving both the problem of the "general" form of trend and the one of the circularity. Consequently, we consider now a two-step approach (see Alg.1): first, the additive trend t(x) is estimated using the backfitting algorithm, and then Simple Kriging is applied to (Zt) with covariance parameters estimated on the basis of those residuals, by likelihood maximization or other.

Algorithm 1 A first two-step approach to fit a Kriging with additive trend 1: Estimate the trend t by backfitting 2: Estimate the covariance parameters and fit a SK model on the basis of the residuals at X

Unfortunately, there are significant drawbacks in the latter procedure, mainly related to the uncontrolled trade-off between deterministic and stochastic parts. Hence, the whole uncertainty reduces here to the Kriging variance estimated on the residuals; there is indeed no global uncertainty on the trend unless we use only splines in the AM. This is likely to cause a large underestimation of the process variance associated with the model. Furthermore, these residuals may be not very well suited to estimate the Gaussian process part: the additive model is constructed to fit z accurately at the design -possibly leading to overfitting-, thus the residuals at X are likely to vary with a smaller magnitude than in prediction. Since we look for a model with reasonable generalization properties, it seems necessary to find an alternative way of estimating the covariance parameters.

We propose here a sequential estimation technique for combined Kriging models like (eq.10). It is based on the idea that when the trend is non-linear, the parameters of the GP model should be estimated on a validation set rather than on the set at which the trend is fitted.

Algorithm 2 An alternative two-step approach to fit a Kriging with additive trend 1: Consider two designs X 1 and X 2 ⊲ possibly obtained by splitting X 2: Estimate the trend t by backfitting, based on the data (X 1 , z(X 1 )) 3: Estimate the SK covariance parameters on the basis of the residuals {t(x)z(x)} x∈X 2 4: the SK model on the basis of all residuals ⊲ with parameters estimated at the previous step

A 3-dimensional application of Kriging with

Additive Trend (KAT)

The previous approach is applied to a 3-dimensional example from an industrial test case. The data are obtained with a flow simulator and the numerical response z, standing for the outcome of interest, is studied as a function of three physical parameters characterizing the porous media and denoted by x 1 , x 2 and x 3 ∈ [-1, 1]. The response is simulated at 1331 locations corresponding to a 11-level full factorial design, denoted by "F" in the sequel. Our goal is to provide a surrogate of the simulator on the basis of a poor design of experiments. The metamodel should interpolate the data (to respect the determinism of the underlying simulation) and provide a prediction uncertainty that allows statistical-based exploration, for instance to solve optimization problems. Furthermore, it should take into account a prior knowledge inherited from a previous study: the phenomenon is almost additive in its parameters.

Our initial design, "X 1 ", is a 20-elements Hammersley sequence. We first perform a graphical analysis (fig. 5) of the response at X 1 , discuss the hypothese of additivity, and propose several kinds of linear and additive trends to model the data. Algorithm 1 is tested with the design X 1 (Table 4). Then a second design, "X 2 ", is used for an intermediate validation of the covariance parameters of the model previously obtained. X 2 is made of 14 points taken from a 40-elements D-optimal design (see fig. 6). Algorithm 2 is then performed by re-estimating the covariance parameters of the previous SK model on the basis of the residuals at X 2 . An original estimation method is proposed, which differs from the traditional MLE: the process variance σ 2 is fixed such that the standardized residuals have most of their values between -2 and 2 [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and the range parameter p is chosen in order to minimize the ISE at the design X 2 (fig. 7). The full factorial design F is finally used for a phase of model validation (fig. 8).

A graphical analysis of the coplots at X 1 does not reject the prior belief of additivity. A first additive decomposition is then estimated using splines in all directions (referred to as "GAM splines" in the following). We observe that we might take a linear trend in the directions of x 1 and x 3 , and a non-linear trend in x 2 without loosing much accuracy (see Table 4 for a quantitative validation). Hence we choose to fit an additive model with mixed trends, called "GAM mixed" in the sequel.

Different Krigings with external trend are fitted to the observed data at the design X 1 . In all cases, the SK part has a structure of isotropic GP with Gaussian covariance (see section 2). We focus on the two additive trends defined above and on two additional linear trends: a first and a second order regression polynomial. For each model, we fit the trends respectively by OLS and backfitting, and we measure their relevance using indicators computed with the residuals at the design X 1 (residuals deviance and p-values when available). Then we fit a Kriging to the residuals, as explained in Algorithm 1. For each Kriging, we store the maximum reached value of the log-likelihood and the corresponding range and variance values. The results are listed in (Table 4).

These results support the belief that a general additive trend is adapted for these data: both the variance of residuals and the values of their likelihood (compared to the 2nd order linear model, which uses more degrees of freedom) indicate their good fit to the data.

In practice, however, we care more about the model's abilities to make correct predictions at new points than about its mean squared error at the design. Hence, model validation should not be blindly supported by the indicator R 2 or the likelihood of the residuals at X 1 . First, we should consider the number of degrees of freedom of the model. Second, it may be worth validating the model outside of the design. Indeed, the residuals drawn from (fig. 5) are computed at the same locations as those used to fit the trend.

Concerning the first point, we can compare the degrees of freedom of both "2nd order linear" and "GAM mixed": respectively 10 and 7. Regarding the second point, we conduct a validation test on some additional data, inspired by the cross-validation procedure. Following Algorithm 2, X 2 is used to valid and update the parameters associated with the model fitted at X 1 . Coplots representing both X 1 (dots) and X 2 (triangles) designs in projection on all pairs of coordinates. The three graphics illustrate the spacefilling behaviour of X 1 and the D-optimal nature of X 2 . X 2 also appears to be reasonably disconnected from X 1 .

The points of X 2 are used to test the validity of the covariance parameters of the model "GAM mixed", previously estimated by ML. Figure 7 shows the associated residuals standardized by the ML variance (left), and the behaviour of the ISE at X 2 as a function of p (right). We recall that the residuals should satisfy the assumption of normality in order to get relevant Kriging variances, insuring correct statistical predictions. Figure 7 (right) shows that the ISE at the validation design X 2 can be significantly reduced by increasing the range p. Following Algorithm 2, we reestimate the covariance parameters based on these residuals at X 2 . Instead of using ML however , we prefer to directly use the work done hereabove to compute the ISE as a function of the range. It appears indeed that the optimal range to accurately fit the residuals at the validation design is given by p 2 = 1.4. Concerning the variance, we observe more satisfying standardized residuals with σ 2 2 = (2 × σ M L ) 2 . So we keep σ 2 2 .

Remark: A ML estimation with the residuals at X 2 delivers p = 0.97.

We finally test the model of Algorithm 2 at the design F (fig. 8). The standardized residuals (with the variance σ 2 2 ) and the ISE as a function of p validate our empirical decisions made on the basis of the intermediate design X 2 (note that ML on X 2 -see remark hereabove-gives also better results than ML on X 1 but the cross-validating strategy minimizing the ISE at X 2 remained the best). To conclude with, the algorithm investigated performed well on this example: Simple Kriging seems to constitute a good complement to additive models in an intent to interpolate data and also possibly explain some non-additive part. The method we use here allows inference of covariance parameters with values suited for a correct quantification of uncertainty. This seems encouraging to develop further "cross-validation-like" methods for the combination Additive model + Kriging.

Conclusions and perspectives

We observed in a one-dimensional frame that MLE could behave very differently from Fisher asymptotical results when n is small. This result should be kept in mind when dealing with higher dimensions, and further studies have to be done in this latter context. Since it relies on the simulation of Gaussian vectors, the experimental approach presented here can easily be transposed in a higher dimensional framework. Perspectives include the empirical comparison of ML and penalized ML [START_REF] Li | Analysis of computer experiments using penalized likelihood in gaussian kriging models[END_REF] when using classical designs.

Further experiments on the topic of trend selection illustrated the fact that the likelihood cannot be considered as only criterion when comparing different functional families. This is suggesting methods penalizing complexity (like in AIC and BIC). But we mainly wish to emphazise on the risks took when predicting with trended Kriging: in higher dimensions, we will allways be in an extrapolation situation. Choosing a trend with the help of a small design then seems very risky. This is an argument to consider Ordinary Kriging in the cases where no prior information on the trend is available.

In other respects, we proposed a model combining an additive model and Simple Kriging. The application to a simple industrial test case confirmed that directly kriging the residuals by ML gives a poor result. Our attempt to adapt a method inspired by cross-validation with a single test set gave here a Kriging with different features from ML, apparently accounting well for the non-additive part of the response. However, the question of the robustness to a change of design has not been raised yet. This is a subject to be treated in further works.

Aknowledgements

All the computations have been performed using R [START_REF] Development | R: A Language and Environment for Statistical Computing[END_REF] and the packages geoR [START_REF] Ribeiro | geoR: A package for geostatistical analysis[END_REF], RandomFields [START_REF] Schlather | Simulation and analysis of random fields[END_REF], gam and gstat. This work was conducted within the frame of the DICE (Deep Inside Computer Experiments) Consortium between ARMINES, Renault, EDF, IRSN, ONERA, and Total S.A.

Figure 1 :

 1 Figure 1: Two realizations of the random function -2 ln L(.; Z) corresponding to two simulated response values Z, both with a Gaussian covariance kernel (c g defined hereafter) and covariance parameters ψ = (5, 0.3). Left: ML estimates are close to the actual parameters (bold dot) : ψ ≈ ψ. Right: ML fails to locate the actual parameters: ψ = ψ.

Figure 2 :

 2 Figure2: Left: Comparison between the experimental law (gray boxplots) and the asymptotic law (black lines) for the scale parameter for increasing size of the design. The boxplots for the experimental laws have been done using 1000 simulations, with Gaussian covariance function of parameters ψ = (5, 0.5). For the asymptotic law the median is represented in continuous line and the first and third quartiles in dashed lines. Right: Histogram of the logarithm of relative errors obtained when estimating ψ 1 by ML using 1000 GP realizations. The shape of this histogram suggests that the distribution of relative errors is far closer to a lognormal law than to a Gaussian.

Figure 3 :

 3 Figure3: One realization (dashdot line) of a GP with linear trend is interpolated by Ordinary Kriging, based on 3 designs. The OK mean and 95% confidence intervals are represented by bold lines and dotted lines, respectively. The first design (left) is a regular grid; the associated OK prediction seems satisfying, even if the trend model is misspecified. The second design (center) is formed by 6 points concentrated at the boundaries of the domain; The Kriging predictor fails to capture the shape of the realization at the center of the domain. The third design is made of three points clustered at the center of the domain; OK automatically comes back to the mean value outside of the design and dramatically miss the actual trend.

Figure 4 :

 4 Figure4: One realization (dashdot line) of a GP with linear trend is interpolated by three different Krigings, based on a regular grid (5 points evenly spaced between -1 and 1). The GP is here represented between -2 and 2, such that all three cases can be referred to as extrapolations. UK with affine trend (center) gives accurate results on the whole domain. On the contrary, both OK (left) and UK with quadratic trend (right) give good results between -1 1 but dramatically fail in extrapolation.

Figure 5 :

 5 Figure5: Coplots of z on the Hammersley design X 1 (left) and summary of the additive components and the residuals obtained after application of the backfitting algorithm (right). The additive model is here chosen with a linear function in both directions of x 1 and x 3 , and a smoothing spline in x 2 's direction.

Figure 6 :

 6 Figure6: Coplots representing both X 1 (dots) and X 2 (triangles) designs in projection on all pairs of coordinates. The three graphics illustrate the spacefilling behaviour of X 1 and the D-optimal nature of X 2 . X 2 also appears to be reasonably disconnected from X 1 .

Figure 7 :

 7 Figure7: Standardized residuals at the validation design X 2 (left) and variation of the ISE with respect to the covariance parameter p (right). One can observe that the value of p found by ML at X 1 (p 1 = 0.16) is clearly suboptimal to get an accurate Kriging predictor when extrapolating to X 2 .

Figure 8 :

 8 Figure 8: Left: Histogram of the standardized residuals at the test design F with the model previously obtained by Algorithm 2 (GAM mixed, σ 2 = σ 2 2 , p = p 2 ). Right: Variation of the ISE with respect to the covariance parameter p: the value p 2 = 1.4 previously chosen at X 2 is almost optimal again.

Table 1 :

 1 ML and ISE values on 1000 simulated realizations of GPs with Gaussian covariance function, for relative parameters ψ rel i

Table 2 :

 2 MLE and ISE measures on 1000 simulated GP realizations with Gaussian covariance kernel, for X = X 10 . The aproximation based on Fisher's Information Matrix is still underestimating the estimation variances but is less unprecise than in table[START_REF] Abt | Estimating the prediction mean squared error in gaussian stochastic processes with exponential covariance structure[END_REF].

	ψ 5 0.3	E ψ rel i -0.054 0.012	i	V ar ψ rel i 0.432 0.105 i 0.105 0.085	asymptotical V ar ψ rel i 0.297 0.057 0.057 0.033	i	E[ISE] 0.177
	5 0.4	-0.042 -0.019		0.424 0.058 0.058 0.024		0.340 0.044 0.044 0.014	0.009
	5 0.5	-0.067 -0.013		0.46 0.051 0.051 0.013		0.362 0.036 0.036 0.008	0.0004
	5 0.6	-0.075 -0.007		0.728 0.059 0.059 0.012		0.375 0.032 0.032 0.005	4.e-05
	ψ = (10, 0.4) and #X = 5 we obtain V ar ψ rel i	i =	43.555 3.242 3.242 0.971

Table 3 :

 3 Comparison of minimum -2 log(L) values obtained by fitting three different Kriging models (OK, UK affine, UK quadratic) to three GP realizations. Each realization is drawn from the GP underlying one of the Kriging models. The design is a regular grid on [-1, 1]. The results illustrate that adding degrees of freedom to a Kriging model always lead to a larger value of the maximum likelihood.

			GP	GP +linear t	GP+quadratic t
	kriging type	p	-2 ln(L(p))	p	-2 ln(L(p))	p	-2 ln(L(p))
	OK	0.4082	32.07	0.4445	36.90	0.4595	38.80
	UK, linear t	0.4085	31.89	0.4085	31.89	0.4387	35.80
	UK, quadratic t 0.4084	31.89	0.4084	31.89	0.4084	31.89

Table 4 :

 4 Optimal loglikelihood values and estimated covariance parameters associated with the residuals provided by Algorithm 1 at X 1 with different trend structures. The R 2 values are computed by comparing the residual deviance after fitting the trend only, to the total variance of the response at the design X 1 .

	Model	Loglikelihood Range	σ 2	R 2 adj	R 2	p-value
	1st order Linear + SK	-121.91	1.04 26101.89 0.78 0.82 4.03e-06
	2st order Linear + SK	-100.69	0.048 1381.13 0.97 0.98 6.44e-08
	GAM splines + SK	-76.01	0.048	117.10	-	0.99	-
	GAM mixed +SK	-80.62	0.16	185.71	-	0.99	-

Mind the fact that by negative relative bias we understood an overestimation of ψ.


INTRODUCTION

Quantitative methods of analysis based on statistical models have been used on dam monitoring data for a long time, as reported in particular at the ICOLD Congress. The Hydrostatic-Season-Time method proposed by Electricité de France for analysing pendulums in 1958 (Ferry & Willm) and in 1967 in a fundamental study (Willm & Beaujoint) has proved to be a powerful tool for interpreting the behaviour of concrete dams in particular (ICOLD 1985(ICOLD , 1989(ICOLD , 2000)).

The statistical models involving the use of regression techniques to find correlations between causes and quantified effects account efficiently for mechanical behaviour, but they have several weaknesses :

-theses methods cannot be used to predict variations in the pore pressure rates ; -they do not include the effects of rainfall because they do not take the history of the loading into account. The analysis of ageing behaviour requires more accurate methods of non-ageing behaviour analysis.

To understand the long term behaviour of a dam, it is essential to carefully interpret the dissipative effects which tend to occur due to seepage, in order to be able to distinguish between the effects of factors such as drift, irreversible events and the ageing of the dam and the effects of other factors not involving ageing processes.

THE HYDROSTATIC-SEASON-TIME MODEL

The main components of the basic Hydrostatic-Season-Time (HST) model are as follows:

-the effect H of the reservoir level Z, which is given by a fourth-degree polynomial, H n =a 1 Z n +a 2 Z 2 n +a 3 Z 3 n +a 4 Z 4 n

(1)

-the seasonnal effect S, represented by the sum of sine functions with one-year and six-month period ( =2¡ /365) S n = b 1 sin( t n )+b 2 cos( t n )+ b 3 sin 2 ( t n )+b 4 sin( t n )cos( t n )

-the irreversible effect T (time drift), which is expressed in terms of monotonic time functions. The HST model is extremely robust and always yields satisfactory results. It gives in a simple, immediately usable form an instantaneous non linear function of the reservoir level H and a periodic function S, which is comparable to a delayed response to the annual and half-yearly cycles in which the total external loads occur.

One of the main gaps in this model is the lack of physical information provided by the parameters.

When the model H is linear (a 2 =a 3 =a 4 =0), the coefficient a 1 is equal to the coefficient ¢ :

H n =¢ Z n (3)
If the value of a is close to zero, this means that the reservoir level will have no instantaneous effect Delayed response analysis of dam monitoring data S. Bonelli

LMA-CNRS, Marseille, France

P. Royet

Cemagref, Aix-en-Provence, France ABSTRACT: Long series of monitoring data are obtained during the routine operation of a dam. The Statistical Hydrostatic-Season-Time models which have been widely used for analysing and interpreting these large sets of data are not suitable for hydraulic measurement analysis, as they do not account for the delayed responses of piezometers or cells. This paper describes a model for performing delay analysis on pore pressure measurements. The method is based on Darcy's law and Richard's equation of seepage and involves the use of a linear dynamic system accounting for the contribution of non-ageing factors, i.e., the reservoir level and rainfall events, to the pore-pressure variations.

hal-00270192, version 1 -4 Apr 2008

on the variations. If a is close to one, this means that the instantaneous variations in the dam monitoring data will be of a similar amplitude to those of the reservoir level.

The three parameters (a 2 ,a 3 ,a 4 ) involved in the explanatory variable H and the four parameters involved in the explanatory variable S are never interpreted in the usual hydraulic data analysis.

The polynomial expression for the effects of the water level (1) was originally based on a mechanical analysis of the effects of the water level on the displacement of an arch dam, based on the resistance of the materials.

This explanatory variable is often used by default in hydraulic data analysis, but a fourth-degree polynomial relationship between a piezometric head and the water level is not mechanically justified.

The only advantage of the previous method is that it provides a means of explaining the readings taken when the explanatory variable is not Z-Z 0 but (Z-Z 0 )E(Z-Z 0 ), where E(Z-Z 0 )=1 if Z£ Z 0 and 0 otherwise, and Z 0 is the threshold piezometric level at which the effect occurs, which depends on the permeability and the slope of the ground and the position of the instrument bottom of the stand pipe piezometer or quotation of pose of the cell.

Seasonal factors are known to affect arch dams : the temperature variations occurring between cold and hot seasons are closely correlated with the downstream drift of the dam, for instance. As it is extremely difficult to use and even to obtain accurate temperature measurements, a periodic seasonal law based on average temperature variations is used.

The hydraulic measurements made at earth dams used for water supply or irrigation purposes also show the effects of the latter factor. The variable S is in fact simply the sum of the two first terms of a Fourier series development (to within the nearest constant):

S n = A 1 sin(¤ (t n +d 1 ))+A 2 sin(2¤ (t n +d 2 )) (4)
The rainfall should be taken into account when interpreting the readings obtained with piezometers, which are instruments used on dams of all kinds, especially for foundation monitoring. There exists a simpler method consisting of taking the rainfall during the last ten days [START_REF] Crépon | An analytical approach to monitoring[END_REF], but this is a purely statistical approach.

Figures 1 and2 give an example of piezometric data obtained on the downstream toe of a homogeneous earthdam.

The variations of piezometric level occurring during the first few years seem at first sight to be in parallel with the water level (fig. 1). This would give a straight line or at least a cloud of aligned points on the graph giving the variations in the piezometric level vs the water level, but it is not in fact the case (fig. 2). Even a high order polynomial would obviously not account for this result.

This example, which we will examine in greater detail below, shows what delayed responses consist of and why neither model (1) nor model (3) can account satisfactorily for these responses.

Delayed effects are due to dissipative behavior (viscoelasticity, seepage ,etc.), and are therefore irreversible. By including the impulse response of a non ageing visco-elastic materials in the statistical analysis, it is possible to analyze the creep deformations occurring in a concrete dam (Dobosz, in ICOLD 1985).

The impulse response of a semi-infinite porous medium given by the Boussinesq equation can also be used to analyze the flow through an earth dam [START_REF] Brunet | Analyse retard des mesures de pression interstitielle et estimation des perméabilités d'un noyau de barrages en terre[END_REF][START_REF] Fabre | Méthodes d'analyse et de surveillance des déplacements observés par le moyen de pendules dans les barrages, R.118, Q.21[END_REF]. At a more general level, one can take the derivative of the loading as the explanatory variable [START_REF] Crépon | An analytical approach to monitoring[END_REF].

It can be seen from Figure 2 that a cycle in which the water level rises and falls is a dissipative one (it involves hysteresis) : the path taken as the water level rises (phase 3) and falls (phase 5) is not the same. For this reason, some measurements can indicate that the pore pressure has increased while the water level was decreasing, and vice-versa.

This well-known paradox is due to the presence of air trapped inside the body of the dam. This situation has been observed in situ, and has also been found to occur under laboratory conditions (Windish & al, 2000).

DESCRIPTION BY A STATIONARY LINEAR DYNAMIC SYSTEM

Seepage processes occur in response to continuous reservoir level variations and rainfall events. Richard's equation shows that the pore pressure does not depend on the instantaneous value of the loading, but on the convolution integral of an impulse response (which still remains to be identified) and the loading conditions (reservoir level, rainfall).

In attempting to describe the variations in the pore pressure measured in situ, this expression led us to search for a more external description in terms of a stationary linear dynamic system, for the following reasons :

-linear: the class of structure in question is that of dams during a routine period of operation undergoing the specific loading levels for which they were designed; external constraints are therefore assumed to result in reversible deformations which are highly unlikely to be large enough to affect the stability or the resistance of the dam; this justifies the use of a linear approach in the first stage.

However, a linear approach is obviously more justifiable under saturated or quasi-saturated than under saturated conditions ; -dynamic: the two main constraints affecting the flow in a hydraulic dam are the variations in the water level and the rainfall; in view of the dissipative nature of the seepage processes, the history of these events has to be taken into account in order to explain the water levels existing at a given instant, and not just the values of these two parameters at the same instant. What is required here is a dynamic description in terms of differential equations ; -stationary: to be able to quantify the drift with time, i..e , the changes occurring under constant conditions (including the ageing), it is necessary to first quantify the stationary changes, i.e., those resulting from external factors regardless of the time; an invariant stationary system, i.e., one which is independent of the time origin, was therefore adopted.

Here the time plays only the kinetic role characteristic of dissipative phenomena (in the sense that the successive events occur at a certain speed), but it has no geological significance (in the sense that the time origin, and hence the age of the system, is not taken into account); here one can refer to the theory of rheology for an exact definition of aging and non ageing behavior. This approach has already yielded some results. The response obtained is the sum of a transient term giving the initial conditions, the response to the variations in the water level and the response to the rainfall events.

The term "ageing" seems to be more suitable than the term "irreversible", which is often used to characterize the effects of time (the age of the dam) on hydraulic measurements. Seepage is by nature an irreversible phenomenon.

Accommodation is one of the characteristics of a stationary linear dynamic system, and it renders the use of the explanatory variable S pointless.

In a permanently operating regime, the response to a periodic signal (harmonic test) will also be periodic (definition of accommodation), and the period will be identical, whereas the amplitude will be different and the response will be delayed.

If the input is harmonic

Z(¥ )=sin(¦ §¥ ) (5) the response will be h Z (¨,¥ ,© )*Z(¥ )= gsin(¦ (¥ +d)) (6)
where is the static structural decrease in damping of the amplitude, g is the dynamic decrease in the amplitude, and d is the delay. The accommo-dation is precisely what the seasonal variable S was designed to model (4).

One might naturally expect an explanatory variable expressed in dynamic terms to be an improvement over methods of interpreting hydraulic measurements purely in terms of seasonal effects.

AN ORDER ONE DELAY MODEL

The delayed effect H of the water level is proportional to the convolution of the impulse response of the dam structure and the water level Z.

The simplest stationary linear dynamic system is an order one system.

Approximating the impulse response h Z by taking a system defined by a characteristic time T Zx leads to

H(Z,x,t)= x Z * (x,t,Z) (7) Z * (x,t,Z)= 1 T Zx 0 t exp - t-t' T Zx Z(t')dt' (8)
The coefficient 0 ! x 1 is the static structural damping of the amplitude. It reflects the efficiency of the drainage system or grout curtain, as well as the position of the measuring point in relation to the reservoir level.

A coefficient x approximately equal to unity will mean either that the instrument has been placed close to the reservoir level or that the actual drainage outlet is far away, which may constitute a weakness of the drainage system. The characteristic time T Zx integrates several items of information about the zone situated between the reservoir surface and the drainage point : the efficiency of the drainage system or the grout curtain, and the diffusive properties of the materials, which depend on the permeability and the state of saturation of the ground and the compressibility of the water.

A very large characteristic response time indicates either that the ground is not saturated (S<85%), the degree of permeability is very low or the drainage distance is very long.

When the characteristic response time is very short, the process is taken to be an instantaneous one corresponding to (3).

The delayed effect P of the rainfall is proportional to the convolution of the impulse response of the dam structure and the rainfall Q.

Approximating the impulse response in terms of an order one dynamic system defined by a characteristic time T Qx gives

P(Q,x,t)=T kx Q * (Q,x,t) (9) Q * (Q,x,t)= 1 T Qx" # $ 0 t exp% & ' )( 0 1 - t-t' T Qx Q(t')dt' (10)
The coefficient T kx is the drainage time in the part of the dam in which piezometers have been installed. The characteristic time T Qx integrates several kinds of information about the zone situated between the surface of the ground and the actual drainage outlet point, such as the efficiency of the drainage system and the diffusive properties of the materials.

The following are suitable expressions for the discrete convolutions involved here :

Z * n+1 =Z * n +2 Z n+1 + (Z n -H n -2 Z n+1 T Zx 2 t n )(1-e -3 t n /T Zx ) (11) Q * n+1 =Q * n +(Q n+1 -Q * n )(1-e -3 t n /T Qx ) ( 12
)
where 2 Z n+1 = Z n+1 -Z n and 2 t n+1 = t n+1 -t n .

VALIDATION

An example confirming the validity of the model ( 11) as compared to an exact solution [START_REF] Carslaw | Conduction of heat in solids[END_REF] of the diffusion problem is given in figures 3, 4, 5 and 6. Figures 3 and4 show the response to a stepwise increase in the reservoir level. Figure 3 gives the spatial profile of the response to one stepwise increase in the water level at several successive instants. Figure 4 gives the pattern of response with time. These two figures show the static damping of the amplitude.

Figures 5 and6 give the responses to a harmonic variation of loading imposed by the reservoir level. They show: 1) the delayed effect and the dynamic and static decrease in the amplitude, 2) the fact that the delayed response model (11) accounts efficiently for periodic loading.

Similar results have been obtained with this delayed response model ( 12) in a flow loading situation [START_REF] Bonelli | Solutions exactes et approximées de problèmes de diffusion, application à l'analyse différée des mesures hydrauliques de barrage[END_REF].

By studying the exact solution, it is possible to specify T Zx and 4 x , given what is known about the characteristics of the dam :

4 x =1-x L , T Zx =T Z 1-4 2 x 6 , T Z = L 2 D ( 11
)
where x is the distance between the measuring instrument and the upstream face, L is the distance between the drainage system and the upstream face (the length of the drainage path), T Zx is the characteristic response time of the measuring instrument, T Z is the characteristic response time in the region to be instrumented and D is the diffusion coefficient characteristic of the region to be instrumented (fig. 7).

It is now possible to assess the consistency of the responses recorded by the instruments placed within the same zone. The permeability can also be assessed by quantifying the capacity coefficient, based on the amount of trapped air present.

APPLICATION TO DETERMINING THE EFFECTS OF THE WATER LEVEL ON PORE-PRESSURE CELLS RESPONSE

The analysis of the readings obtained with a pore pressure cell placed in the body of an earth dam (location CV13 in fig. 10) makes it possible to account accurately for the behavior of the dam (figs. 8 and 9).

Figure 8 shows the delay in the responses. Figure 9 shows the hysteretic characteristic of the dissipative, and therefore delayed, behavior observed. The static HST model cannot account for behavior of this kind. Measurements obtained with the six cells placed inside the structure of the Alzitone earthdam (France) (fig. 10) were analyzed and interpreted (table 1).

The permeability of the region within which the water level fluctuates (10 -8 m/s) was found to be one order of magnitude lower than that of the permanently saturated zones (10 -7 m/s).

The permeabilities deduced from laboratory tests, in situ Lefranc tests and finite element simulations range between 10 -6 and 10 -7 m/ s. Measurements obtained with six cells set within the structure of the Chamboux earthdam (France) (fig. 11) were analyzed and interpreted (table 2). The effects of the water level were instantaneous on the three cells located at the interface with the foundations, giving a permeability value greater than 10 -8 m/s, whereas the three cells located inside the structure of the dam detected a delayed effect. Since these cells were placed within the water level fluctuation zone, the permeability in the permanently saturated region can be taken to have been one order of magnitude less, namely 10 -7 m/s.

The projected permeability (laboratory tests) yielded a vertical permeability of 10 -8 m/s.

In both cases studied here, the accuracy of the results obtained with the delayed response model was all the more impressive as this was not actually the purpose for which the model was designed. 16).

During the filling of the dam, 120 measurements were carried out during a period of 150 days. The water level chart included 212 measurements recorded during a period of 329 days, including 179 days during which measurements were possible before the start-up. After adjustment of the model ( 9), we obtained T Zx =83 days, 7 x =0.44 and T Z =613 days (figs. 12 and 13).

During the operating phase, 570 measurements were carried out during a period of 2280 days. This gave T Zx =159 days, 7 x =0.22 and T Z =1000 days (figs. 14 and 15).

The simulation obtained with the dynamic model is remarkably accurate, given the extreme simplicity of the model. The results were compared with those obtained with the static HST model (3) (which was found to be ineffective) and with the measured data.

The results of the analysis of the readings obtained with the piezometers located in the body of the Alzitone dam (fig. 16) are summarized in table 3. The results of the cell reading analyses are also included in this table, and the great consistency of these results is worth noting. The data obtained with four piezometers located in abutments of La Verne earthdam (France) were analyzed in order to assess the rainfall effect model ( 12). Two of the instruments were placed on the right bank (PZ17 and PZ14) and two on the left bank (PZ18 and PA4).

The effect of the water level was found to be instantaneous (T z8 0). This finding still remains to be interpreted in the light of the structural and operational data available about this dam.

The model used accounted only satisfactorily for the variations observed (figs. 17 and 19). Since the drainage path was not known, it was impossible to interpret the parameters reflecting the effects of the rainfall.

When only the effects of the water level was considered, peaks were observed in the piezometric level corresponding to rainfall events (fig. 19). The graph on which corrected water level values were plotted shows that the rainfall model needs to be improved, since it showed only some of the peaks and hollows which occurred (fig. 20).

Although the delayed model ( 12) for the effects of the rainfall is still an improvement over the previously available models, the results are not as outstanding as those obtained in the case of the water level.

There are two possible reasons for this difference in the efficiency of the model between the two cases tested here: 1) the rain-water is not entirely infiltrated: the quantities of infiltrated and streaming rain-water depend the amount of rainfall, and on the slope and the permeability of the ground; 2) the path taken by the infiltrated water begins at the surface, and crosses an unsaturated zone, which it is difficult to account for using a linear model with a constant diffusivity coefficient. This table gives some orders of magnitude. The rainfall certainly affected the readings obtained on some of the instruments. The response times to the rainfall were distinctly longer than the response times to variations in the water level. 11 CONCLUSION Some important conclusions can be drawn from the results presented above. The delayed response model accounted satisfactorily for seasonal effects on all the piezometric and pore pressure data analyzed. The delayed response model always yielded greater static damping of the amplitude @ x than the static HST model. This is consistent with the previous comment, since the present model takes the variations in the water level into account, whereas the seasonal variations in the piezometric and pore pressure data are not attributed to the variations in the water level in the HST analysis.

Underestimating @ x is not conducive to safety, and the use of the delay model is therefore a must from this point of view.

All instrument readings, whether they are obtained using cells set inside the structure of the dam or piezometers sounding the outside of the dam, the foundation, the abutments and the banks, are liable to exhibit delayed response effects.

The results obtained with the delayed response model were found to constitute an improvement over the previously available methods as far as the effects of the rainfall were concerned, but the results obtained on the effects of the water level were even more satisfactory.

The delayed response model can therefore be used to perform mechanical analyses of the data obtained using instruments set inside the structure of dams where the drainage path can be determined and where the rainfall has negligible effects.

Figure 1 .

 1 Figure 1. Example of piezometer data.

Figure 2 .

 2 Figure 2. Piezometric head v.s. water level. The path taken as the water level rises and is not the same.

Figure 3 .

 3 Figure 3. Validation of the delay model as compared to an exact solution, spatial profile of the response at several instants.

Figure 4 .

 4 Figure 4. Validation of the delay model as compared to an exact solution, pattern of response with time.

Figure 5 .

 5 Figure 5. Validation of the delay model as compared to an exact solution, pattern of response to a harmonic variation of loading imposed by the reservoir level.

Figure 6 .

 6 Figure 6. Validation of the delay model as compared to an exact solution, harmonic test, response v.s. solicitation.

Figure 7 .

 7 Figure 7. The parameters T Zx and A x can be specified for an instrument placed in the body of the dam.

Figure 8 .

 8 Figure 8. Example of data for a cell located in the body of the dam, cell level v.s. time. The delay is about 23 days.

Figure 9 .

 9 Figure9. Example of data for a cell located in the body of the dam, cell level v.s. water level. The path taken as the water level rises (phase 2) and falls (phase 1) is not the same.

Figure 10 .

 10 Figure 10. Location map of the cells, longitudinal profil from downstream (Alzitone dam).

Figure 11 .

 11 Figure 11. Location map of the cells, cross section (Chamboux dam)

Figure 13 .

 13 Figure13. Delay analysis of P1 piezometer during the impounding phase, piezometric head v.s. water level. The path taken as the water level rises and falls is not the same.

Figure 14 .

 14 Figure 14. Delay analysis of P1 piezometer during the exploitation phase, piezometric head v.s. time.

Figure 15 .

 15 Figure15. Delay analysis of P1 piezometer during the exploitation phase, piezometric head v.s.water level. The path taken as the water level rises and falls is not the same.

Figure 16 .

 16 Figure 16. Location map of piezometers in the plan of the dam site (Alzitone dam).

Figure 17 .

 17 Figure 17. Delay analysis of PZ17 piezometer during exploitation, piezometric head v.s. water level (La Verne dam).

Figure 18 .

 18 Figure18. Delay analysis of PZ17 piezometer during exploitation without rainfall effect, piezometric head v.s. water level. Peaks were observed in the piezometric level corresponding to rainfall events.

Figure 19 .

 19 Figure 19. Delay analysis of PZ17 piezometer during exploitation, rainfall effect v.s. time

Table 1 .

 1 Effects of the water level, results of delayed response analysis of data obtained with cells (Alzitone dam) Cell

		5 x T Zx (days)	T Z (days)	L (m)	D (10 -5	k (10 -8 m/s)
						m 2 /s)
	Elevation 52, unsaturated zone (water level fluctuation)
	CV33	0.27	33	211	27	4	1.5
	CV23	0.38	33	228	32	5	2.3
	CV13	0.39	23	160	33	8	2.7
	Elevation 40, saturated zone
	CV32	0.47	25	191 103	60	14
	CV22	0.49	33	260 106	50	13
	CV12	0.38	10	73	88	100	21
	Foundation
	CV31	0.43	28	210 141	100	24
	CV21	0.43	56	409 140	60	11
	CV11	0.51	32	259 163	100	14

Table 2 .

 2 Effects of the water level, results of delayed response analysis of data obtained with cells (Chamboux dam) Cell

		6 x T x (days)	T (days)	L (m)	D (10 -5	k (10 -8 m/s)
						m 2 /s)
	Elevation 488.50, dam/foundation interface
	C1	0.44	0	-	80	-	-
	C2	0.28	0	-	78	-	-
	C3	0.07	0	-	75	-	-
	Elevation 497, body of the dam
	C4	0.25	35	225	17	1.5	0.9
	C5	0.21	35	215	19	2.0	1.3
	C6	0.06	43	256	20	1.8	1.2
	8 APPLICATION TO DETERMINING THE
	EFFECTS OF THE WATER LEVEL ON
	PIEZOMETERS RESPONSE
	The piezometer mentioned above was analyzed in
	two phases (piezometer P1 placed on the down-
	stream toe of the Alzitone earth dam, evolution figs.
	1 and 2, situation fig.

Table 3 .

 3 Effects of the water level, results of delayed response analysis of the piezometer data and some of the data obtained with cells (Alzitone dam)

	Instru-ments	6 x T Zx (days)	T Z (days)	L (m)	D (10 -5	k (10 -8 m/s)
						m 2 /s)
	Left side
	CV33	0.27	33	211	27	4	1.5
	CV32	0.47	25	191 103	60	14
	PID9	0.42	45	325 111	40	14
	PID10 0.23	45	285 100	40	17
	Middle
	CV23	0.38	33	228	32	5	2.3
	CV22	0.49	33	260 106	50	13
	PID8	0.29	30	200	76	30	63
	PID3	0.21	31	193	90	50	2.4
	Right side
	CV13	0.39	23	160	33	8	2.7
	CV12	0.38	10	73	88	100	21
	PID7	0.28	18	119	75	50	14
	PID6	0.53	22	183 150	100	37
	PID5	0.16	91	558 107	20	8.9
	9 APPLICATION TO DETERMINING THE
	EFFECTS OF RAINFALL ON PIEZOMETERS
	RESPONSE

Table 4 .

 4 Effects of the rainfall: results of delayed response analysis of piezometer data. The effects of the water level were instantaneous (La Verne dam)

	Pie-zome-	Water level	Rainfall effect	9	T Zx (days)	T Kx (days)	T Qx (days)
	ter	effect
	PZ 17	60%	24%	0.32	0	273	38
	PZ 14	54%	23%	0.14	0	147	53
	PZ 18	62%	19%	1.00	0	815	39
	PA 4	21%	52%	0.46	0	1016	32
	10 APPLICATION TO DETERMINING THE
	EFFECTS OF THE WATER LEVEL AND
	THE RAINFALL ON PIEZOMETERS
	RESPONSE
	The results of the analysis of the piezometer data
	obtained beneath the Chamboux dam (fig. 20) are
	summarized in table 5. Since the drainage path was
	not known, it was impossible to interpret the pa-
	rameters.

Table 5 .

 5 Effects of the water level and rainfall: results of delayed response analysis of piezometer data (Chamboux dam)

	Pie-zome-	Water level	Rainfall effect	x 9	T Zx (days)	T Kx (days)	T Qx (days)
	ters	effect
	PD1	87%	6%	0.77	5	77	39
	PD2	53%	6%	0.33	4	33	19
	PD3	21%	6%	0.21	3	24	13
	PD4	17%	45%	0.04	0	54	52
	PD5	12%	54%	0.05	3	76	34
	PD6	24%	46%	0.09	0	106	54
	PD7	86%	10%	0.85	7	156	66
	PD8	21%	49%	0.20	30	198	65
	PG1	88%	6%	1.00	3	109	63
	PG2	73%	9%	0.30	5	34	16
	PG3	23%	36%	0.07	3	39	22
	PG4	16%	49%	0.05	11	54	40
	PG5	32%	31%	0.71	3	608	83