Ekaterina Sergienko

Paul Lemaître

Aurélie Arnaud

Daniel Busby

Fabrice Gamboa

Reliability sensitivity analysis based on probability distribution perturbation with application to CO2 storage

Keywords: sensitivity analysis, reliability analysis, uncertainty analysis, failure probability

Introduction

Carbon Capture and Storage (CCS) stands for the collection of CO 2 from industrial sources and its injection into deep geological formations for a permanent storage. There are three possible sites for injection: unmined coalbed formations, saline aquifers and depleted oil and gas reservoirs [START_REF] Polson | Process for tracking the evolving perception of risk during CO 2 storage projects[END_REF]. Nevertheless, the following principal environmental question arises: what is the probability that CO 2 will remain underground for hundreds to thousands of years after its capture and injection into a storage formation?

The primary risk of CO 2 geological storage is unintended gas leakage from the storage reservoir [START_REF] Le Guen | CO 2 storage: managing the risk associated with well leakage over long timescales[END_REF][START_REF] Bowden | Assessing reservoir performance risk in CO 2 storage projects[END_REF]. In this work, we focus on a leakage from the storage formation through a fault or a fracture. This can happen when the reservoir pressure is higher than the caprock fracture pressure. Numerical modelling and simulation has become an integral component of CO 2 storage assessment and monitoring. The reservoir simulation models are constructed based on the reservoir and production data. They are used to predict and analyze CO 2 plume distribution and the reservoir pressure development during the injection and the storage periods. With the help of the numerical simulation models it is possible to forecast CO 2 storage performance and to evaluate the risks of a possible leakage.

Risk and uncertainty analysis has been recognized as a principal part of safety and risk assessment. Generally speaking, when the main sources of uncertainty have been identified, Uncertainty Analysis (UA) is focused on quantifying the uncertainty in the model output resulting from uncertainty in the model inputs. At the same time, Sensitivity Analysis (SA) aims to identify the contributions of each model uncertain input to the variability of the model output. Uncertainty analysis may be equally performed to assess the reliability of the system. A typical example of a failure probability estimation in the CO 2 storage risk analysis is the estimation of the probability of exceeding the caprock fracturing pressure during the CO 2 injection phase. If we denote P reservoir the reservoir pressure and P f racture the caprock fracturing pressure, then we consider the following failure probability: p f = P(P reservoir ≥ P f racture ).

In practice, uncertainty and sensitivity analysis require a large number of reservoir simulator runs to explore all the input variables space. However, higher accuracy of a simulator usually results in a higher simulation time. One simulator run can take from few minutes up to several hours or even days. Therefore, when the simulation time becomes too high, uncertainty analysis may become unfeasible. For this reason, in this work, in order to estimate the failure probability p f for the expensive reservoir simulator, we use Gaussian Process (GP) response surface model, also known as kriging [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Welch | Screening, predicting, and computer experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF]. The kriging method was originally introduced in the field of geostatistics by Krige in the 1950's [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF] and formalized in 1960's by Matheron [START_REF] Matheron | Principles of geostatistics[END_REF]. In [START_REF] Sacks | Design and analysis of computer experiments[END_REF] Sacks et al. proposed the statistical approach to uncertainty analysis of complex computer codes referred to as the Design of Computer Experiments. In a nutshell, the approach consists in building an approximation of the reservoir simulator input/output relationship starting from a set of simulation runs at a carefully chosen input variables configurations referred to as experimental design or training set. The obtained response surface model can then be used to predict the model output for a new non simulated input with a negligible computational time. Therefore, uncertainty and sensitivity analysis become affordable.

One of the most challenging problems in risk analysis is to identify the failure region and to compute the failure probability. However, the failure probability will usually depend strongly on the probability distribution of the input variables. Reliability Sensitivity Analysis can help in understanding the relationship between each input variable uncertainty and the failure probability. The problem is to identify the set of input variables that need to be well characterized in terms of uncertainty distribution.

For the time being, numerous approximation and simulation methods are available for estimating the failure probability such as First / Second Order Reliability Methods, Monte Carlo sampling, importance sampling, directional sampling, subset simulation, etc. [START_REF] Lemaire | Structural reliability[END_REF][START_REF] Melchers | Importance sampling in structural systems[END_REF][START_REF] Ditlevsen | Directional simulation in gaussian processes[END_REF][START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. However, there are few sensitivity analysis methods developed for the failure probability analysis. The widely used methods for sensitivity analysis are based on a variance decomposition of the output. Given the probability distribution of the input variables, Sobol indices are expressed by the ratio of the variance due to a given input on the total output variance [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Knowing the probability distribution of the input variables and the output, we can so define sensitivity indices for each of the input variables. Variance based indices are usually interesting for measuring the input/output sensitivity, however they can be poorly relevant to our problem of evaluating the impact on failure probability.

There have been few attempts to develop sensitivity analysis methods well suited for reliability analysis. First, as complementary results of the First Order Reliability Method, sensitivity to the distribution of the input variables can be obtained. Sensitivity is expressed as the partial derivative of the reliability index β [START_REF] Lemaire | Structural reliability[END_REF]. Another approach was proposed by Morio [START_REF] Morio | Influence of input pdf parameters of a model on a failure probability estimation[END_REF]. Therein, the author uses the variance decomposition and Sobol' sensitivity indices to study the rate of change in the failure probability due to the changes of the input distribution density parameters. Borgonovo et al. [START_REF] Borgonovo | Moment independent importance measures: New results and analytical test cases[END_REF] suggested some moment independent importance measure in the reliability analysis. This measure does not involve the variance. For a fixed variable x i it quantifies the effect of knowing x i by computing the L 1 norm between the unconditional joint density f x(•) and the conditional density

f x|x i (•).
In this paper, we study a moment independent approach for sensitivity analysis of a failure probability [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]. The influence of the input variables on the failure probability is obtained by perturbing the prior probability density function f x( x). In particular, we estimate the effect of the perturbation on the value of the failure probability p f . Here, we propose to distinguish distributions classes by their supports. For the case of a bounded support, such as the uniform or the triangular distributions, the main source of uncertainty is about the boundaries of the support. On the other hand, in the case of infinite support, such as normal or log-normal distributions, the main source of uncertainty comes from the distribution parameters, such as mean and variance. The estimation method has the advantage of being very efficient in terms of number of simulator calls. In order to estimate the sensitivity indices for all the input variables, the performance function is evaluated only once on a Monte Carlo sample used to estimate the reference failure probability p f . Our paper is organized as follows. First, we introduce some density perturbations for different families of distributions. Later, we introduce the technique to compute a perturbed failure probability using the same Monte Carlo sample. This is based on an inverse importance sampling technique [START_REF] Hesterberg | Estimates and confidence intervals for importance sampling sensitivity analysis[END_REF]. Finally, we present the formulation for the moment independent sensitivity indices and demonstrate its applicability on an analytical and a CO 2 storage reservoir case examples.

Density perturbation influence to failure probability

Let us denote by g( x) the performance function of the system, x ∈ Ω ⊂ R d is a set of independent input variables with the joint density f x( x) = d i=1 f x i (x i ). The failure probability is expressed as:

p f = P(g( x) ≤ 0) = E f x I g( x)≤0 = Ω f f x( x)d x,
where

Ω f = { x ∈ Ω : g( x) ≤ 0} is the failure region.
In general, the distribution density f x(•) is provided by experts on the basis of some indirect measurements or some limited observation data [START_REF] Lemaître | Hiérarchisation des sources d'incertitudes vis à vis d'une probabilité de dépassement de seuil[END_REF]. Here, we study how a perturbation of the original probability density f x(•) affects the failure probability of the system p f . We assume that the input variables x i , i = 1, . . . , d are independent random variables with marginal densities

f x i . So that, f x( x) = d i=1 f x i (x i
). This work is inspired by previous work [START_REF] Lemaître | Hiérarchisation des sources d'incertitudes vis à vis d'une probabilité de dépassement de seuil[END_REF][START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]. Originally the method was applicable only to the normal distributions. Here, the method is extended to the more general case of an exponential family. The objective of this study is to estimate the influence of the input random variable from the vector x to the failure probability of the system p f . The proposed approach consists in perturbing the original density for a given fixed variable x i while keeping constant the probability density functions for all the other variables x-i = (x 1 , . . . , x i-1 , x i+1 , . . . , x d ). Then, a new value for the failure probability is computed. If this new value p iδ differs significantly from the reference value p f , it means that this selected input variable x i is influential. Conversely, if the new failure probability p iδ is close to p f , then the input x i has low influence on the failure probability.

Density perturbation for an exponential family

To define the probability density perturbation, first recall the definition of an exponential family. An exponential family is a set of distribution having density function that can be expressed in the form of:

f x (x|θ) = h(x) exp (η(θ)T(x) -A(θ)) , (1)
where θ is a vector of distribution parameters, η(θ) is a natural parameter, T(x) is a vector of sufficient statistics, h(x) is an underlying weight function and A(θ) is the cumulant generating function [START_REF] Barndorff-Nielsen | Information and exponential families in statistical theory[END_REF]. The cumulant generating function ensures that the distribution integrates to one, i.e.:

A(θ) = log h(x) exp(η(θ)T(x))µ(dx),
where µ is the reference measure (for example Lebesgue measure). It could be continuous or discrete. In this paper, we mainly consider continuous measure. However, all the calculations are valid for a discrete measure as well. Exponential family contains most of the standard discrete and continuous distributions that we use for practical modelling, such as the normal, Poisson, Binomial, exponential, Gamma, multivariate normal, etc.

To define the density perturbation for this family of distributions, we use some ideas coming from information theory [START_REF] Cover | Elements of information theory[END_REF]. Kullback-Leibler (KL) divergence is used to measure the magnitude of a perturbation. KL divergence quantifies the "closeness" of two probability distribution P and Q. Suppose that P and Q are continuous probability distributions with densities p( x) and q( x) (with respect to Lebesgue measure). Then, the KL divergence between P and Q is given by:

D KL (P, Q) = ∞ -∞ p( x) log p( x) q( x) d x (2)
For δ > 0, let us denote for the variable x i the new (perturbed) density as f x iτ (•),. We select f x iτ (•) in such way that:

D KL ( f x iτ , f x i ) = δ. (3)
Possible values of the perturbation δ may be restricted by some inequalities on Kullback-Leibler divergence [START_REF] Dragomir | Some inequalities for the Kullback-Leibler and χ 2 -distances in information theory and applications[END_REF]. If we define the function r( x) := p( x) q( x) , x ∈ Ω ⊂ R d and assume that 0 < r < r( x) < R for all x ∈ Ω ⊂ R d . Then, according to [START_REF] Dragomir | Some inequalities for the Kullback-Leibler and χ 2 -distances in information theory and applications[END_REF], we have:

D KL (p, q) ≥ 0 D KL (p, q) ≤ (R -r) 2 4rR = δ max
According to these inequalities, we choose δ ∈ [0, δ max ], where δ max = (R-r) 2 4rR can be computed precisely. Let us consider the original density f x i (x) = h(x) exp(η(θ)T(x) -A(θ)) from exponential family. In order to stay in the same family, we propose to restrict the choice of possible perturbations among the following class of densities:

f x iτ (x) = exp(τT(x) -ψ(τ)) f x i (x) = h(x) exp (T(x)(η(θ) -τ) -(A(θ) + ψ(τ))) , (4)
Here, τ is a constant depending on δ (it is chosen under the condition ( 3)). The function ψ(τ) is a normalization function and it may be expressed as:

ψ(τ) = log ∞ -∞ exp(τT(t)) f x i (t)dt .
It is the cumulant generating function for the perturbed probability distribution f x iτ (x). Moreover, if µ and σ 2 are the mean and the variance of the original probability distribution, then:

ψ(0) = 0 ψ (0) = µ ψ (0) = σ 2 (5)
We aim to perturb f x i in such a way that the KL divergence between the original density f xi and the perturbed density

f x iτ is equal to δ. Notice that: D KL ( f x iτ , f x i ) = ∞ -∞ f x iτ (t) log f x iτ (t) f x i (t) dt = ∞ -∞ f x iτ (t) (τT(t) -ψ(τ)) dt (6)
= τψ (τ) -ψ(τ).

Hence, τ should satisfy the equation:

τψ (τ) -ψ(τ) = δ. (7)
Let τ * = τ(δ) be one solution of [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF]. We use this parameter in order to define the perturbed density modification f x iτ (•) defined by (4). Now, let us consider the function:

G(τ) = τψ (τ) -ψ(τ) -δ.
This function has a global minimum at τ = 0:

G (τ)| τ=0 = 0, G (τ)| τ=0 = ψ (τ) > 0 and G(0) = -δ < 0 for δ > 0. Moreover, G (τ) = τψ (τ): G (τ) < 0, (τ < 0) and G (τ) > 0, (τ > 0)
. Thus, G(τ) is strictly decreasing for τ < 0 and G(τ) is strictly increasing for τ > 0. Hence, the function G(τ) has not more than two zeros τ 1 < 0 and τ 2 > 0, if both of them τ 1 and τ 2 fall into domain of the function ψ(τ).

For every fixed level of δ, we can study two possible effects of the perturbation (4). We denote the corresponding perturbed densities by f x iτ 1 and f x iτ 2 . Then, the joint perturbed probability density is expressed as:

f xiτ j ( x) = f x iτ j d k=1,k i f x k (x k ), j = 1, 2.
The corresponding value of the perturbed failure probability p iδ j ( j = 1, 2) can be computed as the following integral:

p iδ j = E f xiτ j I g( x)<0 = I g( x)<0 f xiτ j d x, j = 1, 2. (8)
In the same way, the interaction effect can be estimated by perturbing two variables x i and x j at the same time by δ 1 and δ 2 respectively. Suppose that τ(δ 1 ) = (τ 1δ 1 , τ 2δ 1 ) and τ(δ 2 ) = (τ 1δ 2 , τ 2δ 2 ) are the solutions of equation [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF], where δ 1 and δ 2 are the perturbations of KL divergence (3) for the variables x i and x j , respectively. The new joint probability density function is:

f xi j,τ(δ 1 ),τ(δ 2 ) ( x) = f x iτ(δ 1 ) f x jτ(δ 2 ) d k=1,k i, j f x k (x k ). ( 9
)
The corresponding value of the perturbed failure probability is estimated in the same way by putting in [START_REF] Matheron | Principles of geostatistics[END_REF] the new joint probability density f xi j,τ(δ 1 ),τ(δ 2 ) (•).

In the next section, we introduce a method to estimate efficiently the perturbed failure probability p iδ using the same Monte Carlo sample as for the estimation of the original failure probability p f . First, we study the effect of the perturbation for different probability distributions.

Resulting distributions

Here, we provide a summary table for the considered distributions from exponential family with the resulting perturbed distribution. For the case when the natural parameter η(θ) is a vector of functions (like normal and lognormal distributions), we propose to analyze the perturbation effect separately for each of the components of η(θ). For example, for the normal distribution:

T(x) = x x 2 , η(µ, σ) = µ/σ 2 -1/2σ 2 .
Then we propose to analyze two different density perturbations (4) by taking:

τ = τ 0 or τ = 0 τ .
Table 1 provides the results obtained for some distributions from the exponential family: normal, log-normal, exponential and Poisson (as discussed all the calculations are valid for a discrete probability measure as well). As it can be seen, the new perturbed density f x iτ is still in the same family of distributions with the perturbed distribution parameters. For the case of normal and log-normal distributions when we study the effect of perturbation of the second component of the natural parameter, there is only one solution for τ: 1 -2τσ 2 > 0. This solution could not be found analytically but with the help of a numerical solver. For the exponential and the Poisson distributions the solutions for τ is expressed with the Lambert W function. This stands for the the multivalued inverse relation of the function f (w) = w exp(w), where w is complex. We denote by W 0 (x) the upper real branch of the Lambert function on the interval [-1/e, 0] and by W -1 (x) the lower real branch on the same interval.

We will analyze the effect of these perturbations on an analytical example in Section 6.

Distribution Natural Parameter Sufficient Statistics

New cumulant function τ 1 (δ) and Resulting distribution η(θ)

T(x) ψ(τ) τ 2 (δ) Normal N(µ, σ 2 ) µ/σ 2 -1/2σ 2 x x 2       µτ + τ 2 σ 2 2 -µ 2 τ 1-2τσ 2 + 1 2 log(1 -2τσ 2 )       τ 1,2 = ± √ 2δ σ
Numerical solver

N(µ + τσ 2 , σ 2 ) N( µ 1-2τσ 2 , σ 2 1-2τσ 2 ) LogNormal logN(µ, σ 2 ) µ/σ 2 -1/2σ 2 log x (log x) 2       µτ + τ 2 σ 2 2 -µ 2 τ 1-2τσ 2 + 1 2 log(1 -2τσ 2 )       τ 1,2 = ± √ 2δ σ Numerical solver log N(µ + τσ 2 , σ 2 ) log N( µ 1-2τσ 2 , σ 2 1-2τσ 2 ) Exponential Exp(λ) -λ x log λ λ-τ τ 1,2 (δ) = λ(W-1,0(-e -1-δ )+1) W-1,0(-e -1-δ ) Exp(λ -τ) Poisson Pois(λ) log(λ) x λ(e τ -1) τ 1,2 (δ) = W -1,0 -λ-δ eλ + 1 Pois(λ exp(τ))
Table 1: Resulting distributions.

Bounded support distribution

Now, let us consider x i ∼ U[a, b], the uniform distribution on the interval [a, b]. The density is expressed as:

f x i (x) = 1 b -a I x∈[a,b] (x), (b > a).
The uniform distribution does not belong to the family of exponential distributions. It has a limited and a bounded support [a, b]. If we apply the same perturbation as for exponential families, the normalization function becomes:

ψ(τ) = log e τb -e τa τ(b -a) , τ ∈ R. ( 10
)
Then the equation for τ(δ) is: τbe τbτae τae τb + e τa e τbe τa log e τbe τa τ(ba) = δ.

This equation has no explicit solutions for τ. The solutions can be found using a numeric solver. Suppose that τ * = τ(δ) is a solution of equation [START_REF] Melchers | Importance sampling in structural systems[END_REF]. Then, the perturbed density is:

f x iτ = τ * e τ * x e τ * b -e τ * a I x∈[a,b] (x).
Therefore, the new perturbed variable x iτ is no longer uniform on [a, b]. This density modification for a = -1, b = 1 and δ = 0.5 is displayed in Figure 1(a).

Notice that working with uncertain variables defined on a compact support, the main source of uncertainty is on the boundaries of the support. For such distributions with a bounded support, we propose to apply another density perturbation. The idea consists in perturbing the original boundaries by τ = ±δ. In the same way as with infinite support we consider the effect of positive or negative perturbation. For example, consider

f τ1 (x) = 1 b -a -δ I x∈[a+δ,b] (x) or f τ2 (x) = 1 b -a -δ I x∈[a,b-δ] (x).
The same perturbation may be applied to a triangular or a trapezoidal distribution. It can be also applied to the truncated Gaussian distribution if one is interested about the boundary influence on the failure probability p f . In this case the density function should be corrected for the new boundaries. Next, we explain how to estimate efficiently a perturbed failure probability p f δ with no additional CPU cost.

Inverse importance sampling and sensitivity analysis

Monte Carlo sampling is one of the most popular simulation methods to estimate a failure probability. We consider the input variables space Ω ∈ R d . Recall that all the input variables are independent and that f x( x) = d k=1 f x k (x k ) is the joint density of the input variables. Let X N = { x1 , . . . , xN } i.i.d.

∼ f x(•) be a sample of size N. Then, the estimation of the failure probability p f is given by:

p f = 1 N N k=1 I g( xk )≤0 . (11)
Now, assume that f x iτ is a new perturbed density for the input variable x i . Then the new joint density is

f xiδ = f x 1 • • • f x i-1 f x iτ f x i+1 • • • f x d = f x iτ (x i ) d k=1,k i f x k (x k ).
The corresponding failure probability p iδ is defined as an expectation of the indicator function:

p iδ = E f xiδ I g( x)≤0 = Ω I g( x)≤0 f xiδ ( x)d x
Here, we propose to apply the technique used in the Importance Sampling (IS) simulation method. We multiply the integrand function by 1 = f x ( x) f x ( x) . Both density functions f x( x) and f xiδ are the products of the density functions of the independent variables x = (x 1 , . . . , x d ) ∈ Ω ⊂ R d with the only difference for the variable x i . Therefore, we obtain:

p iδ = Ω I g( x)≤0 f xiδ ( x)d x = Ω I g( x)≤0 f x iτ (x i ) f x i (x i ) f x( x)d x = E f x I g( x)≤0 f x iτ (x i ) f x i (x i ) .
By doing so, we do not need to throw a new sample according to the unknown density function f xiδ ( x). We are working in the same probability space integrating the function I g( x)≤0

f x iτ (x i ) f x i (x i ) .
So that, to estimate the perturbed failure probability p iδ we keep the same sample points from the failure region: X N f = { x ∈ X N : g( x) ≤ 0} that provide non-zero values of the indicator function I g( x)≤0 . The estimation of the failure probability for the perturbed density is expressed as:

p iδ = 1 N N k=1 I g( x)≤0 f x iτ ( xki ) f x i ( xki ) . ( 12
)
If we are interested in the interaction effects, we perturb the probability densities for the variables x i and x j simultaneously. Then, the new joint density is:

f xi j,τ(δ 1 ),τ(δ 2 ) ( x) = f x iτ(δ 1 ) (x i ) f x jτ(δ 2 ) (x j ) d k=1,k i, j f x k (x k ).
Therefore, in this case the new failure probability can be estimated by:

p i j,δ 1 ,δ 2 = 1 N N k=1 I g( x)≤0 f x iτ(δ 1 ) ( xki ) f x jτ(δ 2 ) ( xk j ) f x i ( xki ) f x j ( xk j ) .
This describes the interaction effects of two variables x i and x j on the failure probability p f . The proposed reliability sensitivity analysis is based on the analysis of the value of the perturbed failure probability. In order to clearly differentiate the magnitude of the influence to the failure probability, we propose a sensitivity measure for the input variables. We provide the indices formulation in the next section and we study their statistical properties in Appendix B.

Sensitivity indices formulation

There are many possible choices of sensitivity indices. In this paper, we propose the sensitivity indices based on the difference p iδp f with the original failure probability p f . It is expressed as ratio:

S iδ = p iδ -p f p f (13)
Another possible formulation could be found in Lemaitre and Sergienko et. al, 2012 [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF].

The support of S iδ formulated in ( 13) is [-1, +∞). A negative value of S iδ means that the proposed density modification reduces the failure probability. Conversely, a positive value of this index means an increase in the failure probability. Zero value of S iδ means that the variable x i has no impact on the failure probability.

In practice, p f is estimated by p f with the Monte Carlo simulation method. In the same way, according to (12) p iδ is estimated by p iδ . The estimator of the indices S iδ can be expressed as:

S iδ = p iδ -p f p f .
This estimation provides an asymptotic unbiased estimation of S iδ . Moreover, according to the Central Limit Theorem (CLT) and the ∆ -method we have:

1 VAR f x S iδ S iδ -S iδ N→∞ → N(0, 1)
The proof and an expression of the asymptotic variance can be found in Appendix B.2.

The above basic indices are computed separately for every fixed value of the density perturbation δ. It allows to study the effect of perturbation by varying the value of δ. It can help in reliability design optimization by adjusting the distribution parameters of the input variables x = (x 1 , . . . , x d ) in order to achieve the lowest failure probability p f .

Analytical function example

To investigate the previous indices, we consider a linear function of three independent normally distributed variables:

x 1 , x 2 , x 3 ∼ N(0, 1): g(x 1 , x 2 , x 3 ) = 3 -0.1x 1 -0.5x 2 -1.0x 3 . ( 14
)
As a linear combination of independent Gaussian random variables, g( x) is distributed as N(µ, σ), where µ = 3 and σ 2 = 0.1 2 + 0.5 2 + 1.0 2 = 1.26. The original failure probability can be explicitly computed:

p f = 1 -Φ( µ σ ) = 3.7 × 10 -3 .
The estimation provided by a Monte Carlo sample of size N = 10 6 yields p f = 3.69 × 10 -3 . We use the same Monte Carlo sample to estimate the perturbed failure probabilities p iδ and the sensitivity indices S iδ . Recall, that the natural parameter of the normal distribution is given by the vector:

η(µ, σ) = µ/σ 2 -1/2σ 2 .
First, we consider the effect of the perturbation of the first component of natural parameter η. This perturbation refers as a shifting of the original mean of the distribution (see Table 1). Figure 2 depicts sensitivity indices calculated for δ ∈ [0, 1] with the positive (2(b)) and the negative (2(a)) values of perturbation parameter τ. It can be clearly observed that the highest impact on the failure probability is due to the variable x 3 and that the variable x 1 has the lowest impact for both cases of τ. Moreover, the higher the value of δ, the higher the influence of the variables. Now, we study the effect of the perturbation of the second component of the vector of natural parameters -1/2σ 2 . In this case, we perturb both the mean and the variance of the original distribution (see Table 1). We use a numerical solver to find the solutions for τ. In this case we consider the only possible effect of perturbation when the normalization function ψ(τ) is defined and 1 -2τσ 2 > 0. Figure 3 displays the estimated sensitivity indices in this case. We can observe the same relationship: the variable x 3 has the highest influence on the failure probability. However, the magnitude of the influence differs from the case where we only consider a mean shifting. The variable x 1 still has almost negligible impact.

CO2 storage case example

In this section, we consider a CO 2 storage reservoir simulation model. As discussed in the introduction, subsurface CO 2 storage is always associated with an excess reservoir pressure. One of the primary environmental risks is a pressure-driven leakage of CO 2 from the storage formation.

In order to assess the risk of CO 2 leakage through the cap rock we consider a synthetic reservoir model. The structure of the reservoir is reduced to its simplest expression. The model is made up of three zones (Figure ( 4)): a reservoir made of 10 layers, a cap-rock made up of 1 layer, a zone-to-surface composed of 9 layers. The XY size of the grid is set at 10 km total length representing 26x26x20 model grid. Each layer is 5m thick, including the cell above the cap-rock. The zone above the cap-rock (up to the surface) is currently set to 1 layer. The salinity of the water is 35gm/l. The temperature of the reservoir is set to 60C and the initial pressure is hydrostatic. The injection bottom rate is set to 10 6 tons/year. The fracture pressure is estimated by geomechanical experts to be P f racture = 122 bars. Exceeding this value during the injection can lead to a leakage. The simulation period is 55 years that include an injection period of 15 years followed by 40 years of storage. In this study we analyze the possibility of a leakage through a cap rock. Therefore, we consider pressure in the storage reservoir at the last year of injection as an objective function.

The uncertain variables selected for this study characterize the reservoir and the fluid properties. It implies different CO 2 flowing possibilities between the reservoir layers. Table (2 For sake of clarity we transform the original intervals into [-1, 1]. In this section, we assume the truncated standard normal distribution for all the input variables N [-1,1] (0, 1). The performance function for this example can be formulated. Suppose that P reservoir ( x) is a function of the reservoir pressure depending on the input variables configuration x ∈ Ω ⊂ R d . Then, the performance function defining the event of a gas leakage can be expressed as:

g( x) = P f racture -P reservoir ( x).
The reservoir pressure P reservoir ( x) is computed with a complex dynamic reservoir simulator. For this reason, we use the Gaussian process based response surface model approximation P( x) [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Welch | Screening, predicting, and computer experiments[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF]. By approximating the function of the reservoir pressure, we can quantify the risk and estimate the reliability of the system.

The original reference value of the failure probability computed by GP model approximation with Monte Carlo sample of size N = 10 6 provides an estimation: p f = 2.26 × 10 -4 . We keep this sample to estimate the perturbed failure probability p iδ and the sensitivity indices S iδ .

There are two possible ways to perturb the truncated Gaussian distribution. We can use the perturbation defined for the exponential family or we can study the effect of the perturbation of the distribution support boundaries. In this example, we compare the results for both cases. We start with the sensitivity indices calculated by the mean shifting. Figure 5 displays the evolution of the sensitivity indices for δ ∈ [0, 1] for negative (Figure 5(a)) and positive (Figure 5(b)) values of shifting. The variables ranking depends on the sign of τ. When τ > 0 (i.e. the positive mean shifting) the porosity variable has the highest impact on the failure probability. It means that increasing the mean value of the reservoir porosity PORO leads to increasing the failure probability. On the contrary, increasing the mean value of the reservoir permeability KSAND and the water relative permeability end-point KRSAND has a negative effect on the risk of leakage. For the negative mean shifting, the variables KSAND and KRSAND have the highest influence on the failure probability. Reducing the reservoir permeability and the end-point water relative permeability impedes the gas flow in the reservoir. It increases the risk of the leakage. When τ < 0 (i.e. the resulting distribution is N [-1,1-δ] (0, 1)) the porosity variable has the highest impact on the failure probability. It means that by decreasing the maximum value of the reservoir porosity PORO the failure probability decreases. It is also shown by Figure 5(a). It can be also observed that for this variable the new perturbed failure probability p iδ is equal to zero when δ > 0.3. When τ > 0 the resulting distribution is N [-1+δ,1] (0, 1). For this case, increasing the reservoir permeability KSAND and the water relative permeability end-point KRSAND reduces the failure probability p iδ . Both methods provide comprehensive and complementary results. If the main uncertainty is about the boundaries the one can start with the boundaries perturbation. By moving the boundaries of the original distribution, the one can determine the safe intervals for the input variables by detecting the value of δ: p iδ = 0 or S iδ = -1. After that, the effect of the mean perturbation can be studied.

Conclusions

In this paper, we have studied and adapted a recently introduced approach to the reliability sensitivity analysis. Currently the majority of the methods for reliability analysis is based on the variance decomposition and Sobol' sensitivity indices. We present a moment independent sensitivity measure. The method is based on a perturbation of the original probability distribution of the input random variables. We can analyze the a priori assumption about the input distributions and measure the effect of some possible deviations from this assumption. In particular, we select the Kullback-Leibler divergence as a measure of the perturbation.

We have provided different possible density perturbations with the resulting distributions for the exponential family of distributions. We have also studied the distributions with a bounded support and the effect of the boundaries perturbation.

Considering a proposed perturbation for an input variable, we present an effective method to estimate the corresponding perturbed failure probability. The method is based on a technique coming from the importance sampling simulation method. It allows to estimate the new failure probability without supplementary performance function evaluations. The new sensitivity indices formulation describes the relationship between the new failure probability p iδ and the original failure probability of the system p f . By varying the value of the perturbation δ, we can study how the positive or negative probability density perturbation affects the failure probability. If the model has controllable input variables, the method can help improving the system reliability and the design optimization. The presented analysis on the statistical properties of the proposed estimators for the perturbed failure probability p iδ and the sensitivity indices S iδ shows asymptotic normality of the estimators.

We investigated the method on an analytical and a CO 2 storage reservoir cases. The method provides promising results and can be applied in the reliability sensitivity analysis.

Appendix A.1. Derivatives • ψ (τ) = E[x τ ], τ ∈ D ψ (τ) = d dτ ∞ -∞ exp(τx) f (x)dx ∞ -∞ exp(τx) f (x)dx = ∞ -∞ x exp(τx -ψ(τ)) f (x)dx = ∞ -∞ x f τ (x)dx = E[x τ ] • ψ (τ) = E [x τ -E(x τ )] 2 , τ ∈ D ψ (τ) = d dτ ∞ -∞ x exp(τx -ψ(τ)) f (x)dx = = ∞ -∞ x(x -ψ (τ)) exp(τx -ψ(τ)) f (x)dx = = ∞ -∞ x 2 f τ (x)dx -ψ (τ) ∞ -∞ x f τ (x)dx = ∞ -∞ x 2 f τ (x)dx -[E[x τ ]] 2 = = E [x τ -E(x τ )] 2 = VAR(x τ ) • ψ (τ) = E [x τ -E(x τ )] 3 , τ ∈ D E [x τ -E(x τ )] 3 = ∞ -∞ x -ψ (τ) 3 f τ (x)dx = = ∞ -∞ x 3 f τ (x)dx -3ψ (τ) ∞ -∞ x 2 f τ (x)dx + 3ψ (τ) 2 ∞ -∞ x f τ (x)dx -ψ (τ) 3 = ∞ -∞ x 3 f τ (x)dx -3ψ (τ) ∞ -∞ x 2 f τ (x)dx + 2ψ (τ) 3 = ψ (τ) Appendix B.

Statistical properties of the indices estimator

Here, we will study statistical properties of the estimator of the perturbed failure probability p iδ and corresponding estimator of the sensitivity indices S iδ = p iδp f p f . We will start with studying the properties of p iδ . Appendix B.1. Estimator of the perturbed failure probability Suppose, f x (x) = d i=1 f x i (x i ) is the input joint density and f x iτ is a perturbed probability density for the variable x i . Recall that for a sample of size N: {x 1 , . . . , x N } i.i.d.

∼ f x (x), the estimation of p iδ is computed by:

p iδ = 1 N N k=1 I g(x k <0) f x iτ (x k i ) f x i (x k i ) .
First, we study the expectation and the variance of this estimator.

1. E f x p iδ = p iδ 2. VAR f x p iδ = 1 N VAR f x I g(x)<0 f x iτ (x i )

f x i (x i ) = 1 N I g(x)<0 f 2 x iτ (x i ) f 2 x i (x i ) f x (x)dx -p 2 iδ
This variance tends to 0 when N → ∞. Furthermore, by the Central Limit Theorem (CLT): 1 VAR f x p iδ p iδp iδ N→∞ → N(0, 1)

Note that the covariance between the estimator p f and p iδ does not vanish. Indeed, we use the same sample to estimate p f and p iδ . We can compute this covariance:

COV( p f , p iδ ) = E f x ( p f p iδ ) -E( p f )E( p iδ ) = = N N 2 I g(x k <0) f x iτ (x i ) f x i (x i ) f x i (x i )dx -p f p iδ = = 1 N p iδ (1 -p f )
The value of this covariance decreases when the sample size N increases.

Appendix B.2. Sensitivity indices estimator

Recall, that the first sensitivity index is:

S iδ = p iδ -p f p f = p iδ p f -1, i = 1, . . . , d.
We estimate this value with the Monte Carlo method and the importance sampling by estimating consistently p iδ and p f . The estimator of this index is:

S iδ = p iδ p f -1.
Here, we will study some proprieties of this estimator.

It is not straightforward to compute directly E f x S iδ = E f x p iδ p f -1 and VAR f x S iδ = VAR f x p iδ p f . We propose to use the Delta Method to approximate these values [START_REF] Van Der | Asymptotic statistics[END_REF].

Let us recall the Taylor expansion with integral form for the remainder. Let φ be a two times differentiable function on [t 0 , t], then: φ(t) := φ(t 0 ) + φ (t 0 )(tt 0 ) + t t 0 (1u)φ (u)du (B.1)

We will define a function:

φ(t) = y(t) x(t) -1,
where x(t) = (1t)p f + t p f and y(t) = (1t)p iδ + t p iδ . For this function: φ(0) = S iδ and φ(1) = S iδ . Following the Taylor expansion (B.1) we will expand φ(t) with t = 1 and t 0 = 0. First, we will compute the derivatives. Knowing the variance of the estimator, the confidence region for the indices may be computed.

  x i ∼ U[a, b] to be uniformly distributed on the interval [a, b]. Then, in order to stay inside the support the perturbed random variable x iτ is uniformly distributed either on U[a + δ, b] or on U[a, b -δ]. The corresponding density for perturbed uniform distribution can be expressed as:

Figure 1 :

 1 Figure 1: Uniform density perturbation.

Figure 2 :

 2 Figure 2: Basic sensitivity indices example. Mean shifting.

Figure 3 :

 3 Figure 3: Basic sensitivity indices example. Mean and variance shifting.

Figure 4 :

 4 Figure 4: Reservoir model.

  ) represents the variables description with their range of minimum and maximum values.

Figure 5 :

 5 Figure 5: Mean shifting. Now, we consider the boundaries perturbation. It means that we are moving one of the the distribution boundaries by δ in positive or negative directions keeping the values of mean and variance unchanged.Figure 6 depicts the sensitivity indices for δ ∈ [0, 1] for negative (Figure 6(a)) and positive (Figure 6(b)) values of τ = ±δ.

Figure 6

 6 depicts the sensitivity indices for δ ∈ [0, 1] for negative (Figure 6(a)) and positive (Figure 6(b)) values of τ = ±δ.

Figure 6 :

 6 Figure 6: Boundaries shifting.

p 3 f 1 p 2 fVAR f x p iδ - p 2 iδ ( 1 -p f ) N p 3 f. 1 VAR

 3122131 COV( p iδ , p f ) = Therefore, the variance of the indices estimator tends to 0 when N → ∞. With some extra computations we can show that: f x S iδ S iδ -S iδ N→∞ → N(0, 1).

Table 2 :

 2 Uncertain variables.

  p f ) p f p iδp f p iδ x 3 (t) . t)( p fp f ) p f p iδp f p iδ (1t)p f + t p f 3 dt = ( p fp f )(p f p iδp f p iδ ) p f )(p f p iδp f p iδ )Now, we can approximate the mean and the variance of S iδ .1. E f x S iδ ∼ S iδ 2. VAR f x S iδ ∼ VAR f x

				p iδ p f + VAR f x	p f p iδ f p 2	-2p iδ
	Then,			φ (t) =	y (t)x(t) -x (t)y(t) x 2 (t)	=	p f p iδ -p f p iδ (1 -t)p f + t p f	2
				φ (t) t=0 =	p f p iδ -p f p iδ p 2 f
	The second derivative is:
	φ (t) = 2( p f Therefore, the reminder is: 2x(t)x (t) p f p iδ -p f p iδ x 4 (t) = 2x (t) p f p iδ -p f p iδ x 3 (t) =
		0	1	2(1 p 2 f p f	.
	So that, by the Taylor expansion (B.1) we obtain:
				S iδ = S iδ +	p f p iδ -p f p iδ p 2 f	+	( p f p 2 f p f	.
	The last term R =	( p f -p f )(p f p iδ -p f p iδ ) p 2 f p f	is the remainder. This remainder is bounded and we can neglect it for the
	approximation.					S iδ ≈ S iδ +	p f p iδ -p f p iδ f p 2	.

1. x (t) = p fp f 2. y (t) = p iδp iδ

Acknowledgments

This work has been partially supported by the French National Research Agency (ANR) through COSINUS program (project COSTA-BRAVA n ANR-09-COSI-015). The authors also would like to thank Nicolas Maurand and Dan Bossie-Codreanu for the presented CO 2 reservoir model. We thank Sébastien Da-Veiga, Nicolas Bousquet and Bertrand Iooss for fruitful discussions.

Appendix A. Properties of the normalization function

For a given random variable x with probability density f (x). We propose the density modification x τ ∼ f τ (•) as follows:

, where ψ(τ) is a normalization function given by:

Let us define D = {τ ∈ R : ψ(τ) < +∞}. D defines the interior of D. We will also suppose, that ∃ε : D ⊃] -ε, ε[. Here, we will study the properties of this normalization function.


Ekaterina Sergienko

Fabrice Gamboa

Daniel Busby

Shape invariant model approach for functional data analysis in uncertainty and sensitivity studies

Keywords: Kriging, Functional data analysis, Semiparametric model

Dynamic simulators model systems evolving over time. Often, it operates iteratively over fixed number of time-steps. The output of such simulator can be considered as time series or discrete functional outputs. Metamodeling is an effective method to approximate demanding computer codes. Numerous metamodeling techniques are developed for simulators with a single output. Standard approach to model a dynamic simulator uses the same method also for multi-time series outputs: the metamodel is evaluated independently at every time step. This can be computationally demanding in case of large number of time steps. In some cases, simulator outputs for different combinations of input parameters have quite similar behaviour. In this paper, we propose an application of shape invariant model approach to model dynamic simulators. This model assumes a common pattern shape curve and curve-specific differences in amplitude and timing are modelled with linear transformations. We provide an efficient algorithm of transformation parameters estimation and subsequent prediction algorithm. The method was tested with a CO 2 storage reservoir case using an industrial commercial simulator and compared with a standard single step approach. The method provides satisfactory predictivity and it does not depend on the number of involved time steps.

Introduction

Simulation models are used nowadays in many industrial applications to predict and analyze the behaviour of complex systems. A simulator is a complex computer code embedding the physical laws governing the physical system under investigation. The input of such simulators can be adjustable or uncontrollable parameters which are only partially known and thus are affected by uncertainty. Uncertainty analysis of numerical experiments is used to assess the confidence of the model outcomes which are then used to make decisions [START_REF] De Rocquigny | Uncertainty in industrial practice: a guide to quantitative uncertainty management[END_REF]. Here, we focus on a particular type of dynamic simulators that are used to make predictions in the future. These simulators are based typically on finite element/finite difference codes used for instance for the simulation of flows and transfers in porous media. Industrial applications using this type of simulators are for instance hydrocarbons reservoir forecasting and carbon dioxide (CO 2 ) underground storage [START_REF] Sergienko | Optimal well placement for risk mitigation in co2 storage[END_REF][START_REF] Busby | Combining probabilistic inversion and multi-objective optimization for production development under uncertainty[END_REF]. Such applications involve very complex numerical codes with a large number of inputs and with high uncertainty derived from an incomplete knowledge of subsurface formation [START_REF] Subbey | Prediction under uncertainty in reservoir modeling[END_REF]. The uncertainty on the simulator output is usually assumed to be mostly due to the propagation of the uncertainty on the input. Nevertheless, modelling errors can also play a role.

The simulator models a multi-phase 3-D fluid flow in heterogeneous porous media, operating over fixed number of time-steps. The typical output of such simulators consists of a sequence of outputs at different time-steps. Therefore, it represents time series related, for instance, to a recovery rate for a given well or a group of wells. It can be also a spatial output such as a pressure map or a saturation map also evolving with time. Here we focus on 1D time series output which can be typically measured in a well.

Let us consider the output of a simulator as a deterministic function Y(t) = F( x, t), where x ∈ Ω ⊂ R d is a configuration of preselected input parameters and 0 < t < T refers to a time step. Y(t) is a time dependent output, e.g. oil production rate or reservoir pressure.

The function F : R d × [0, T ] → R models the relationship between the input and the output of the numerical simulator. Methods such as Monte Carlo ones can be used to propagate uncertainty. However, each model evaluation being generally very time-consuming this approach can be unpractical in industrial applications. Moreover, the probability distribution of the input may vary as new data comes in or within a process of dynamic data assimilation. Therefore, a response surface approach (also referred to a metamodel approach) can be more useful in such applications [START_REF] Busby | Combining probabilistic inversion and multi-objective optimization for production development under uncertainty[END_REF][START_REF] Busby | Adaptive design of experiments for bayesian inversion an application to uncertainty quantification of a mature oil field[END_REF][START_REF] Busby | Hierarchical nonlinear approximation for experimental design and statistical data fitting[END_REF].

The advantage of a metamodel is that is fast to evaluate and it is designed to approximate the complex computer code based on a limited number of simulator evaluations. These evaluations of the simulator are taken at some well chosen input configurations also called the training set or the experimental design. Numerous experimental designs have been proposed and many of them are quite sophisticated. Here, we use Latin Hypercube designs [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] for their good space filling properties. Usually, it can be coupled with other criteria such as the maximin design (maximum minimum distance) [START_REF] Sacks | The design and analysis of computer experiments[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. In this work, we focus on Gaussian process (GP) based metamodels also known as kriging [START_REF] Matheron | Principles of geostatistics[END_REF][START_REF] Sacks | The design and analysis of computer experiments[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

The aim of this work is to propose a new method to address multi-time series outputs. For such dynamic simulator a standard approach assumes a single step GP model for each time step. This basic approach can be computationally intensive for a high number of time steps and when the construction of a single GP model is time consuming (i.e. when the size of the training set and the number of variables is large). Therefore, the procedure can become unpractical in some industrial applications.

The problem of metamodeling for dynamic simulators was recently investigated by different authors and several principal approaches can be distinguished. A first possible approach is GP metamodelling considering the time steps as the model additional input parameter [START_REF] Busby | Combining probabilistic inversion and multi-objective optimization for production development under uncertainty[END_REF][START_REF] Conti | Bayesian emulation of complex multi-output and dynamic computer models[END_REF][START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF]. This approach is easy to implement, however if we want to take into account all the information from an experimental design at every time steps, the size of new experimental design is multiplied by the number of time steps. It results to matrices of high dimensions that can lead to some numerical problems in case of large size of original experimental design or in case of high density of the steps in the time scale. [START_REF] Conti | Gaussian process emulation of dynamic computer codes[END_REF] [START_REF] Conti | Gaussian process emulation of dynamic computer codes[END_REF] developed an iterative approach to build a model of dynamic computer code, assuming that the model output at a given time step depends only on the output at the previous time step. To reduce the dimensionality of the problem, we can also represent the given observations as truncated functions by some selected basis functions with the following GP modelling for the coefficient of the selected representation. [START_REF] Bayarri | Computer model validation with functional output[END_REF] [START_REF] Bayarri | Computer model validation with functional output[END_REF] introduced wavelet decomposition. [START_REF] Campbell | Sensitivity analysis when model outputs are functions[END_REF] [START_REF] Campbell | Sensitivity analysis when model outputs are functions[END_REF], Higdon et al. (2008) [START_REF] Higdon | Computer model calibration using high-dimensional output[END_REF] and [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] suggested application of principal component decomposition. Auder et al. (2010) [START_REF] Auder | Screening and metamodeling of computer experiments with functional outputs. application to thermal-hydraulic computations[END_REF] extended this approach by preliminary classification.

In this work, we propose a new functional approach involving a combination of Shape Invariant Model (SIM) approach and Gaussian Process modelling. Considering J time steps and a training set X n = { x1 , .., xn } the simulator output consist then in a set of curves:

Y n = Y i, j = F( xi , t j ), 1 ≤ i ≤ n, 1 ≤ j ≤ J .
In our practical example (the CO 2 storage case study) we have observed that usually these curves have quite similar behaviour. So that, we assume that there is a mean pattern from which all the curves can be deduced by a proper parametrical transformation. The shape invariant model assumes that the curves have a common shape, which is modelled nonparametrically. Then curve-specific differences in amplitude and timing are modelled with the linear transformations such as shift and scale [START_REF] Lawton | Self modeling nonlinear regression[END_REF][START_REF] Brumback | Self modeling with flexible, random time transformations[END_REF][START_REF] Izem | Analysis of nonlinear modes of variation for functional data[END_REF]. Hence, we consider the following model:

Y i, j = α * i f (t j -θ * i ) + v * i + σ * i ε i j , 1 ≤ i ≤ n, 1 ≤ j ≤ J (1)
where t j stands for observation times and

θ * = {θ i , 1 ≤ i ≤ n}, v * = {v i , 1 ≤ i ≤ n}, α * = {α i , 1 ≤ i ≤ n} ∈ R n
are vectors of transformation parameters corresponding to horizontal and vertical shifts and scales of the unknown function f . Without loss of generality we assume that the simulator provides the data at some constant time intervals. Time t is then equispaced in [0, T] with J time steps. The unknown errors σ * i ε i, j are independent zero-mean random variables with variance σ * i 2 . Here, we can also assume that σ * i 2 = 1 without loss of generality. The approach proposed in this work for the functional outputs modeling combines an efficient estimation of the transformation parameters (α * , θ * , v * ) and a subsequent GP modeling for these parameters that can be used to predict new curves without running the simulator.

The paper is organized as follows. Section 2 presents the basics of GP modeling and the model validation criteria. Section 3 describes the method of efficient estimation of transformation parameters in the shape invariant model. The method is illustrated with an example on an analytical function. In Section 4 we present the forecast algorithm for dynamic simulators. Section 5 provides the practical application of the algorithm with a CO 2 storage reservoir case.

GP based metamodeling

The method proposed in this paper is based on a combination of a shape invariant model and a Gaussian process metamodel. In this section, we recall the basics of GP modeling or kriging.

The idea of modeling an unknown function by a stochastic process was introduced in the field of geostatistics by Krige in the 1950's [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF] and formalized in 1960's by [START_REF] Matheron | Principles of geostatistics[END_REF] [START_REF] Matheron | Principles of geostatistics[END_REF]. Later [START_REF] Sacks | The design and analysis of computer experiments[END_REF] [START_REF] Sacks | The design and analysis of computer experiments[END_REF] proposed the use of kriging for prediction and design of experiments. The theory and the algorithms are formalized in [START_REF] Sacks | The design and analysis of computer experiments[END_REF], [START_REF] Welch | Screening, predicting, and computer experiments[END_REF] and [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

Consider the output of a simulator as an unknown deterministic function F( x) ∈ R, where x ∈ Ω ⊂ R d is a specified set of selected input parameters. The function F is only known in predetermined design points:

X n = {( xk , F k = F( xk )) , 1 ≤ k ≤ n}.
The objective is to predict the function F 0 = F( x0 ) for some new arbitrary input x0 . The function is modeled as a sample path of a stochastic process of the form:

F( x) = m j=1 h j ( x) • β j + Z( x) = β h( x) + Z( x) (2)
where:

• β h( x) is the mean of the process and corresponds to a linear regression model with preselected given realvalued functions h = {h i , 1 ≤ i ≤ m}. Here, we only consider the case h = 1.

• Z( x) is a centered Gaussian stationary random process. It is defined by its covariance function: Numerous families of correlation functions have been proposed in the literature. We use here Gaussian correlation function, the special case of the power exponential family. The power exponential correlation function is of the following form:

C( x, ȳ) = E Z( x)Z(ȳ) = σ 2 R( x, ȳ). R( x,
R( x, ȳ) = exp          - d j=1 x j -y j p j θ j          = exp         - d j=1 d p j j θ j         (3)
where d j = |x jy j |, 0 < p j ≤ 2 and θ j > 0. The hyperparameters (θ 1 , .., θ d ) stands for correlation lengths which affect how far a sample point's influence extends. A high θ i means that all points will have a high correlation (F(x i ) being similar across our sample), while a low θ i means that there are significant difference between the F(x i )'s [START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF]. The parameters p j corresponds to the smoothness parameters. These parameters determine mean-square differentiablity of the random process Z(x). For p j = 2 the process is infinitely mean-square differentiable and the correlation function is called Gaussian correlation function. Hence, Gaussian correlation function is infinitely mean-square differentiable and it leads to a stationary and anisotropic process Z(x) [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF][START_REF] Welch | Screening, predicting, and computer experiments[END_REF]. Regardless the choice of a correlation function, the estimation of hyperparameters (θ 1 , .., θ d ) is crucial for reliable prediction. We are using maximum likelihood estimation algorithm [9] that we will discuss later.

The experimental design points are selected in order to retrieve most information on the function at the lowest computational cost. The number of design points for a reliable response surface model depends on the number of inputs and on the complexity of the response to analyze [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Latin Hypercube Designs (LHD) provides a uniform coverage of the input domain. If we wish to generate a sample of size n, first, we partition the domain of each variable in n intervals of equal probability. Then, we randomly sample x1 according to the distribution of each of the n intervals. Further, for each of the n values for x1 , we randomly select one interval to sample for x2 , so that only one sample of x2 is taken in each interval. We continue the process of a random sampling without replacement until all the variables have been sampled. As a result we generate a sample where each of d inputs is sampled only once in each of N intervals. Latin hypercube designs have been applied in many computer experiments since they were proposed by [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF].

In this work, we use modified version of LHD -maximin LHD. It is based on maximizing a measure of closeness of the points in a design

D n : max design D n min x1 , x2 ∈D n d( x1 , x2
) It can guarantee that any two points in the design are not "too close". Hence, the design points are uniformly spread over the input domain.

Consequently, when we have the experimental design X n = ( x1 , .., xn ) and the observation data Y n = (F( x1 ), .., F( xn )) the multivariate distribution according to the model (2) for the Gaussian correlation function can be expressed as:

Y 0 Y n ∼ N 1+n h ( x0 ) H β, σ 2 1 r ( x0 ) r( x0 ) R , where R = (R( xi , x j )) 1≤i, j≤n ∈ R n×n is the correlation matrix among the observations; r( x0 ) = (R( x1 , x0 ), .., R( xn , x0 )) ∈
R n is the correlation vector between the observations and the prediction point;

h ( x0 ) = h j ( x0 ) 1≤ j≤m
∈ R m is the vector of regression function at the prediction point x0 and H = h j ( xi )

1≤i≤n,1≤ j≤m ∈ R n×m is the matrix of regression functions at the experimental design. The parameters β and σ are unknown.

Considering the unbiasedness constraint, the parameter β is replaced by the generalized least squares estimate β in (2). Here, β is of the following form:

β = H R-1 H -1 H R-1 Y n .
Assuming that the correlation function is the Gaussian correlation function, the prediction is therefore given by:

F( x0 ) = h ( x0 ) • β + r( x0 ) R-1 F n -H • β (4)
The hyperparameters θ = (θ 1 , .., θ m ) and the process variance σ 2 are estimated by Maximum Likelihood (MLE). Using the multivariate normal assumption, the MLE for σ 2 is:

σ2 = 1 m Y n -H β R -1 Y n -H β (5)
Knowing the estimations for β and σ2 , the coefficients θ are estimated by maximizing the log likelihood:

l β, σ2 , θ = - 1 2 m log σ2 (θ) + log(det(R(θ)) + m (6)
The function [START_REF] Busby | Hierarchical nonlinear approximation for experimental design and statistical data fitting[END_REF] depends only on θ.

After having estimated all the model parameters, we now need to validate the model. In this work for estimation of prediction accuracy of the model, we use the predictivity index, Q 2 , and root mean squared error, RMS E . The predictivity index is calculated basing on cross validation and it has the following form:

RMS E := n i=1 ( Ŝ X/x i -F(x i )) 2 n Q 2 := 1 - n i=1 Ŝ X/x i -F(x i ) 2 n i=1 F(x i ) -F 2 (7)
Ŝ X/x i is the kriging model computed using all the design points X n excepting xi and F is the mean of F( xi ).

The closer Q 2 is to 1 or RMS E is to 0, the higher is the model predictivity. These criteria can be also calculated on a separate validation test data by performing additional simulations. It provides higher accuracy measure though it requires additional time costs.

Shape Invariant Model

In this section, we discuss the shape invariant model representation and the procedure for efficient parameters estimation.

The shape invariant model was introduced by Lawton et al. (1972) [START_REF] Lawton | Self modeling nonlinear regression[END_REF]. The model assumes that we are working with a set of curves that have a common shape function that is modeled nonparametrically. The deformation of this function is modeled parametrically by choosing a proper parametrical transformation. We consider a class of linear transformations only. These parameters can be normally interpreted as shift in timing or scale in amplitude. For this reason, shape invariant model is widely applied to study periodic data such as temperature annual cycle [START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF] or traffic data analysis [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF]. Indeed, in these cases there is always some variability in time cycles or amplitude. The model has also been used to study biological data , where the departure from the pattern can be caused by a different environmental conditions [START_REF] Izem | Analysis of nonlinear modes of variation for functional data[END_REF][START_REF] Brumback | Self modeling with flexible, random time transformations[END_REF]. In this work, we use the model to propagate uncertainty on a reservoir simulator output data. For example, we consider a reservoir pressure during CO 2 injection, then shift in time is caused by different moment of stopping injection. So that, by applying shape invariant model, we can study the influence of the model input parameters on the overall shapes of the selected output.

In figure [START_REF] De Rocquigny | Uncertainty in industrial practice: a guide to quantitative uncertainty management[END_REF] we display three possible transformations that we consider in our work: horizontal shift 1(a), vertical shift 1(b) and vertical scaling 1(c). The bold line represent original pattern shape. We are interested in the common shape as well as in the efficient estimation of transformation parameters. So that, we can reproduce any curve and align it to the pattern. Moreover, we can make a prediction of a possible new curve for an input configuration x0 by modeling the transformation parameters.

As already mentioned, for an experimental design X n = { x1 , .., xn } we have a set of observations:

Y n = Y i, j = F( xi , t j ), 1 ≤ i ≤ n, 1 ≤ j ≤ J
, where xi , 1 ≤ i ≤ n is the set of preselected input parameters and t j , 1 ≤ j ≤ J refers to the time sample. So that, Y i, j is the j th observation on the i th experimental design unit, with 1 ≤ i ≤ n and 1 ≤ j ≤ J. Thereby, the idea is to find a general pattern curve that makes possible subsequent transformation of any curve to this selected pattern with properly adjusted transformation parameters.

We focus here on linear transformations. Thus, the model structure may be written as:

Y i, j = α * i • f (t j -θ * i ) + v * i + σ * 2 • ε i j ( 8
)
where t j , 1 ≤ j ≤ J are observation times which are assumed to be known and equispaced in the interval [0, T ]. The vector of parameters: (α * , θ * , v * ) = (α 1 , .., α n , θ 1 , .., θ n , v 1 , .., v n ) is unknown as well as the pattern function f . The errors ε i j are i.i.d. with a normal distribution, (i, j) ∈ {1, .., n} × {1, .., J}. It characterizes observation noise. Without loss of generality we may assume that σ * 2 = 1. The variance does not affect the parameters estimation procedure and the method still works for a variance σ * 2 . The function f is assumed to be 2π-periodic [START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF][START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF].

In this section we will provide the algorithm for efficient estimation of transformation parameters under unknown function pattern f . Since the functional pattern f is unknown, the pattern is replaced by its estimate. So that, it seems natural to study the problem (8) in a semi-parametric framework: the transformation shifts and scales are the parameters to be estimated, while the pattern stands for an unknown nuisance functional parameter. We use an Mestimator built on the Fourier series of the data. Under identifiability assumptions it is possible to provide a consistent algorithm to estimate (α * , θ * , v * ) when f is unknown. The algorithm is described in details in the following subsection with an illustration on an analytical function example.

Model Assumptions

Consider Y n = F(X n , t) = {F( xi , t j ), 1 ≤ i ≤ n, 1 ≤ j ≤ J} is (n × J) matrix of observations. We model these observations in the following way:

Y i, j = α * i • f (t j -θ * i ) + v * i + ε i j , 1 ≤ i ≤ n, 1 ≤ j ≤ J (9)
where f : R → R is an unknown 2π-periodic continuous function,

θ * = (θ * 1 , ..., θ * n ), α * = (α * 1 , ..., α * n ), v * = (v 1 , ..., v n ) ∈ R n are
unknown parameters, ε i j is a Gaussian white noise with variance equal to 1. The time period is translated in such a way that: [0, T[→ [0, 2π[, therefore t i = j-1 J 2π, j = 1, .., J are equispaced in [0, 2π[. The objective is to estimate the horizontal shift θ * = (θ * 1 , ..., θ * n ), the vertical shift v * = (v 1 , ..., v n ) and the scale parameter α * = (α * 1 , ..., α * n ) without knowing of the pattern f . Fourier analysis is well suited for the selected structure of the model. Indeed, this transformation is linear and shift invariant. Therefore, by applying a discrete Fourier transform to (9) and supposing J is odd, the model becomes:

d kl = α * k e -ilθ * k c l ( f ) + w kl , 1 ≤ k ≤ n, 0 < |l| ≤ (J -1)/2 α * k c 0 ( f ) + v * k + w k0 , 1 ≤ k ≤ n, l = 0 ( 10
)
where c l ( f ) = 1 J J m=1 f (t m )e -2iπ ml J , (|l|≤(J-1)/2) are the discrete Fourier coefficients and w kl , (1≤k≤n, |l|≤(J-1)/2) is a complex white Gaussian noise with independent real and imaginary parts.

We also notice that in order to ensure the identifiability of the model [START_REF] Matheron | Principles of geostatistics[END_REF] we are working in the parameter space:

A = (α * , θ * , v * ) ∈ [-π, π[ 3×n : α 1 = 1, θ 1 = 0, v 1 = 0 .
To summarize, in this section we estimate the transformation parameters (α * , θ * , v * ) without prior knowledge of the function f . The estimation depends on the unknown functional parameter (c l ( f )) |l|≤(J-1)/2 , the Fourier coefficients of the unknown function f . So that, we consider a semi-parametrical method based on an M-estimation procedure. M-estimation theory enables to build an estimator defined as a minimiser of a well-tailored empirical criterion that is given in the following subsection.

Parameters estimation procedure

The goal is to estimate the vector of parameters (α * , θ * , v * ) that depends on the Fourier coefficients of the unknown function f . We consider a semi-parametric method based on an M-estimation procedure [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF].

To construct an M-estimator, we define the rephased (untranslated and rescaled) coefficients for any vector (θ, α, v) ∈ A:

ckl (α, θ, v) =      1 α k e ilθ k d kl , 1 ≤ k ≤ n, 0 < |l| ≤ (J -1)/2 1 α k (d kl -v k ) , 1 ≤ k ≤ n, l = 0
and the mean of these coefficients:

ĉl (α, θ, v) = 1 n n k=1 ckl (α, θ, v), |l| ≤ (J -1)/2
Therefore, for (α * , θ * , v * ) we obtain:

ckl (α * , θ * , v * ) = c l ( f ) + 1 α * k e ilθ * k w kl 1 ≤ k ≤ n ĉl (α * , θ * , v * ) = c l ( f ) + 1 n n k=1 e ilθ * k • w kl α * k . Hence, |c kl (α, θ, v) -ĉl (α, θ, v)| 2 should be small when (α, θ, v) is closed to (α * , θ * , v * ).
Now, consider a bounded positive measure µ on [0, T ] and set

δ l := T 0 exp 2iπl T ω dµ(ω), (l ∈ Z) (11)
Obviously, the sequence (δ l ) is bounded. Without loss of generality we will assume that δ 0 = 0 and |δ l | > 0, l 0. Assume further that l |δ l | 2 |c l ( f )| 2 < +∞. So that f * µ is a well defined square integrable function:

f * µ(x) = f (x -y)dµ(y)
We construct the following empirical contrast function [START_REF] Qian | Gaussian process models for computer experiments with qualitative and quantitative factors[END_REF]:

M n (α, θ, v) = 1 n • n k=1 J-1 2 l=-J-1 2 |δ l | 2 |c kl (α, θ, v) -ĉl (α, θ, v)| 2 (12)
The random function M n is non negative. Furthermore, its minimum value is reached closely to the true parameters (α * , θ * , v * ). We define the estimator by:

( α, θ, v) n = arg min α,θ,v∈A M n (α, θ, v)
The proof of convergence ( α, θ, v) n P -→ n→+∞ (α * , θ * , v * ) and asymptotic normality of the estimators can be found in [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF][START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF]. The weights δ l in [START_REF] Conti | Bayesian emulation of complex multi-output and dynamic computer models[END_REF] are chosen to guarantee the convergence of the contrast function to a deterministic function M n and to provide the asymptotic normality of the estimators. Moreover, the convergence can be speeded up by proper selection of weights. The analysis of covergence at different weights is presented in [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF]. In this work we used the weight δ l = 1/|l| β with β = 1.5.

The computation of the estimators is very fast since only a Fast Fourier algorithm and a minimization algorithm of a quadratic functional are needed. The procedure is summarized in Algorithm ( 1

M n (α, θ, v) = 1 n • n k=1 J-1 2 l=-J-1 2 |δ l | 2 |c kl (α, θ, v) -ĉl (α, θ, v)| 2 Compute ( α, θ, v) = arg min α,θ,v∈A M n (α, θ, v) ∈ R 3×(n-1) Return: ( α, θ, v) ∈ R 3×(n-1)

Analytical Function Example

In this section the shape invariant model and the efficient parameters estimation are presented on an analytical function. The minimization algorithm used in the estimation procedure is a Newton-type algorithm.

Let us consider the following function:

f (x) = 20 • (1 -x/(2π)) • x/(2π)
Simulated data are generated as follows:

Y i j = α * i • f (t j -θ * i ) + v * i + ε i j , 1 ≤ i ≤ n, 1 ≤ j ≤ J
with the following choice of parameters: J = 5, N = 101, t j = j-1 J 2π, 1 ≤ j ≤ J are equally spaced points on [0, 2π[. Transformation parameters (θ * , α * , v * ) are uniformly distributed on ]0, 1], where θ * 1 = 0, v * 1 = 0, α * 1 = 1; the noise ε ji , j = 1, .., J, i = 1, .., N are simulated with a Gaussian law with mean 0 and variance 0.5.

Results are displayed in Figure 2. The function f is plotted by a solid red line in Figure 2(d). This analytical example shows that the method is effective in estimating the transformation parameters of the shape invariant model. In the next section, we will explain how this model can be applied in reservoir engineering forecast problems.

Functional data approximation

To apply the shape invariant model approach to approximate functional data from a dynamic simulator, firstly we have to modify the parameters estimation procedure for large number of curves. When we are working with uncertainty modeling, we always start from an experimental design X n and a set of observation Y n . As we have mentioned, the number of design points depends on the number of inputs and on the complexity of the response. So that, some optimization problems could arise when we compute the contrast function with a large number of curves. It can be time consuming and the results can be inaccurate. Therefore, we propose the following modification to the Algorithm (1): the original observation data is split into blocks and the optimization procedure is then performed on each of the blocks. Following the identifiability condition, we are working on a compact set

A = (α * , θ * , v * ) ∈ [-π, π[ 3×n : α 1 = 1, θ 1 = 0, v 1 = 0 .
This condition should be satisfied on every block optimization by adding the first reference curve to the block. With this procedure, we do not have limitations on experimental design of any size. As soon as we have estimated the parameters for every curve from observation data set, we can formulate the prediction algorithm. Instead of reproducing the simulator output for a prediction point x0 at every time step, we model the whole output curve with the transformation parameter. This curve will provide the approximation of the output for the selected input configuration x0 at each of considered time steps {t j , j = 1, .., J}. The transformation parameters for the input x0 are evaluated with the Gaussian process response surface modeling. The model is based on the experimental design and the set of evaluated transformation parameters calculated for the observation data curves. The prediction framework for an arbitrary input configuration x0 is presented by the following Algorithm [START_REF] Busby | Combining probabilistic inversion and multi-objective optimization for production development under uncertainty[END_REF].

Summing up, with this proposed algorithm the problem of response surface modeling for dynamic simulators is reduced from single step GP modeling for each of J time steps to an optimization problem and a GP modeling for the transformation parameters. However, it is worth to mention that before performing the algorithm it is important to analyze the curves behaviour for the observation data set. Studying the curves characterization, probably we may fix for example vertical shifts v or horizontal shifts θ at zero. Also, if the curves have significantly different behaviour at different time intervals we can split the observation data in time as well in order to achieve higher prediction accuracy.

The next section presents an application with a dynamic reservoir simulator case.

CO 2 storage reservoir case

Carbon Capture and Storage technology stands for the collection of CO 2 from industrial sources and its injection underground. Carbon dioxide is stored in a deep geological formation that is sealed on a top by a low permeability cap rock formation. Subsurface storage of CO 2 is always associated with an excess pressure in the reservoir and one Algorithm 3 Prediction algorithm for dynamic simulator Input: Dynamic simulator Y = F( x, t) with t ∈ {t j , j = 1, .., J} and prediction point x0 Output: Prediction Y 0 = F( x0 , t j ) for all j = 1, .., J Generate an experimental design X n = ( x1 , .., xn ) to span the space of interest Evaluate Y n = F(X n , t j ) at every time step t j , 1 ≤ j ≤ J Generate a set of discrete curves {Y i, j }, i = 1, .., n; j = 1, .., J Estimate the (α, θ, v) ∈ (R n ) 3 with Algorithm 2 Construct new experimental design for the function of parameters:

(X n , θ(X n )), (X n , α(X n )) and (X n , v(X n ))
Estimation of hyperparameters for GP models of transformation parameters α( x0 ), θ( x0 ) and v( x0 ) are approximated with corresponding GP models Reproduce:

F( x0 , t j ) = α( x0 ) f (t j -θ( x0 )) + v( x0
) for all {t j , j = 1, .., J} Return: Discrete time series Y 0 = F( x0 , t) with t ∈ {t j , j = 1, .., J} of the primary environmental risks is a pressure-driven leakage of CO 2 from the storage formation. In order to assess the risk of CO 2 leakage through the cap rock we consider a synthetic reservoir model. The model is made up of three zones (3):

• a reservoir made of 10 layers

• a cap-rock made up of 9 layers

• a zone-to-surface composed of 1 layer The XY size of the grid is set at 10 km total length. Each layer is 5m thick, including the cell above the cap-rock. The total number of cells is 13520 (26x26x20 model grid). The structure of the reservoir is reduced to its simplest expression. The zone above the cap-rock (up to the surface) is currently set to 1 layer. The salinity of the water is 35gm/l. The temperature of the reservoir is set to 60C and the initial pressure is hydrostatic. The injection bottom rate is set to 10E+06 tons/year. The fracture pressure is estimated by geomechanical experts to 158 bars. Exceeding this value of reservoir pressure can lead to a leakage. The simulation period is 55 years that include an injection period of 15 years followed by 40 years of storage. In this study we analyze the possibility of leakage through a cap rock. Therefore, we consider pressure in the storage reservoir as an objective function to be approximated.

10

The uncertain parameters selected for this study characterize the reservoir and the fluid properties. It implies different CO 2 flowing possibilities between the reservoir layers. The distribution law for the parameters is uniform. As we are sure here that the model parameters are efficiently estimated, we can proceed with the next step of prediction algorithm (3). The next step is to build Gaussian process response surface models for the transformation parameters basing on the estimations.

In figure [START_REF] Busby | Adaptive design of experiments for bayesian inversion an application to uncertainty quantification of a mature oil field[END_REF] we display the model validation criteria: Root Mean Square Error (RMSE) and predictivity indices (Q2), calculated separately for every year. The criteria were computed with the help of additional confirmation test data. The low predictivity in the first and last years is caused by low variance of data in that period. In general, the method provides reliable level of predictivity. It is also reflected by crossplots of test and predicted data. Table [START_REF] Sergienko | Optimal well placement for risk mitigation in co2 storage[END_REF] compares the simulation CPU time for both methods. SIM method provides suffucuently accurate results (although less accurate than the single step GP approach), but with a CPU reduction of a factor five.

SIM method

Single step GP Optimization Parameters modeling Total CPU time 00:00:35 00:00:15 00:00:50 00:04:28 It is worth to mention, that in this study we consider a simple model with only 3 uncertain parameters. So that, to estimate the function with GP model at every single step takes approximately 10 seconds. For more complex functions and more input uncertain parameters involved a single step model evaluation can take up to 10-20 minutes. So that, for the same simulation period of 55 years CPU time can increase to 10-20 hours. Whereas SIM approach does not depend on number of time steps and the method always conclude only a single optimization problem and as maximum 3 GP models for the transformation parameters.

Conclusion

This paper focuses on two general problems. First, we introduce a shape invariant model approach and we provide an efficient algorithm for estimation of transformation parameters of the model with the method specification for large sets of curves. Second, we suggest the application of this approach to model the time series outputs from a dynamic simulator. The proposed method reduces the problem of functional outputs modeling to one optimization problem and three GP response surface models. We have tested the method with a CO 2 storage reservoir case. We have also compared the method with the standard single-step approach. Presented numerical results show that the method provides satisfactory and comparable predictivity at lesser CPU time. The method also does not depend on the number of the involved time-steps. It can be very advantageous when we are working with a model involving large number of times steps such as CO 2 storage when the reservoir model simulation period can include up to hundreds or thousands time steps. However, if the set of output curves have significantly different behaviour, preliminary curves classification may be required.

  ȳ) is the correlation function and σ 2 = E[Z( x) 2 ] denotes the process variance. Stationarity condition assumes: R( x, ȳ) = R(| x -ȳ|), where | x -ȳ| denotes the distance between x ∈ Ω and ȳ ∈ Ω.

Figure 1 :

 1 Figure 1: Parametrical transformation examples

) Algorithm 1

 1 Parameters estimation procedure Input: Input set of curves from experimental design Y n = {Y i j , i = 1, .., n; j = 1, .., J} Output: Transformation parameters estimation (α * , θ * , v * ) Define the identifiability condition: A = (α * , θ * , v * ) ∈ [-π, π[ 3×n : α 1 = 1, θ 1 = 0, v 1 = 0 Compute the matrix of discrete Fourier coefficients D = {d kl , k = 1, .., n; |l| ≤ (J -1)/2} Compute the matrix of rephased Fourier coefficient C = {c kl , k = 1, .., n; |l| ≤ (J -1)/2} Compute the vector of mean of rephased coefficients C = { c l , |l| ≤ (J -1)/2} Choose the weight sequence δ l Define

Figure 2 (

 2 a) shows the original simulated noisy data Y i, j . The cross-sectional mean function of these data is presented in Figure 2(d) by black dotted line. Figure 2(b) plots estimated transformation parameter versus the originally simulated parameters. As it can be seen, the estimations are very close to the original parameters. The inverse transformation using the estimated parameters is displayed in Figure 2(c) and the mean function of restored curves is displayed in Figure 2(d) by blue dashed line. Figure 2(d) compares the cross-sectional mean of inversely transformed data and the cross-sectional mean of originally simulated data. Despite the noise, it is noticeable that the data after inverse transformation are much more closer to the original function f than the original cross-sectional mean function.

Figure 2 :

 2 Figure 2: Analytical Example

Algorithm 2

 2 Parameters estimation procedure for large n Input: Input set of curves from experimental design Y n = {Y i j , i = 1, .., n; j = 1, .., J} Output: Transformation parameters estimation (α * , θ * , v * ) Split the observation data into N b blocks of (K + 1) curves for m = 1, .., N b do Define block curvesY K+1 = {Y 1 , Y (m-1)(K+1)...,Y mK } = {Y i j , i = 1, .., K + 1, j = 1, .., J} Perform Algorithm (1) Compute ( α, θ, v) = arg min α,θ,v∈A M n (α, θ, v), where ( α, θ, v) ∈ R 3×Kend for Return: ( α, θ, v) ∈ R 3×(n-1)

Figure 3 :

 3 Figure 3: Reservoir Model

  We start from the observation data Y n -[30; 55] matrix of simulator outputs. By means of Algorithm (2) we provide the transformations parameters estimations. Figure (4(a)) provides the original set of curves and Figure (4(b)) represents the same set after inverse transformation. The pattern curve is differentiable after the inverse transformation with the estimated parameters.

Figure 4 :

 4 Figure 4: Original observation data and data after inverse transformation

  Figure (6(a)) is based on new proposed approach and Figure (6(b)) corresponds to single step GP modelling. Both methods provide a good level of approximation.

Figure 5 :Figure 6 :

 56 Figure 5: Predictivity Indices and RMSE

Table ( 1

 ( ) represents the parameters description with their range of minimum and maximum values.

	Name	Description	Min Max
	PORO	Reservoir Porosity	0.15 0.35
	KSAND	Reservoir Permeability	10	300
	KRSAND Water relative permeability end-point 0.5	1.0

Table 1 :

 1 Uncertain Parameters

Table 2 :

 2 CPU time comparison

Acknowledgements

This work has been partially supported by the French National Research Agency (ANR) through COSINUS program (project COSTA-BRAVA n ANR-09-COSI-015). I also would like to thanks Nicolas Mourand and Dan BOSSIE CODREANU for the presented CO 2 reservoir model.