Numerical Method for Control of Piecewise-Deterministic Markov Processes

The aim of this chapter is to present the computational method developed in [DES 10] for the control of a Piecewise-Deterministic Markov Process (PDMP) (X t ) t 0 and to explain how similar ideas can be used for statistical inference. Most technical details are omitted to focus on the practical application of the procedure: we state in details the class of PDMPs it can be applied to and give all the algorithms necessary to its implementation. The main material of this chapter was originally published as [DES 10] and [BRA 12b]. Additional details and examples can also be found in [DES 15].

Introduction

Roughly speaking, a control problem is an optimization problem where a controller acts on a process, either continuously on the flow or jump rate or punctually at jump times on the jump kernel Q or even by creating new jumps and selecting new post-jump locations, in order to maximize a gain or minimize a cost. In this chapter, we deal with the most simple form of control, namely optimal stopping. The controller only chooses an admissible time to stop the process in order to maximize a reward and the dynamics of the process before stopping is not controlled.

We are interested in the optimal stopping problem with random horizon T N , where N 0 is a nonnegative integer and T N the N -th jumps time of the PDMP. The family M N of admissible times are the stopping times with respect to the natural filtration (F t ) t 0 of the PDMP (X t ) t 0 that are bounded by T N . The reward is the expectation Chapter written by Benoîte DE SAPORTA and François DUFOUR.

of some function g of the process at the chosen stopping time. The value function v is the best possible performance. It depends on the deterministic starting point x 0 of the process

v(x 0 ) = sup ⌧ 2MN E x0 [g(X ⌧ )]. [6.1]
An optimal stopping time is a stopping time ⌧ ⇤ 2 M N reaching this supremum:

v(x 0 ) = E x0 [g(X ⌧ ⇤ )].
In general, optimal stopping times do not exist and the next best thing is a so-called ✏-optimal stopping time ⌧

⇤ ✏ 2 M N satisfying v(x 0 ) ✏  E x0 [g(X ⌧ ⇤ ✏ )]  v(x 0 ).
Our objective is to propose a numerically tractable approximation of the value function v, and to construct a simulatable ✏-optimal stopping time.

Optimal stopping problems for PDMPs have been extensively studied in the literature from the theoretical point of view, see e.g. [COS 88, COS 00, DAV 93, G ' AT 91, GUG 86, LEN 85]. In [GUG 86] the author defines a dynamic programming operator related to the first jump time of the process, and shows that the value function of the optimal stopping problem is a fixed point for this operator. The basic assumption in this case is that the final cost function is continuous along trajectories, and it is shown that the value function will also have the same property. In [G ' AT 91, LEN 85] the authors adopt some stronger continuity assumptions and boundary conditions to show that the value function of the optimal stopping problem satisfies some variational inequalities, related to integro-differential equations. In [DAV 93], M.H.A. Davis assumes that the value function is bounded and locally Lipschitz along trajectories to show that the variational inequalities are necessary and sufficient to characterize the value function of the optimal stopping problem. In [COS 00], the authors weakened the continuity assumptions of [DAV 93, G ' AT 91, LEN 85]. To the best of our knowledge, the only paper addressing computational issues for optimal stopping of PDMPs is [COS 88]. The numerical approach uses a fixed discretization of the state space based on the kernel Q. The authors derive a convergence result for the approximation scheme but no estimation of the rate of convergence is provided. They also provide ✏-optimal stopping times based on the level sets of the approximate value function, but they require a continuous-time minimization and are not fully numerically tractable, except for simple special cases. Although the literature on numerical methods for PDMPs is still quite scarce, that for diffusion processes is especially rich. We focused our attention on the optimal quantization approach developed for numerical probability in the series of papers [PAG 98, BAL 03, PAG 04b, PAG 04a, BAL 05]. Optimal quantization consists in replacing a continuous state-space random variable by the best discrete one in the sense that the L 2 -norm of the difference is minimum. The difference is given by Zador theorem [ZAD 82] and goes to zero as the number of points K in the support of the discrete distribution goes to infinity, at a speed proportional to K 1/d where d is the dimension of the original random variable. In reasonably small dimension, this result is especially appealing as it can lead to explicit error bounds for approximation schemes based on quantization. In addition, there exist simulation-based algorithms providing the discrete quantized variable.

Our approach is the following. We start form the dynamic programming operator defined by [GUG 86], and rewrite it as an expectation of a functional of the PDMP. Then, we use the key property that all the randomness of the PDMP is contained in its pots-jump locations and inter-jump times Markov chain (Z n , S n ) n 0 , see Subsection 1.2.1 in Chapter 1. The dynamic programming operator can thus be rewritten as a conditional expectation involving this Markov chain. The next step is to replace the continuous state space chain (Z n , S n ) n 0 by its quantized approximation, thus turning conditional expectations into finite weighted sums. The discretized operator can now be iterated numerically, and we estimate the difference between the original value function and this approximation. A by-product of the discretization is a numerically tractable stopping rule that is an ✏-optimal stopping time for the original PDMP.

The main idea of this approach is first to rewrite the quantity to be computed as an expression of the embedded Markov chain (Z n , S n ) n 0 and second to replace it by its quantized approximation. It is our belief that many problems related to PDMPs, including statistical estimation problems, can be numerically approximated this way. In particular, this idea was successfully used in other control problems, such as optimal stopping problems under partial observation [BRA 13], impulse control problems [DES 12a] or more recently change-point detection problems [CLE 17], but also to approximate expectations of functionals of PDMPs [BRA 12a] or the distribution of exit times [BRA 12b]. We detail the latter problem here to emphasize the versatility of the approach that we propose. This chapter is organized as follows. As our approach is simulation based, we start by recalling how to simulate a PDMP in Section 6.2. In Section 6.3, we build the numerical approximation of our optimal stopping problem, as well as the proposed ✏-optimal stopping time. In Section 6.4, we explain how similar ideas can be used to estimate the distribution of an exit time. Finally, in Section 6.5, both problems are illustrated on an industrial example of corrosion provided by Astrium.

Simulation of Piecewise-Deterministic Markov Processes

We will assume in the sequel that the flow can be computed analytically, otherwise, one can use any suitably precise differential equation solver. Still, it is worth noting that very few results in the literature study the additional error caused by the use of a numerical differential equation solver into the global simulation. The only source of randomness in a PDMP comes from the jumps. In order to simulate a PDMP, one must thus be able to simulate its jump times and its post-jump locations. As in most examples simulating post-jump locations is straightforward, we concentrate on simulating the inter jump-times, which boils down to simulating a random variable with time-dependent intensity, as seen from Eq. [1.2] in Chapter 1.

As always, if the cumulative distribution function of the inter-jump time can be inverted analytically, it is preferable to simulate the inter-jump time by direct inversion. This can be done for instance when the intensity is fixed (exponential distribution), or for the Weibull distribution.

If the cumulative distribution function cannot be inverted analytically, one can use the thinning algorithm 6.1 based on the Poisson distribution that simulates a random variable with time dependent intensity. It is especially suitable when one just wants to simulate the embedded Markov chain (Z n , S n ) n 0 instead of the whole continuoustime trajectory of the PDMP. It only requires a local maximum of the intensity. If the flow must be solved numerically, Gillespie-type algorithms are usually more efficient, see e.g. [LEM 17].

Optimal Stopping

We start with stating the main regularity assumptions that we require in Subsection 6.3.1, then in Subsection 6.3.2 we give the exact form of the dynamic programming operators. In Subsection 6.3.3 we introduce the quantized approximation of the value function as well as an ✏-optimal stopping time.

Assumptions and Notation

In order to fully use the properties of optimal quantization, we make the main assumption that the driving parameters of the PDMP as well as the reward function of the optimal stopping problem are bounded and Lipschitz continuous.

Let L be the subset of C b of Lipschitz-continuous functions, and for any function w 2 L, denote by C w its upper bound and by L w its Lipschitz constant

C w = sup x2X |w(x)|, L w = sup x6 =x 0 2X |w(x) w(x 0 )| |x x 0 | .
Denote by L the set of functions that are Lipschitz continuous along the flow i.e. the real-valued, bounded, measurable functions w defined on X and satisfying the following conditions: -for all x 2 X, the map w( (x, •)) : [0, t + (x)) ! R is continuous, the limit lim t!t + (x) w( (x, t)) exists and is denoted by w (x, t + (x)) , -there exists [w] 1 2 R + such that for x, y 2 X and t 2 [0, t + (x) ^t+ (y)], one has

|w( (x, t)) w( (y, t))|  [w] 1 |x y|; -there exists [w] 2 2 R + such that for all x 2 X and t, t 0 2 [0, t + (x)], one has |w( (x, t)) w( (x, t 0 ))|  [w] 2 |t t 0 |; -there exists [w] ⇤ 2 R + such that for all x, y 2 E, one has |w( (x, t + (x))) w( (y, t + (y)))|  [w] E ⇤ |x y|.
Note that any function w in L is also in L with L w  [w] 1 as for all x 2 X, x = (x, 0). Conversely, if the flow and the exit time t + are bounded and Lipschitz continuous:

2 L(X ⇥ R + ) and t + 2 L then any function w 2 L is in L with [w] 1  L w L , [w] 2  L w L and [w] ⇤  L w L (L t + + 1).
We can now state the main assumptions that will be in force throughout this chapter. ASSUMPTION 6.1.-Assume that:

-the deterministic time to reach the boundary of the state space is bounded and Lipschitz-continuous: t + is in L;

-the jump intensity is bounded by C and there exists [ ] 1 2 R + such that for all x, y 2 X and t 2 [0, t + (x) ^t+ (y)], one has

| ( (x, t)) ( (y, t))|  [ ] 1 |x y|; -the Markov kernel Q is Lipschitz in the following sense: there exists [Q] 2 R + such that for all functions w 2 L , 1) for all x, y 2 X and t 2 [0, t + (x) ^t+ (y)), one has |Qw( (x, t)) Qw( (y, t))|  [Q][w] 1 |x y|,
2) for all x, y 2 X , one has

|Qw( (x, t + (x))) Qw( (y, t + (y)))|  [Q][w] ⇤ |x y|; -the reward function g is in L .
In most practical applications, the physical properties of the system ensure that either t + is bounded, or the problem has a natural deterministic time horizon T . In the latter case, there is no loss of generality in considering that t + is bounded by this deterministic time horizon. The other assumptions are technical and easy to check on most applications.

We finish this section with some generic notation. For a, b 2 R, a _ b = max(a, b) is maximum of a and b. For a Markov kernel P on (X , B(X )) and functions w and w 0 in C b , set

P w(x) = Z X w(y)P (x, dy), (wP w 0 )(x) = w(x) Z X w 0 (y)P (x, dy),
for any x 2 X .

Dynamic Programming

We now recall the dynamic programming equations derived in [GUG 86]. They were originally given in an analytical form depending on the local characteristics , and Q of the PDMP, but we give their form as expectations as it is the suitable form for our intended discretization. For all x 2 X , w 2 L and t 0, one defines

Hw(x, t) = w (x, t ^t+ (x)) P S n+1 t ^t+ (Z n )|Z n = x , Iw(x, t) = E ⇥ w(Z n+1 )1 {Sn+1<t^t + (Zn)} |Z n = x ⇤ , Jw(x, t) = Hg(x, t) + Iw(x, t), Kw(x) = E ⇥ w(Z n+1 )|Z n = x ⇤ , Lw(x) = sup ut + (x) Jw(x, u) _ Kw(x).
As (Z n , S n ) n 0 is a time-homogeneous Markov chain, the quantities above do not depend on n. Roughly speaking, operator L selects the best compromise between waiting for the next jump (operator K) or stopping at the best position along the flow (sup J). The value function v can be recursively constructed by iterating operator L as follows. Set, for 0  n  N ,

v N = g, v n = L(v n+1 , g).
The last term v 0 of this recurrence is then exactly the value function v. In order to obtain a numerical approximation of v, it is thus sufficient to be able to discretize and iterate operator L. This is not straightforward as L involves a continuous supremum as well as conditional expectations. Our first step is to replace the recursion on functions v n by a recursion on random variables. Set

V N = v n (Z n ) = g(Z n ) and, for 0  n  N , V n = v n (Z n ) = sup ut + (Zn) n E ⇥ V n+1 1 {Sn+1<u} + g (Z n , u) 1 {Sn+1 u} Z n ⇤ o _ E ⇥ V n+1 Z n ⇤ .
It is a highly remarkable property of this problem that the dynamic programming equation involving a sequence of functions (v n ) can simply be rewritten as a selfcontained dynamic programming equation involving only the random variables (V n ).

As random variables are much more tractable numerically than functions, this is a crucial step towards a numerical approximation of the value function.

Quantized Approximation

The next step is to replace the Markov chain (Z n , S n ) n 0 by its optimal quadratic quantization. This can be done for instance by using the Competitive Learning Vector Quantization (CLVQ) algorithm described in Algorithm 6.2, see e.g. [BAL 03, Section 3.3] for details. Let ( n ) 0nN be the optimal quantization grids for (Z n , S n ) 0nN , and p n the nearest-neighbor projection from X ⇥ onto n . One thus defines the quantized approximation ( Ẑn , Ŝn ) 0nN of (Z n , S n ) 0nN by

( b Z n , b S n ) = p n (Z n , S n ).
We also denote by Z n the projection of n onto X .

For z 2 X , set (z) 2 (0, t + (z)). Define n(z) = ⌅ t + (z) (z)
⇧ 1, where bxc denotes the greatest integer smaller than or equal to x. The set of points (t i ) 0in(z)} with t i = i (z) is denoted by G(z). This is the grid associated to the time interval [0, t + (z)].

The process ( Ẑn , Ŝn ) 0nN is not a Markov chain and especially not time-homogeneous. Therefore we define a sequence of approximate operators. For

n 2 {1, . . . , N}, w 2 L( Z n ), z 2 Z n 1 , and s 2 R + , set b J n (w, g)(z, s) = E ⇥ w( b Z n )1 {b Sn<s} + g (z, s) 1 {b Sn s} b Z n 1 = z ⇤ , b K n (w)(z) = E ⇥ w( b Z n ) b Z n 1 = z ⇤ , b L d n (w, g)(z) = max s2G(z) b J n (w, g)(z, s) _ b K n (w)(z).
From these operators, we define a new sequence of random variables b

V n = b v n ( b Z n ) where b v N (z) = g(z) for z 2 Z N and b v n 1 (z) = b L d n (b v n , g)(z), [6.2]
for 1  n  N and z 2 Z n 1 . Note that the random variables b V n are discrete and can be recursively computed as finite weighted sums. The last term of the recurrence b V 0 is our chosen approximation for the value function v at the starting point x 0 of the PDMP. It is detailed in Algorithm 6.3. It also contains additional information related to which operator won the maximization at each step that will be useful to construct an ✏-optimal stopping time. We can also obtain an upper bound for the error in terms of the regularity constants and the quantization error.

for 0  n  N do for 1  i  K do (W n ) i w i,n /N R end end for 1  n  N do for 1  i  K do for 1  j  K do (P n ) i,j p i,j,n /N R end end end return: ( NR n ), W n , 0  n  N , P n , 1  i  N . end Algorithm 6.3: Approximation of the value function v(x 0 ) input : Quantization grids ( ⇥ n ) 0nN , Transition matrices (P n ) 1nN , cost function g, time discretization step output : b v n , indices for ✏-optimal stopping time: s ⇤ n , i n , 1  n  N 1 begin 2 for z 2 Z N do 3 b v N (z) g(z) 4 end 5 for n N 1 to 0 do 6 for z 2 Z n do 7 b Kb v n+1 (z) X (z 0 ,s 0 )2 n+1 P n+1 (z, (z 0 , s 0 ))v n+1 (z 0 ) 8 for u 0 to t ⇤ (z) by step: 9 do 10 b J b v n+1 (z, u) X (z 0 ,s 0 )2 n+1 P n+1 (z, (z 0 , s 0 )) v n+1 (z 0 )1 {s 0 <u} + g( (z, u))1 {s 0 u} 11 end 12 b v n (z) max u { b J b v n+1 (z, u)} _ b Kb v n+1 (z) 13 s ⇤ n+1 (z) min u ⇤ ; b J b v n+1 (z, u ⇤ ) = max u { b J b v n+1 (z, u)} 14 if b Kb v n+1 (z) > max u { b J b v n+1 (z,
More importantly for applications, from this procedure one can also construct an ✏optimal stopping time, see Algorithm 6.4. Its rigorous construction and the proof that it is a stopping time for the original PDMP and ✏-optimal can be found in [DES 10]. An example of application of this procedure is given in Section 6.

7 if i n+1 (z) = 1 then 8 r t ⇤ (Z n ) 9 else 10 r s n+1 (z)1 {sn+1(z)<t ⇤ (Zn)} + (t ⇤ (Z n ) )1 {sn+1(z) t ⇤ (Zn)}

Exit Time

Our approach to approximate the optimal stopping problem can be summarized as follows. First, we obtain a recursive construction of the quantity of interest, namely the dynamic programming equations for the value function. Hopefully, this can generically be done using the Markov property. Second, we rewrite the recursion operator in terms of the embedded Markov chain (Z n , S n ) n 0 . Again, this can be generically done as this chain is the only source of randomness in a PDMP. Third, we replace the continuous state space chain (Z n , S n ) n 0 by its quantized approximation ( b Z n , b S n ) n 0 leading to finite weighted sums. This way, one obtains a fully computable approximation, and one can hope to obtain error bounds and even convergence rates providing the underlying functional are Lipschitz-continuous. It is our belief that this program can be successfully applied to other problems related to PDMPs, and in particular to statistical estimation problems.

To illustrate this potential, we present in this section the exit time distribution estimation procedure detailed in [BRA 12b] or [DES 15, Chapter 5]. Our aim is not to provide all the technical details, but only to detail the class of PDMPs it can be applied to and the algorithms. We define the exit time problem in Subsection 6.4.1, then we give the recursive formulation of the problem in Subsection 6.4.2.

Problem Setting and Assumptions

Let (X t ) t 0 = (Y t , M t ) t 0 be a PDMP on the state space X = {(y, m), m 2 M, y 2 Y m } with mode set M. For all mode m 2 M, let U m be an open subset of Y m and set U = {(y, m), m 2 M, y 2 U m }. We deal with the exit time ⌧ of (X t ) t 0 from U ⌧ = inf{s 0 : X s 6 2 U }.

Denote by µ the distribution of the initial state of the process X 0 = Z 0 . To make our exit problem relevant, we assume that the process starts in U and eventually leaves it almost surely i.e. the support of µ is included in U and P µ (⌧ < +1) = 1. The aim of this chapter is to provide an approximation scheme for the survival function of ⌧ .

As our approach is based on a recursive computation using the underlying discretetime Markov chain of the PDMP, we actually use a computation horizon that equals a jump time of the process: we will study ⌧ ^TN rather than ⌧ itself. We assume that the PDMP is non explosive, i.e. lim n!1 T n = +1, so that as N goes to infinity, one gets

⌧ ^TN ! ⌧.
We assume that N is chosen such that P µ (⌧ > T N ) is small enough.

Let us define u ⇤ (x) for all x 2 U as the time for the flow starting from the point x to exit from U u ⇤ (x) = inf{s 0 : (x, s) 6 2 U }, where u ⇤ will play a similar role as t + in the previous section. We make the following assumptions. ASSUMPTION 6.2.-Assume that:

-the function u ⇤ is in L with upper bound C u ⇤ and Lipschitz constant L u ⇤ ; -for all m 2 M, the set U m is convex;

-for ↵ > 0, set U ↵ = x 2 X : min{|x y| : y 2 @U }  ↵ .
There exist C U > 0 and > 0 such that for all 0

 n  N , P µ (Z n 2 U ↵ )  C U ↵ .
The last two assumptions are technical and ensure that the difference between Z n and its quantized approximation b Z n can be controlled even though they are close to the boundary of U . Finally, we assume that once the process has exited U , it cannot get back in U . If it is not the case naturally, one may just kill the process as it exits U .

Recursive Formulation

On can obtain a recursive construction of the survival function of the truncated exit time as follows. For all s > 0, define the sequences (p n (s)) n2N , (q n ) n2N and (r n (s)) n2N by

p n (s) = P µ (⌧ > s ⌧  T n ), q n = P µ (⌧  T n ), r n (s) = P µ ({⌧ > s} \ {T n < ⌧  T n+1 }).
The conditional probability p n (s) does not exist when q n = 0. We then choose to extend the sequence by setting p n (s) = 0 in this case. Our objective is to approximate p N (s). The following proposition provides a recursion for the sequence (p n ) nN , pointing out that p N may be computed as soon as the sequences (q n ) nN and (r n ) nN 1 are known. PROPOSITION 6.1.-For all n 2 N, s > 0, one has p 0 (s) = 0 and

p n+1 (s) = p n (s)q n + r n (s) q n+1 1 {qn+16 =0} .
Algorithm 6.5: Approximation of (q n ) and (r n (s))

input : Quantization grids n , 0  n  N , Weight vectors W n , 0  n  N , Transition matrices P n , 1  n  N output: (b q n ) and (b r n (s)) 1 begin 2 b q 0 0 3 for n 0 to N 1 do 4 b q n+1 X (z,t)2 n+1 z / 2U W n+1 (z, t) 5 b r n (s) X (z,t)2 n z2U X (z 0 ,t 0 )2 n+1 z 0 / 2U 1 {(t+u ⇤ (z))^t 0 >s} W n (z, t)P n ((z, t), (z 0 , t 0 )) 6 end return: (b q n ), (b r n (s)) 7 end
It is important to note that both sequences (b q n ) and (b r n (s)) may be computed easily from the output of the quantization algorithm as shown in Algorithm 6.5. The approximation of p n is then readily obtained by plug-in. For all s > 0 and for all n  N 1, let b p 0 (s) = 0 and

b p n+1 (s) = b p n (s)b q n + b r n (s) b q n+1 1 {b qn+16 =0} .
We can prove that this approximation scheme converges, see [BRA 12b] or [DES 15, Chapter 5] for the proof. THEOREM 6.2.-For all n  N and for almost every s > 0, one has

b p n (s) ! p n (s),
when the quantization error goes to zero.

Numerical Example

We study the evolution of the thickness of an aluminum metallic structure subject to corrosion. This model was provided by Astrium and is fully described in [DES 12b]. It concerns a small homogeneous structure within a strategic ballistic missile. The missile is stored successively in three different environments, the workshop, the submarine in operation and the submarine in dry-dock. Then it goes back to the workshop and so on. The missile stays in each environment during a random duration with exponential distribution. Its parameter depends on the environment. At the beginning of its service time, the structure is treated against corrosion. The period of effectiveness of this protection is also random, with a Weibull distribution. The thickness loss only begins when this initial protection is gone. The degradation law for the thickness loss then depends on the environment through two parameters, a deterministic transition period and a random corrosion rate uniformly distributed within a given range. Typically, the workshop and dry-dock are the more corrosive environments. The randomness of the corrosion rate accounts for small variations and uncertainties in the corrosiveness of each environment.

Piecewise-Deterministic Markov Model

The dynamics of the thickness loss for the structure can be described by a PDMP. The finite set of modes is

M = {1, 2, 3} ⇥ {0, 1},
where the first coordinate corresponds to the environment: 1 for the workshop, 2 for the submarine in operation and 3 for the dry-dock, and the second coordinate states whether the corrosion has started (1) or not (0). The state space in modes m 2

{1, 2, 3} ⇥ {0} is Y m = (0, +1) ⇥ R + ,
with boundary @Y m = {0} ⇥ R + , and the Euclidean variables are Y t = ( t , a t ) the remaining time of effectivity of the protection, and the time since the last change of mode. The state space in modes m = (i, 1) 2 {1, 2, 3} ⇥ {1} is

Y m = R + ⇥ [⇢ i , ⇢ + i ] ⇥ R + ,
with empty boundary and the Euclidean variables are Y t = (d t , ⇢ t , a t ) the thickness loss, the corrosion rate of the given environment, and time since the last change of mode.

Originally at time 0, the mode is M 0 = (1, 0) and Y 0 = ( 0 , 0), where 0 is drawn according to a Weibull distribution with parameters ↵ = 2.5 and = 11800 hours 1 . In other words, the process starts in the workshop environment with a random protection that will be effective during 0 hours. The flow in modes

m 2 {1, 2, 3} ⇥ {0} is m (( , a), s) = ( s, a + s),
with t + (m, , a) = , meaning that the boundary is reached when the the protection is gone. In modes m = (i, 1) 2 {1, 2, 3} ⇥ {1}, the flow is

m ((d, ⇢, a), s) = (d + d i (⇢, a + s) d i (⇢, a), ⇢, a + s),
where

d i (⇢, s) = ⇢ s + ⌘ i (e s/⌘i 1) ,
is the dynamics of the thickness loss in environment i. The parameters are given in Table 6.1. All the parameters come from expert opinion. The boundary is empty in this case, so that t + (m, d, ⇢, a) = +1. We will discuss in the next section how to change the value of t + so that our boundedness assumption may hold true.

The jump intensity in environment i is a constant i , which value is given in Table 6.1. The jump kernel is as follows:

-if m = (i, 0), y = ( , a) with > 0, Q((y, m), {( , 0)} ⇥ {m 0 }) = 1 {(0,m 0 =(i+1[3])} , with i + 1[3] = i + 1 if i  2 and 3 + 1[3] = 1; -if m = (i, 0), y = (0, a), and A is in B(R + ) Q((y, m), {0} ⇥ A ⇥ {0} ⇥ {m 0 }) = 1 {m 0 =(i,1)} 1 ⇢ + i ⇢ i Z ⇢ + i ⇢ i 1 A (s)ds,
where the parameters ⇢ ± i are given in Table 6.1; -if m = (i, 1), y = (d, ⇢, a), and A is in

B(R + ) Q((y, m), {d} ⇥ A ⇥ {0} ⇥ {(j, 1)}) = 1 {j=i+1[3]} 1 ⇢ + j ⇢ j Z ⇢ + j ⇢ j 1 A (s)ds.
Figure 6.1 shows some simulated trajectories of the thickness loss along time. The slope changes correspond to changes of environment. The observed dispersion is characteristic of the random nature of the phenomenon, and illustrates the poorness of an average rule to monitor such a process. One needs to make decisions adapted to each specific trajectory. The structure is considered unusable if the loss of thickness reaches 0.2 mm.

Deterministic Time to Reach the Boundary

Because the deterministic exit time is unbounded in modes m 2 {1, 2, 3} ⇥ {1}, the PDMP model described above does not fit our framework. There are two main ways to overcome this difficulty.

The first one consists in killing the process when the thickness lost by corrosion reaches the prescribed level of 0.2 mm. The state space in modes m = (i, 1) 2 {1, 2, 3} ⇥ {1} thus becomes

Y m = [0, 0.2) ⇥ [⇢ i , ⇢ + i ] ⇥ R + , with boundary @Y m = {0.2} ⇥ [⇢ i , ⇢ + i ] ⇥ R + ,
and deterministic time to boundary is

t + ((i, 1), d, ⇢, a) = inf{t > 0 : d + d i (⇢, a + t) d i (⇢, a) = 0.2}.
Thus, t + is uniformly bounded on Y m for m 2 {1, 2, 3} ⇥ {1}. Rigorously, this does not hold true for all points in Y m with m 2 {1, 2, 3} ⇥ {0} as the support of the Weibull law is unbounded. However, the initial protection is drawn only once and for a given initial protection time , the quantities t + (X t ) are bounded uniformly for all t, which is actually all that is required in our proofs.

The second way to bound the deterministic time to the boundary is by considering the time-augmented process ( e X t ) t 0 = (X t , t) t 0 with a finite fixed horizon T for the trajectories. Thus, the state space becomes

Y m = (0, +1) ⇥ R + ⇥ [0, T ), in modes m 2 {1, 2, 3} ⇥ {0}, with boundary @Y m = {0} ⇥ R + ⇥ [0, T ) [ (0, +1) ⇥ R + ⇥ {T } [ {0} ⇥ R + ⇥ {T },
and deterministic time to reach the boundary

t + (m, , a, t) = ^(T t).
The state space in modes m = (i, 1) 2 {1, 2, 3} ⇥ {1} becomes

Y m = R + ⇥ [⇢ i , ⇢ + i ] ⇥ [0, T ) ⇥ [0, T ),
with boundary

@Y m = R + ⇥ [⇢ i , ⇢ + i ] ⇥ [0, T ) ⇥ {T },
and deterministic time to reach the boundary

t + (m, d, ⇢, a, t) = (T t),
as a  t by construction, the boundary may not be reached by a. Thus, t + is bounded by T uniformly in all the modes and state variables.

Quantization

In this example, the main difficulty to implement the CLVQ quantization algorithm lies in the fact that the variables of interest have very different scales: from about 10 6 for the corrosion rate ⇢ to 10 5 for the average time spent in environment 2. This poses a problem as searching for the nearest neighbor and gradient calculations are done in Euclidean norm, regardless of the magnitudes of the components. This defect is corrected by a renormalization of the variables. Therefore a weighted Euclidean norm is used to quantify the Markov chain associated with the corrosion process to rescale the variables. Projections of the quantization grids obtained for 2000 points on the variables inter-arrival times S n (abscissa) and thickness loss at the jump times d Tn (ordinate) are given in Figure 6.2. Notice how different the grids are from one jump time to the other. Note also that the Markov chain does not seem to follow any well-known probability distribution, and that if we had used a cartesian grid on the rectangle, a huge number of points would have been useless. Thus, the quantization approach is especially suitable as it allows for dynamically changing grids along time instead of a single fixed grid.

Optimal Stopping

The reward function g is given in Figure 6.3. It depends only on the thickness loss, and reflects the fact that on the one hand, one should not perform a maintenance too early, when the thickness loss is low (low gain) but on the other hand one should not cross the security threshold of 0.2 mm (zero gain). It also satisfies our Lipschitz continuity requirements. The reward function g is maximum (value 4) at abscissa 0.18 and the thickness loss is continuous and non-decreasing. Thus, it is easy to see that the true value function at our starting point is 4, which is the maximum of the reward function g, and an optimal stopping time is the first moment when the loss reaches 0.18 mm thick.

Our numerical procedure is valid for any sufficiently regular reward function, and the knowledge of the true value function or optimal stopping time shall not be used in the numerical procedure. Besides, recall that the thickness loss is not measured continuously. However, this exact optimal stopping time will serve us as a benchmark to evaluate the performance of our approximation. Figure 6.4 shows two examples of computation of the quasi optimal maintenance time on two specific simulated trajectories. The thick vertical line represents the moment provided by the algorithm to perform maintenance. The other vertical lines materialize the moments of change of environment, the horizontal dotted line the theoretical optimum. In both examples, we stop at a value very close to the optimum value. In addition, the intervention did take place before the critical threshold of 0.2 mm. An approximate value function v was calculated in two ways. The first one is the direct method obtained by the algorithm described above. The second one is obtained by Monte Carlo simulation using the quasi-optimal stopping time provided by our procedure. The numerical results obtained are summarized in Table 6.2. As expected, the greater the number of points in the quantization grid, the better our approximation becomes. Furthermore, the specific form of the function g indicates that at the threshold of 1, the intervention takes place between 0.15 and 0.21 mm and when the threshold increases, this range is narrowed. Our approximation is therefore good even for low numbers of grid points. The last column of the table also shows the validity of our stopping rule. It should be noted here that this rule does not use the optimal stopping time stop at the first moment when the thickness loss reaches 0.18 mm. The method that we propose is general, and implementable even when the optimal stopping time is unknown or does not exist, as will be illustrated in the following sections.

Moreover, one can also construct a histogram of the values of our stopping time, that is to say, a histogram of the values of effective moments of maintenance and compare it with the exact optimal stopping time. Figure 6.5 shows the results for 10 5 Monte Carlo simulations. Image (a) is obtained by our approximation procedure and

Conclusion

We have presented a general framework for obtaining numerical approximations of quantities related to PDMPs. First, one obtains a recursive formulation of the quantity to be computed, thanks to the Markov property. Second, one rewrites the expression using the embedded Markov chain (Z n , S n ). Third, one replaces this Markov chain by its quantized approximation in order to obtain fully computable quantities. This approach has been successfully conducted on an optimal stopping problem as well as on an exit time estimation problem.

Computing the quantization grids may be very time consuming, and will lead to poor results in high dimension. However, in small enough dimension, this computation can be done off-line in advance. The same grids can also be reused to compute different quantities. Thankfully, they require only a simulator of the process to be computed. On could also imagine to use observation data to construct such grids, as in [CHA 15], if enough data is available. This approach can also be used in other statistical inference problems.

Bibliography

[BAL 03] BALLY V., PAGÈS G., "A quantization algorithm for solving multi-dimensional discrete-time optimal stopping problems", Bernoulli, vol. 9, num. 6, p. 1003-1049, 2003.

Algorithm 6. 1 :

 1 Thinning algorithm to simulate a time dependent intensity input : Intensity function (t), positive real number A output: One realization of the random variable with intensity (t) = ; do 5 select L sup ata+A (t) 6 draw n from a Poisson distribution of parameter L ⇥ A 7 if n 6 = 0 then 8 for 1  k  n do 9 draw u k from a uniform distribution on [a; a + A] 10 draw v k from a uniform distribution on [0; L] 11 if v k  (u k ) then 12 S = S [ {u k }

, 1

 1 Under our regularity assumptions, we have the following convergence result. Its proof can be found in [DES 10] or [DES 15, Chapter 7]. THEOREM 6.1.-The difference k V0 V 0 k goes to zero as the number of points in the quantization and time discretization grids go to infinity. Algorithm 6.2: CLVQ algorithm to quantize a Markov chain input : Number of points K in each grid, Number of runs NR, Horizon N , Initial grids ( 0 n ) 0nN with K points each, Sequence ( n ) Simulator of trajectories of target Markov chain output: Optimized grids ( ⇤ n ) and weights (W n ) 0  n  N , Transition probabilities P n , 1  n  N begin w 0 (size K ⇥ N + 1) p 0 (size K ⇥ K ⇥ N ) for m 0 to NR 1 do simulate trajectory (x 0 , x 1 , . . . , x N ) according to law of Markov chain for n 0 to N do competitive phase: select i n in {1, . . . , K} such that y m in is the closest neighbor of x n in m n  i  K} updating weights: w in,n w in,n + 1 if n > 0 then p in 1 ,in,n p in 1 ,in,n + 1

Figure 6 . 1 :

 61 Figure 6.1: Simulated trajectories of the thickness loss along time for the corrosion process.

Figure 6 . 2 :

 62 Figure 6.2: Quantization grids (2000 points) for the corrosion process: inter-arrival times (abscissa) and thickness loss at the jump times (ordinate).

Figure 6 . 3 :

 63 Figure 6.3: Reward function for the maintenance of the corrosion process.

  Figure 6.4: Examples of stopped trajectories for the corrosion process.

Figure 6 . 6 :

 66 Figure 6.6: Survival function of the service time of the corrosion process obtained through Monte Carlo simulations (dashed), quantized approximation (solid) and the error with 2000 points in the quantization grids.

  5, other examples can be found in [DES 15, Chapter 9]. Computation of an ✏-optimal stopping time input : Quantization grids ( n ) 0nN , Outputs of Algorithm 6.3: s ⇤ n , i n , 1  n  N , positive real number , simulator of trajectories of the PDMP, projection operators p n onto the quantization grids

	Algorithm 6.4: output: Stopping date ⌧
	1 begin
	2	n	0
	3	⌧	NaN
	4	while n < N do

5 simulate (Z n , S n , S n+1 ) 6 (z, s) p n (Z n , S n )

Table 6 .

 6 1: Numerical values of the parameters of the corrosion model.

							environment 1 environment 2 environment 3
				i (h 1 )			(17520) 1 (131400) 1	(8760) 1
			⌘ i (h)						30000	200000	40000
			⇢ i (mm•h 1 ) ⇢ + i (mm•h 1 )				10 6 10 5	10 7 10 6	10 6 10 5
	0.25
	0.2
	0.15
	0.1
	0.05
	0
	0	1	2	3	4	5	6	7	8	9
										x 10 5

Table 6 .

 6 2: Numerical results for the calculation of the value function for the corrosion process. Number of points Approximation of the Approximation of the value in the quantization value function by the function by Monte Carlo with the

	grids	direct algorithm	quasi-optimal stopping time
	10	2.48	0.94
	50	2.70	1.84
	100	2.94	2.10
	200	3.09	2.63
	500	3.39	3.15
	1000	3.56	3.43
	2000	3.70	3.60
	5000	3.82	3.73
	8000	3.86	3.75

Proof of Proposition 6.1.-First, recall that T 0 = 0 so that one has p 0 = 0 since the process starts in U . Then, let n 2 N such that q n+1 6 = 0 and note that {⌧  T n+1 } = {⌧  T n } [ {T n < ⌧  T n+1 }. Thus, one has

showing the result. ⇤

It remains to express the above quantities in term of the embedded Markov chain. We work with the post-jump locations and jump times chain (Z n , T n ) n2N rather than the post-jump locations and inter-jump times chain (Z n , S n ) n2N . They play similar roles as S n = T n T n 1 for all n 1 and our recursions already involve (T n ), therefore it is natural to keep using this sequence.

First, note that

as the process cannot go back to U once it exited U . Moreover, on the event (Z n 2 U, Z n+1 6 2 U ), one has

as u ⇤ is the deterministic exit time from U . Thus, one obtains

Numerical Approximation

We now naturally define the quantized approximations of the previous sequences. For all s > 0, define

Exit Time

Considering the approximation scheme for the exit time distribution, one may note that the quantized value b p N (s) is not necessary smaller than 1. Therefore, it appears natural to replace b p N (s) by b p N (s) ^1. This does not change the convergence theorem and can only improve the approximation error. Figure 6.6 presents the survival function of ⌧ obtained through Monte Carlo simulations (the dashed line), through our approximation scheme (the solid line) and the error. Table 6.3 contains the empirical error for different numbers of points in the quantization grids.


E Bas

A Menafoglio

T Okhulkova

Polynomial Chaos Expansion for an Efficient Uncertainty and Sensitivity Analysis of Complex Numerical Models

In this work we address the problem of performing uncertainty and sensitivity analysis of complex physical systems where classical Monte-Carlo methods are too expensive to be applied due to the high computational complexity. We consider the Polynomial Chaos Expansion (PCE) as an efficient way of computing a response surface for a model of gas injection into an incompressible porous media aiming at assessing the sensitivity indices and the main distributional features of the maximal spread of the gas cloud. The necessity of an uncertainty study for such a model arises in case of CO2 storage risk assessment and is here performed by jointly using a numerical scheme to solve the system of partial differential equation (PDE) governing the model and the PCE method to efficiently simulate the physical system response by a meta-model. The performances of the PCE method and a standard MC approach are compared through an extended simulation study showing that the computational gain of the PCE approach is remarkable without significant loss in the precision of the estimates.

INTRODUCTION

The uncertainty quantification for a complex numerical model is an important goal for which significant challenges are yet to be achieved. In general, the uncertainty might be associated to several elements of the system such as input data measurements, parameter values, model structure, model solution algorithms and even human behavior. Therefore, the uncertainty assessment under a given cost might not be straightforward.

Most of the research efforts have been focused so far on the development of efficient algorithms to be used in different applications, assuming ideal input data measurements with precisely defined computational domains: to this end, for instance, Monte-Carlo approaches can be employed. However, the simulation of a large-scale system turns in computationally very heavy operations: thus the need of making the convergence to the solution faster naturally arises. The classical Monte-Carlo approach is used in such cases as the last resort. Another option is the perturbation method where all the stochastic quantities are expanded around the mean via Taylor series, nevertheless, it is limited to small perturbations and does not provide information on high-order statistics of the response. The expansion of the inverse of the stochastic operator in a Neumann series can be considered, but it is again limited to small fluctuations [START_REF] Shinozuka | Response variability of stochastic finite element systems[END_REF].

A more effective approach pioneered by [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF] in the context of finite element method is the Polynomial Chaos Expansion, which is a method based on a spectral representation of uncertainty. This idea comes from the theory about polynomial chaos by [START_REF] Wiener | The Homogeneous Chaos[END_REF][START_REF] Wiener | Nonlinear Problems in Random Theory[END_REF], which was firstly applied in the study of turbulence in the 1960s [START_REF] Meecham | Wiener-Hermite Expansion in Model Turbulence at Large Reynolds Numbers[END_REF][START_REF] Siegel | Wiener-Hermite Expansion in Model Turbulence in the Late Decay Stage[END_REF][START_REF] Meecham | Use of the Wiener-Hermite expansion for nearly normal turbulence[END_REF]. Nevertheless, its convergence rate was judged slow [START_REF] Orszag | Dynamical Properties of Truncated Wiener-Hermite Expansions[END_REF][START_REF] Crow | Relationship between a Wiener-Hermite expansion and an energy cascade[END_REF][START_REF] Chorin | Gaussian fields and random flow[END_REF]): hence, it did not receive much attention until the work of [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF]. A rich state of art with many applications can be found in [START_REF] Sudret | Stochastic finite elements and reliability: a state-of-the-art report[END_REF], Xiu and Karniadakis 2002, Xiu et al. 2002).

However, the availability of solid theoretical basis as well as efficient numerical codes (OpenTURNS 2012) for performing uncertainty analysis has recently contribute in making possible the sensitivity analysis of very complex systems [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Formaggia | Global sensitivity analysis through polynomial chaos expansion of a basinscale geochemical compaction model[END_REF][START_REF] Oladyshkin | An integrative approach to robust design and probabilistic risk assessment for CO 2 storage in geological formations[END_REF].

The main purpose of this work is to show how Polynomial Chaos Expansion (PCE) can be used for uncertainty and sensitivity analysis. In particular, we study as a specific application a model of gas injection into porous media, presented in Section 2. The theoretical basis of PCE is briefly introduced in Section 3. The validity of this approach is proved by simulation in Section 4: here, we first perform the sensitivity analysis on our model, then we compare the results of the analytical solution from available simplified model and the results of the Monte Carlo simulation.

INJECTION OF GAS INTO POROUS MEDIA

The predictive numerical model engaged in this study represents a two-phase flow of immiscible compressible fluids through incompressible porous medium. This model is developed in INERIS in order to estimate the intensity of physical processes in underground compartments in case of CO2 storage in a saline aquifer. The injected gas displaces water (initially in place) and simultaneous flow of the two fluids takes place in the porous medium. The model accounts for gravity, capillarity and fluid compressibility effects, and is formulated in terms of CO2 and water pressures. Empirical laws by [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] are used to close the system with the relations between capillary pressure, relative permeability and saturation.

The corresponding system of differential equations is solved within the software COMSOL (2012) in 2Daxisymmetric mode. The source (injection well) is concentric to the symmetry axis.

The upper and lower bounds of the reservoir are considered impermeable, therefore a "zero flux" boundary condition is used. On the bound along the injection well the debit of the gas phase is fixed, on the other side the boundary condition fluid pressures are fixed. The main configuration and the variable of interest are shown in Figure 1.

In Figure 1, p w and p nw [P a] are water and gas pressures respectively, C p [P a -1 ] the coupling term, θ w and θ nw the water and gas volume fractions, k rw and k rnw [-] the water and gas relative permeabilities; Q nw [m/s] represents the gas injection rate, K int [m 2 ] the intrinsic permeability, p w0 and p nw0 [P a] the initial water and gas pressures, ρ w and ρ nw [kg/m 3 ] the fluid densities, µ w and µ nw [P a • s] the fluid viscosities and, finally, χ w and χ nw [P a] are the fluid compressibilities.

The numerical solution of the system at time t = 10 years is finally shown in Figure 2. The uncertainty analysis as well as the sensitivity analysis which will be described in Section 4 are referred to this solution.

POLYNOMIAL CHAOS EXPANSION FOR SENSITIVITY ANALYSIS

Polynomial Chaos Expansion

Let y be the scalar output of a model that depends on the n-dimensional random vector x = (x 1 , ..., x n ) collecting n random parameters x i , valued in Ω i ⊆ R, i = 1, ..., n:

y = f (x), x ∈ Ω ⊆ R n , (1)
where Ω denotes the cartesian product

Ω 1 × • • • × Ω n .
Suppose x 1 , ..., x n to be statistically independent and the response y to be a second-order integrable variable: the PCE of y [START_REF] Maître | Spectral methods for uncertainty quantification[END_REF][START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF], Xiu and Karniadakis 2002) is the spectral decomposition of f over a set of polynomials Ψ i (x), which are orthogonal with respect to the probability measure on Γ (i.e. the law of x), namely [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]:

y = ∞ k=0 β k Ψ k (x).
(2)

The PCE was first introduced for Gaussian input parameters involving Hermite polynomials (e.g., [START_REF] Wiener | The Homogeneous Chaos[END_REF]). In the presence of non-Gaussian distributed variables, two strategies can be adopted in order to preserve the optimal exponential convergence rate: a) transform the input parameters to obtain Gaussian random entries b) consider other families of polynomials (Xiu andKarniadakis 2002, Soize and[START_REF] Soize | Physical systems with random uncertainties: Chaos representations with arbitrary probability measure[END_REF]. In this work, for every input parameter x i , i = 1, ..., n, a log-normal distribution ρ Γ i (x i ) is considered. In this case, the simplest way to reach the exponential rate of convergence is to use the Hermite basis, by transforming each input parameter x i into a standard normal random variable ξ i :

ξ i = Φ -1 (F x i (x i )), (3)
where Φ and F denotes the cumulative density function (CDF) of a standard Gaussian variable and of the variable x i respectively. Given the vector of transformed (independent) input parameters ξ = (ξ 1 , ..., ξ n ), the PCE can be computed as:

y = f (ξ) = ∞ k=0 β k Ψ k (ξ). (4)
For computation, expansion (4) needs to be truncated, obtaining finally an approximation of the response as:

y ≈ P -1 k=0 β k Ψ k (ξ), (5)
Figure 1: Specification of the problem: the vertical section of the domain and the system of differential equations entered in Comsol Multiphysics.

-----ρ w (C p ∂p nw ∂t + (-C p + θ w K w

) ∂p w ∂t ) + ∇(-

K int k rw ρ w µ w ∇(p w -ρ w gz)) = 0, ρ nw (C p ∂p w ∂t -(C p + θ nw K nw
) ∂p nw ∂t ) + ∇(-

K int k rnw ρ nw µ nw ∇(p nw -ρ nw gz)) = 0.
Zero flux: n • ∇p w = 0, n • ∇p nw = 0.

Zero flux: n • ∇p w = 0, n • ∇p nw = 0.

Hydrostatic pressure p w = p w0 , p nw = p nw0 r = 0 The P coefficients β k , k = 0, ..., P -1 appearing in (5) are assumed deterministic and fully characterize the randomness of the response y. In applications, these coefficients are unknown and thus need to be properly estimated. A number of methods are presented in the literature [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF]; among them, in this work, the regression method is used. It consists in estimating β k by least-squares on the basis of an experimental design {y j , j = 1, ..., N }. The experimental design is built by collecting N exact solutions y j , j = 1, ..., N , which are derived solving the complete mathematical model for N realizations x j , j = 1, ..., N , of the random input x. The least-squares estimate of β = (β k ) ∈ R P is found by solving the minimization problem:

Q nw = 1M t/year, Q w = 0.
β = argmin β 0 ,...,β P -1 ∈R 1 N N j=1 y j - P -1 k=0 β k Ψ k (ξ j ) 2 , (6) being ξ j = (Φ -1 (F x i (x i ))). If the design matrix Ψ = (Ψ j (ξ i )) ∈ R N,P is of full rank, problem (6) admits the unique solution β: β = (Ψ T Ψ) -1 Ψ T y, (7)
where y = (y 1 , ..., y N ). Techniques for an efficient computation of coefficients estimates ( 7) has been presented in the literature (e.g., [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]) and are exploited in the simulations here presented.

The utility of PCE for the uncertainty analysis of the response is two-fold:

• Response mean and variance can be directly obtained from (5) as:

E[y] = β 0 ; Var(y) ≈ P -1 k=0 β 2 k ψ k 2 .
Analogously, the sensitivity Sobol' indices can be computed directly from coefficients β k , k = 0, ..., P -1 (see Subsection 3.2)

• The approximation (5) can be used as metamodel for Monte-Carlo simulations (i.e., the response is simulated through the approximate model instead of solving the PDE system). This allows to perform a much more efficient assessment of the CDF as the approximated model is easier to simulate compared to the full model (1).

Sensitivity Analysis through Sobol Indices

Sobol' indices [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]) will be used in Section 4 to perform the sensitivity analysis on the response. The total variance of the response can be decomposed into the effect of each random parameters x j and the joint effects of all the subsets x i 1 , ..., x is of the input parameters x 1 , ..., x n . As stated in [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF][START_REF] Formaggia | Global sensitivity analysis through polynomial chaos expansion of a basinscale geochemical compaction model[END_REF], this variance decomposition follows from the Sobol' decomposition of the function f (ξ) which appears in (4):

f (ξ) = f 0 + N i=1 f i (ξ i ) + N i,j=1 f i,j (ξ i , ξ j )+ (8) + • • • + f 1,2,...,n (ξ 1 , ξ 2 , . . . , ξ n ),
where

f 0 = Γ f (ξ)ρ Γ (ξ)dξ, f i (s) (ξ i (s) ) = Γ -i (s) f (ξ)ρ Γ -i (s) (ξ -i (s) )dξ -i (s) - I⊆i (s) f I denoting with i (s) the multi-index {i 1 , . . . , i s }, with ξ i (s) the collection {ξ i 1 , ..., ξ is }, with ρ Γ = n i=1 ρ Γ i the multivariate gaussian density of ξ, Γ -I = × i∈{1,...,N }\I Γ i and ρ Γ -I = i∈{1,...,N }\I ρ Γ i .
Notice that in the considered case Γ i = R, i = 1, . . . , n, but what follows remains valid for any measurable Γ i ⊆ R. In particular, decomposition ( 8) is unique whenever f is integrable over Γ and each term f i (s) in ( 8) is orthogonal with respect to the others.

To define the Sobol' indices, call V f , V f the variance of f , f respectively. The Sobol' index relative to the mixed effect ξ i (s) (or to x i (s) ) is defined as:

S i (s) = 1 V f Γ i (s) f 2 i (s) (ξ i (s) )ρ Γ i (s) dξ i (s) , (9)
where

Γ i (s) = Γ i 1 × • • • × Γ is . Hence, S i (s)
represents the proportion of the total variance explained by ξ i (s) (or x i (s) ) and thus the Sobol' indices sum to 1.

Concerning the total effect of a parameter ξ i (x i ) (or of x i ) on the total variability, it can be computed directly from Sobol' indices as:

S T i = i (s) :i∈i (s) S i (s) . (10)
The computational effort needed to compute the sensitivity indices ( 9) and (10) consists mainly in the cost of calculating (numerically) 2 n -1 integrals. To this end, Monte-Carlo quadrature schemes can be used but, due to its low rate of convergence, the computational cost might become unaffordable. However, this problem can be faced and solved by introducing the PCE of the response y, since Sobol' indices can be directly computed from coefficients β k , k = 0, ..., P -1. Indeed, expansion (5) can be reordered so that it is equivalent to (8):

f (ξ) = i (s) ⊆{1,...,n} k∈K i (s) β k Ψ k (ξ), (11)
where

K = {0, . . . , P -1}, K i (s) = {k ∈ K| Ψ k (ξ) = Ψ k (ξ = ξ i (s) )} (Crestaux et al. 2009).
From the equivalence between ( 4) and ( 11) the bijective relation existing between Sobol indices and coefficients β k , k = 1, ..., n can be derived. Indeed, the PC-based Sobol' indices can be obtained as:

S i (s) = k∈K i (s) β 2 k Ψ k , Ψ k P -1 k=0 β 2 k Ψ k , Ψ k
, being i (s) ⊆ {1, . . . , n}, while the total indices can be computed from (10).

SIMULATION AND RESULTS

After a preliminary study on the analytical model, the following four input parameters of the model have been considered: the maximal relative gas permeability kr max [-], the reservoir thickness H[m], the reservoir total porosity φ [-] and the irreducible water saturation (or maximal gas saturation Snw max [-]). These four parameters form the random vector x appearing in (1): x = (kr max , H, φ, Snw max ). In this study the log-normal probability density function is attributed to each parameter.

The response variable which is of interest for the sensitivity analysis is the natural logarithm of the maximal spread of the gas cloud, ln(r max ). Once a set of input parameters x 0 = (kr 0 max , H 0 , φ 0 , Snw 0 max ) is fixed, the corresponding ln(r 0 max ) can be calculated by numerically solving the system of partial differential equation governing the model (Figure 1). The solution includes the saturation 2D profile for any time step: ln(r 0 max ) is calculated as the maximum horizontal extension of the 0.99 saturation contour. The numerical solution of the system for a given set of parameters requires 5.67 minutes on a personal computer (4 CPU Xeon E31245 3.30GHz, RAM 16GB).

The experimental design In order to determine the least-square estimate of vector β as in ( 6), an experimental design {ln(r j max ), j = 1, ..., N } of size N = 100 has been considered. The number of replications N has been tuned by balancing the computational cost and the precision in least square estimate: the higher N is, the higher either the precision in estimating β or the computational cost are. In particular, the computational cost spent for generating an experimental design of size N on a personal computer with the technical characteristic specified before is about 5.67 • N minutes (O(N )), which corresponds to 9 hours and 27 minutes in the considered scenario.

In order to justify the choice of N , a preliminary study has been performed on the analytical model corresponding to the simplified physical problem presented in [START_REF] Nordbotten | Similarity solutions for fluid injection into confined aquifers[END_REF], for which the exact expression of the response distribution can be computed. In Figure 3 the empirical CDF of the output samples on which the meta-model is built (solid green line) is compared to the exact CDF calculated (c) N = 1000 analytically (dashed blue line), for N = 50, 100, 1000 (subfigures (a), (b) and (c) respectively). Even the meta-model built out of 50 samples reproduces the response variability quite closely; for N = 100 the approximation is very accurate without showing significant gains with respect to N = 1000. Therefore, a reasonable value for the parameter N is N = 100: this choice seems to balance the most the precision in the CDF estimate and the corresponding computational cost.

Supported by these results, the PCE on the gas injection model -detailed in Section 2-has been performed on the basis of an experimental design of size N = 100 (Figure 4). In Figure 4 the N = 100 realizations of the (log-transformed) response r max are represented against the parameters kr max , H, φ and Snw max , superimposing the corresponding least squares lines. The slope of the regression line is positive in the case of parameters kr max , while it is negative for H, φ and Snw max . This result can be expected from the point of view of the physical model: the higher the relative permeability of the gas is, the further it can penetrate, that is why the first correlation is positive. Concerning the other three parameters, the correlation is negative. Indeed, the increase of the reservoir height is not favorable for the horizontal spread. Increasing the porosity and maximal gas saturation leads to a bigger fraction of gas staying in place after the front passage, which means that for the same injected volume the maximal horizontal extension would be smaller.

PCE expansion and sensitivity analysis

The coefficients to be included in PCE has been selected following the sparse PCE approach proposed in (Blatman and Sudret 2010, [START_REF] Blatman | Adaptive sparse polynomial chaos expansions for uncertainty propagation and sensitivity analysis[END_REF]. Having fixed the maximum degree p of Hermite polynomials to p = 4, the 15 coefficients β 0 , β 1 , β 2 , β 3 , β 4 , β 5 , β 7 , β 9 , β 15 , β 18 , β 20 , β 38 , β 40 , β 48 , β 58 , proved to be significant for the expansion. Only the first 5 coefficients are reported in Table 1, since the other coefficients proved unstable when replicating the simulations. Noticed that coefficients β 1 , β 2 , β 3 , β 4 reflects the behav- These estimates are coherent with a standard Monte Carlo simulation: indeed, 1000 Monte Carlo simulation provides an estimate of 7.2197 for the mean (error: 0.05%), and of 0.0863 for the standard deviation (error: 0.9%)).

In order to assess the uncertainty propagation through the system, Sobol' indices have been computed by using ( 11) and ( 10); the results are reported in Table 2. According to Sobol' total indices (Table 2, second column), the parameters which seem to influence the most the response variability are the reservoir total porosity Φ and the maximal gas saturation Snw max , whose variability explains about 31.0% and 29.8%, respectively, of the output variability. The remaining variability is explained in decreasingly order of importance by maximal relative gas permeability kr max (20.2%) and reservoir thickness H (18.2%). Very similar conclusions can be drawn by considering the Sobol' indices reported in the first column of Table 2, meaning that the influence of the input parameters on the output is mainly due to single effects, while mixed effects play a minor role.

To access the uncertainty in the estimated statistics, the simulation has been replicated with the same settings starting from different seeds. The results are 3. It can be noticed that the response mean and standard deviation seem precisely estimated. Instead, the parameters affected by the most significant variability are the Sobol' indices relative to the reservoir thickness H and to the maximal gas saturation Snw max -either the indices relative to the single effects or the total indices-. Notice also that the coefficients β k selected by the sparse PCE approach varied by replicating the simulation, with only the first five coefficients steadily selected. However, both the variability in the PCE estimates and in the selected coefficient do not seem to significantly affect the conclusions previously drawn about the uncertainty propagation through the system.

Assessing the response distribution function Although Sobol' indices provide an important tool for sensitivity analysis, in this study the whole CDF is of interest. Generally, CDF can be estimated in three ways: (a) analytical calculation (b) standard MC simulation (c) meta-model MC simulation -i.e., generation of an experimental design of size N by a standard MC-simulation on the complete model, computation of the β coefficients in approximation ( 5) and finally MC simulation on the meta-model exploiting the same approximation-. The main drawback of method (a) is that an analytic form of response CDF is rarely available; on the other hand the computational cost of a standard MC simulation on a complex model as the one considered here increases fast with the number of samples required (more than 94 hours for 1000 samples). Notice in particular that a preliminary MC simulation is needed when using method (c), but a fair size of the experimental design is typically much smaller than the number of simulations required to assess the complete CDF with method (b).

To evaluate the performance of methods (b) and (c) when the analytical form of the CDF is available, the simplified gas-injection model has been again considered. In Figure 5 the following CDFs are compared: 1) exact CDF -method (a), dashed blue line-; 2) empirical CDF obtained by 1000 standard MC-simulations -method (b), green line-; 3) empirical CDF obtained by 10000 PCE-based meta-model simulations having an experimental design of size 100 -method (c), red line-. The parameters of the latter simulation have been chosen according to the feasibility of an analogous simulation on the complete gas injection model: indeed, the computational cost of the MC-simulation on the complete model soon becomes unaffordable as the sample size increases and the simulation of 1000 samples already requires more than 94 hours 1 .

The empirical CDFs obtained with the two methods are pretty similar to the exact CDF, but more accurate and smooth results are obtained for PCE-based simulations, being based in fact on a rich meta-model MC sample. Figure 6 represents the statistical moments of the response. Notice that the meta-model reproduces exactly the even statistical moments even when it is build on a small size of the experimental design. Theoretically being zero, the odd statistical moments seem to be better approximated for metamodels based on a larger size of the experimental design. Observe that the time consumption required for computing the meta-model estimate is essentially due to the computation of the experimental design. Therefore the very low additional cost required for the PCE step turns in a much more accurate CDF estimate.

Supported by these results, methods (b) and (c) have been applied on the full numerical model of gasinjection detailed in Section 2. The response empirical CDF computed by 1000 standard Monte Carlo simulations (blue dotted line) and with a PCE-based meta-model simulation (red line) are presented in the left panel of Figure 4, superimposed to the results of the replicated simulations (grey lines). It has to be noticed that the results are fairly stable among the replicates, as the lines seems perfectly superimposed. This supports the conclusion previously drawn about the stability of the results with respect to the variability in the selected coefficients and their estimates. Moreover, all the empirical CDF are very close to the MC results obtained from the full model.

In order to appreciate the most the existing differences among the obtained estimates, the right panel of Figure 4 reports the kernel density PDF estimates computed as before from the MC sample, the metamodel sample and the replicates. Replicates have been built either with the same parameters (darker lines), or by fixing the maximum degree of Hermite polynomials to p = 5. Although some differences are now appreciable, all the meta-model curves appear very similar each other and very close to the MC estimates. The main difference is present in the right tail of the distribution, where a slight underestimation occurs. However, the meta-model PDF estimates provides a very good approximation of the MC PDF estimates obtained in less than 10 hours against the 94 hours required by the full MC-simulation.

Finally, Figure 8 shows the statistical moments of the resulting response distributions. All the replicates 1 4 CPU Xeon E31245 3.30GHz, RAM 16GB seem to provide similar results, that are close to the MC estimate. No significant gains are reached if increasing the maximum degree of the Hermite polynomials. The first two moments of the distribution, in particular the mean, are very precisely estimated, as already noticed when analyzed Table 3. Increasing uncertainty is observed when increasing the order of the estimated moments. As in the simplified model, the uncertainty in the estimates is higher for even order moments, while odd moments are more accurately estimated.

CONCLUSION

In this work, the Polynomial Chaos Expansion method as an efficient technique for performing uncertainty and sensitivity analysis of complex numerical models has been investigated. In particular, the potentiality of this approach has been explored in terms of computational efficiency when dealing with a model of gas injection into porous media. The regularity of the response for this particular physical model is favorable for the use of the described method.

The crucial points encountered while performing PCE on a real system have been underlined. For instance, one need to pay close attention when making the choice of the size of the experimental design used for estimating the expansion coefficients and the basis functions form. The minimal size is dictated by the number of coefficients in the decomposition, further validation of the meta-model quality can only be done by experience in every particular case. Here, the size of the experimental design has been fixed by exploiting the information coming from the application of the same technique to a simplified model, but other considerations -e.g., in terms of variance of the estimate-can be made.

A great deal of attention has been focused on the computational gain of the PCE approach, which is very significant for the considered application: while obtaining the response PDF by 1000 standard Monte Carlo simulation takes more than 94 hours, with a PCE approach the computational cost is reduced to less than 10 hours.

However, the robustness of the procedure with respect to the size of the experimental design would be worth being further investigated. Indeed, by repeating the simulations with a different size of the experimental design, the coefficients estimate might significantly vary, with possible effects on the Sobol' indices and on the consequent sensitivity analysis. Although this does not seem to be the case, in order to avoid this problem, a direct control on the estimates variance should be considered, in particular if no 'a priori' information is available on the response distribution.

A further development regarding the considered physical model, which might be of industrial and physical importance, is to consider a space-time varying response, such as the water and gas pressure, while analyzing the propagation of the uncertainty from the input parameters to the space-time random field described by the response. In such a complex setting, the use of PCE could make the computational cost of performing uncertainty and sensitivity analysis affordable, which is a very important goal to be reached for this kind of complex numerical systems.

ACKNOWLEDGMENT

The authors would like to thank the organizers of ISRI seminar. Without their effort for bringing researchers from different disciplines together, this work would not have been possible. The event is organized under the companion of Supélec, École Centrale and Politecnico di Milano. The authors also thank MSS-Mat laboratory for providing the numerical resources and technical support. The work of E. Bas ¸tug has been supported by Alcatel Lucent Chair on Flexible Radio and the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering).

Figure 2 :

 2 Figure 2: Numerical solution of the system for t = 10 years. The color mapping shows the water saturation level. -2350 -2380 -2410 -2440 -2470 1.00 0.75 0.50 0.25 0.00 0 300 600 900 1200 depth (meter) horizontal distance from the injection well (meter)

Figure 3 :

 3 Figure3: Comparison of the meta-model CDF (red line) with the exact CDF (dashed blue line) and with the empirical CDF of the samples used for the meta-model construction (green line). The figures correspond to the analytical solution with varying number of the samples for PCE decomposition, from left to right: 50 samples (first panel), 100 samples (second panel), 1000 samples (third panel).

  Figure 4. Afterward, the estimates of the first two moments of the response probability density function (PDF) have been derived by using (5), obtaining: E[ln(r max )] ≈ 7.2157 sd[ln(r max )] ≈ 0.0871.

Figure 4 :

 4 Figure 4: Experimental design. From left to right: plot of r max against kr max (first panel), H (second panel), φ (third panel) and Snw max (fourth panel).

  Figure 5: Comparison of the meta-model CDF obtained with an experimental design of size N = 100 (red line) with the exact CDF (dashed blue line) and with the empirical CDF obtained with 1000 standard MC-simulations (green line)

Figure 7 :

 7 Figure 7: Response CDF (left panel) and kernel density estimated PDF (right panel) for: 1000 standard Monte Carlo simulation (dotted blue line), meta-model simulation through PCE (solid red line), replicated meta-model simulation (dark grey lines).

Figure

  Figure 8: Comparison of the statistical moments for: 1000 standard Monte Carlo simulation (dotted blue line), meta-model simulation through PCE (solid red line), replicated meta-model simulations when the maximum degree of the Hermite polynomials is p = 4 (light grey lines) and p = 5 (solid grey lines).

Table 1 :

 1 Estimates of the decomposition coefficients

	Parameter Estimate
	β 0	7.2157
	β 1	0.0387
	β 2	-0.0370
	β 3	-0.0485
	β 4	-0.0475
	ior previously noticed by inspecting

Table 2 :

 2 Sobol' indices

	Parameter Sobol' index Sobol' total index
	kr max	0.2019	0.2098
	H	0.1823	0.1884
	Φ	0.3102	0.3177
	Snw max	0.2977	0.2990

Table 3 :

 3 Mean, standard deviation and coefficient of variation relative to the estimated indices.

		Mean	Sd	C.V.
	E[ln(r max )] 7.2151 0.0020 2.8e-4
	sd[ln(r max )] 0.0861 8.5e-4 0.0099
	S krmax	0.2026 0.0087 0.0431
	S H	0.1640 0.0179 0.1091
	S Φ	0.3119 0.0104 0.0335
	S Snwmax	0.3149 0.0214 0.0681
	S T krmax S T H S T Φ S T Snwmax	0.2079 0.0082 0.0395 0.1701 0.0236 0.1384 0.3159 0.0100 0.0318 0.3168 0.0192 0.0608
	reported in Table


N Lafrance

M Souley

C Auvray

O Favreau

G Eo R Esso

V Labiouse

A g e in g o f ch alk rock s in an underground quarry

In the present work, an attempt has been achieved to describe the effects o f ageing and of water on the hydro-mechanical behaviour of chalks from underground quarries of northern Parisian basin. Physical and mechanical properties were obtained as a function of distances from the pillar face (horizontal core drillings, from the wall to the centre) and relative humidity (Hr) of air. It has been stated that physical and mechanical properties change with depth in the pillar: the density and compressive strength of chalk tend to decrease from the pillar face (at least 1 m from the wall) to the centre, contrary to the classical observations on the ageing of geomaterials. Uniaxial and triaxial compression tests under controlled relative humidity put into light the sensitivity to water (and then relative humidity in the quarry) on the mechanical strength and deformability of chalks.

INTRODUCTION

Chalk rocks have been used during centuries as an easily crafted and strong building material. The shal lower underground room-and-pillar quarries where the chalk rocks were extracted are now abandoned and endure the effects o f time and weather, which poten tially increase their risk o f collapse. In order to prevent this risk, INERIS studies the rock behaviour o f two partially flooded -facing water table fluctuationsunderground quarries in northern France (Estreux and St-Martin-le-Noeud) through in situ instrumentation implemented in 2004 and 2009 respectively [START_REF] Gombert | In-situ and labo ratory tests to evaluate the impact o f water fluctuations on stability o f underground chalk mines[END_REF]). The two quarries located at an average depth of 30 m, have been exploited by the room and pillar method with extraction rateof 60-75%, and then abandoned for more than a century. The original goal o f the in situ characterization was to determine the impact o f groundwater table fluctuations on their sta bility. These fluctuations lead to chemical water-chalk interaction and impact the chalk mechanical behaviour as a consequence of climate change action. Laboratory characterization is also undertaken. Nevertheless, its hydro-chemo-mechanical behaviour is not well known despite several previous studies. Due to complex inter locked phenomena acting on the mechanical behaviour o f chalks, it remains necessary to address the issue with multi-scale and multi-physics approaches.

The present work is devoted to the experimental characterization o f the ageing on the hydro mechanical behaviour o f chalk collected in under ground quarries and water-chalk interactions when linked to water flooding/unflooding cycles.

MATERIALS AND METHODS

Tested chalk

The studied samples came from one underground chalk quarry in the Parisian Basin, from Santonian age, in the township of Saint-Martin-le-Noeud near Beauvais (Oise, Picardie, France). The work is achieved by studying samples from one pillar. The studied chalk can be regarded as a dolomitic chalkdolomite proportion occupying nearly 14% o f its mineralogical composition -and is homogeneous and isotropic with a high porosity (approximately 40%).

Sample preparation

All o f the samples were extracted by a corer (oriented horizontally) with water as lubricant and cut once extracted perpendicularly to the stratigraphy by a saw guaranteeing parallel and smooth surfaces. All o f them were sized up by a high-grade electronic calliper and placed at 60°C to constant weight to guarantee no dam age and no more water in the rock. The samples were then instrumented in order to measure axial and lateral strains. Strain gauges were gluedby pairs to the sample using the same material. The gauges -lateral or axialwere positioned opposing each other half the length of the sample. Fully saturated samples were obtained by placing the dried samples under equilibrate water to constant weight. Partially saturated samples were obtained by placing the dried samples under imposed hygrometry using saline solutions. The sampling of the studied pillar could be decomposed into two parts:

• "ageing sampling": constituted by horizontal 120mm-diameter core drillings SCI (z = 2.15 m from floor) and SC3 ( z = 1.70 m from floor). The cylindrical samples (h ^ 40 mm, a ^ 20 mm) were taken perpendicularly to the core drilling and placed in fully saturated and dry states • "water-effect sampling": composed by blocks from the pillar face towards the heart of pillar. The cylindrical samples (h ^ 100mm, a ^ 50 mm) were placed in dry, fully saturated and partially saturated states.

Physical properties

Mechanical behaviour is sorely dependent on its physical properties such as densities.

Dry pi, fully saturated p\, and partially saturated Per density was calculated knowing the volume Vt and the dry ma, fully saturated mh and partially saturated mm masses of the sample.

For the "chalk ageing characterization" two tech niques have been used in order to obtain an idea o f the porosity. The easier way consisted in a mass compar ison between the saturated and the unsaturated states o f a same sample. The second way was to calculate porosity by use o f mercury porosimeter.

CONTROLLED HYGROMETRY TESTING

Methodology and procedure

The aim o f these tests was to determine the influ ence o f water on the mechanical behaviour of the chalk. In order to characterize the water effect on the dolomitic chalk behaviour, mechanical tests were per formed on 1 to 3 samples for several degrees o f relative humidity (of air) ranged from 90% to 98%; values typ ically encountered in situ and for different confining pressure values (0, 1 and 3 MPa). Relative humid ity Hr was enforced using saline solutions: at 20°C zinc sulphate heptahydrate (ZnS0 4 , 7 H2O), potassium sulphate (K2SO4) and copper sulphate pentahydrate (CUSO4, 5H2O) were used to obtain respectively 90%, 95% and 98% degrees of relative humidity. Uniaxial and triaxial compression tests under imposed hygrom etry conditions were carried out on samples equipped with strain gauges. During uniaxial compression tests the sample is placed in a steel-based cell above the salt solu tion (Figure la) and separated from the exterior by plastic film. Triaxial compression tests were made using an imposed hygrometry autocompensated tri axial cell under drained conditions (Figure lb ) and a lOOOkN capacity press. The tests were performed at 48 |xm/min with 1 and 3 MPa as a confining pres sure. Unloading/loading cycling started every 2.5 kN. Finally mechanical tests were performed under con trolled strain rate.

Results

Table 1 summarizes the results about mechanical prop erties o f the chalk depending on the hydric state of samples (from dry to saturated) for uniaxial com pression tests. These properties are then correlated to densities. Table 2-3 summarizes the results about mechanical properties o f the chalk depending on the relative humidity for triaxial compression tests. The tests were carried out at 1 and 3 MPa o f confining pressure (03). These properties are also linked to den sities. To illustrate the typical mechanical behaviour of the dolomitic chalk, Figure 2 presents the stressstrains curves for three uniaxial compression tests for dry, partially (Hr = 90%) and fully saturated samples. Axial strain curve is in black, lateral strain curve in The objectives o f this series o f tests were to eval uate the sensitivity o f the dolomitic chalk when in contact with water. Although further tests are still ongoing at the time of writing, it appears that for uni axial compression tests, the compressive strength of the dolomitic chalk decrease with wetting (Figure 2) and presents a seemingly abrupt transition at the last stages o f wetting to complete saturation (Table 1). This evolution o f Sc according to Hr is identical to those encountered in the literature for some geoma terials, particularly for chalks as shown in Figure 3 (Watelet, 1996). However static Young modulus E doesn't seem to evolve with water content consider ing uniaxial compression tests alone -like Poisson's ratio v. There is no loss o f its stiffness when partially then fully saturated with water. presented herein) -there is no plastic stage between elastic domain and failure (pic strength) (Figure 2a). This characteristic can be linked with its mineralogy, dolomite crystals tougher (in terms o f deformability) than calcite crystals and less sensitive to water. However the range o f this plastic stage -almost non-existent at dry state -increases with wetting. With water, the dolomitic chalk tends to present a more ductile mechanical behaviour, like the previously mentioned chalks (Lixhe or Estreux). Generally speaking, for a given Hr degree (0, 90%, 95%, 98% and 100%), compressive strength increase with the confining pressure (Table 23). Foragivenuj, variation of Sc with respect to Hr is qualitatively identi cal to that discussed under uniaxial conditions (except the test: <73 = 3 MPa and H r= 100%). The exception of test at 3 MPa of confining pressure under saturated conditions can be explained by the fact that water puts pressure onto the porous media and as a result there is no "vacuum" in the sample.

Finally, for the rest of unsaturated state, the mechan ical behaviour under triaxial conditions follows the same trend as the unconfined ones.

AGEING EFFECTS TESTING

Methodology and procedure

In order to characterize the ageing of chalk -in sat urated and unsaturated states -SCI and SC3 core drillings were mapped: each sample was extracted depending of its range inside the pillar. This methodol ogy has been previously applied on gypsum [START_REF] Auvray | Vieillissement el comportement rhéologique du gypse[END_REF][START_REF] Auvray | The aging o f gypsum in underground mines[END_REF]) and iron ore [START_REF] Grgic | Ageing o f Lorraine (France) abandoned iron mines. Rock Mechanics: a chal lenge for society[END_REF][START_REF] Grgic | Ageing o f Lorraine (France) abandoned iron mines. Rock Mechanics: a chal lenge for society[END_REF]. The sampling collection went every 10 cm from the face to the heart o f pillar. Then for each depth, physical properties such as density and matrix structure are obtained. It gave us the distribution of physical and mechanical properties inside the pillar depending on the distance (from the pillar face to the centre) and the stratigraphy.

SEM analysis is performed on mirror polished samples at the depths of 0, 120 and 240 cm for the SC3 core drilling. Mercury porosity is also measured at the same distances.

Physical -such as density -and mechanical prop erties are then presented in a diagram depending on the depth in the pillar.

Because o f the size o f the samples, uniaxial compression tests were carried out using a 5kN capacity press under a constant displacement rate of 30 |i.m/min.

Results

Figures 4-5 show how chalk physical and mechanical properties evolve depending on the depth in the pil lar. Each point represents one sample. The diagramscompressive strength (Figure 4 showing difference between the edge and centre o f pillar) -are presented as follow: i) SCI core drilling samples in dry state: the dolomitic chalk is more resistant at the edge than the inside pillar; ii) SC3 core drilling samples in dry state: it can be subdivided into two parts (depths ranged between 0 and 100 cm where Sc presents a very slight decrease and range from 100 and 300 cm where Sc decrease is more evident); iii) SC3 core drilling sam ples in a fully saturated state: the same as SCI core drilling dry samples.

The same goes for density diagrams. The scattering of measurements is also more significant at the edge than at the inside pillar.

Table 4 summarizes mercury porosity for distances of 0,120 and 240 cm from the pillar edge, in the case of SC3. Similarly porosity increases from the inside pillar to the edge for SC 1. No clear trend shows on Figure 5 for SC3 dry samples density. For SCI dry samples and SC3 saturated samples, a tendency is however clearly shown.

Figure 6 shows the difference between the distances 0, 120 and 240 cm for SC3 regarding the microstructure o f the rock matrix. The scale on each SEM photographs represents 2.5 //m. The crystals are well defined -almost holomorphic -with different sizes at the pillar wall (at the pillar wall) and rougher at the inside pillar (depth o f 240 cm) -as if they were eaten up -with almost identical sizes. The depth of 120 cm can be seen as intermediary.

Discussion

Comparing SEM photographs (Figure 6) with Table 4, the mercury porosity and rock matrix microstructure are well associated. The crystals were dissolved pref erentially at the centre than at the wall where chalk is exposed to air transport. 1 lowever, results on Figures 3 4 according to other geomaterials like gypsum and iron ore, tends to show a very liglil contrary trend. Ageing o f chalk seems to be of no standards. Bui the tendency is not that clear for SC3 dry samples compared to SCI dry and SC3 fully saturated samples. For iron ore ageing is partly shown by redox reactions (Grgic el al. 2013[START_REF] Ci | Mineralogical and microtextural characterisation o f the anthropic-origin age ing o f iron ore in Lorraine (France)[END_REF]. For gypsum [START_REF] Auvray | The aging o f gypsum in underground mines[END_REF], Hoxha et al. 2006), ageing is partly shown by traces of dissolution. Forllic dolomitic chalk dissolution traces (third photo on Figure 5) are well defined but occur al the inside pillar instead o fth ee d g e-asn o t seen for gypsum. The parting in the results shown in Figure 4 for SC3 dry samples may be explain by micro-cracking phenomena with an increase of the porosity.

The compressive strength tends to be weaker where the dolomitic chalk seems to be more compact. This possible contrary trend could be due to more impor tant water (liquid) circulation at the centre than in the vicinity of pillar wall: based on in sirii measure ment [START_REF] Gombert | In-situ and labo ratory tests to evaluate the impact o f water fluctuations on stability o f underground chalk mines[END_REF]) relative humidity cyclically fluctuates at the edge and stabilizes at the centre of the pillar: the edge and the inside pillar are not sub jected to the same hydric conditions. The defraying of electrostatic bondings -connecting crystals to each other -is one o f the numerous phenomena governing the degradation of the mechanical behaviour o f chalks (Shroeder 2002, Shroeder et al. 2003,1 lomand & Shao 2000). It generally occurs when chalk is in contact with the polar molecules of water.

As samples are small, a possible scale effect is not to be spread, especially in order to explain the parting in the results of densily and compressive strength. The tendency is also nol totally clear.

The hypothesis risen to explain the clear distinc tion between edge -regarding the matrix structureis that the main water transport -with the highest velocity -occurring in the pillar is gravity, but il may not exclude the very-low-velocity transport of water occurring transversely. The main mechanism may tend then to be more of the chemo-hydric type Ilian the chemo-mechanical-hydric ones.

CONCLUSION

Relating to water effect on chalk behaviour, uniaxial and triaxial compressive lests under imposed humid ity conditions highlight the sensitivity to water of the dolomitic chalk and ils deformability depending on relative humidity in lire quarry. According to [START_REF] Gombert | In-situ and labo ratory tests to evaluate the impact o f water fluctuations on stability o f underground chalk mines[END_REF], it could be expected a greater impact by climate change with its action on ageing in abandoned mines. Uniaxial and triaxial compression testing with other relative air humidity degree s is ongoing, and goes from degrees 36% to 85% in order to complete the data base of Saint-Martin-le-Noeud site.

Relating to ageing of chalk, although further tests are still ongoing at that time, initial results suggests that water (liquid ) circulation is more important at the centre than in the vicinity of pillar wall. At this step of the study, the possibly contrary tendency shown on the results has yet to be substantiated: complemen tary studies on other core drillings and parameters such as dynamic elastic modulus -with the help of speed velocity -must be characterized to com plete the results. Nevertheless, ageing o f chalk rock in underground quarries seems to be o f the 11011standard-type and sent more apart chalk rocks from other geomaterials.

Figure 1 .

 1 Figure 1. Experimental setup o f controlled relative air humidity for mechanical tests (a) uniaxial compressive test (b) triaxial compressive cell. Temperature at 20°C.

Table 2 .

 2 Results o f triaxial compression and physical prop erties o f chalk ( 0 3 = 1 MPa). E and v are calculated using cycles in the elastic domain and are presented here as average valuesvalues with low standard deviation. dark grey and volumetric strain curve in light grey. The scale was standardized to compare all the curves. 3.3 Discussion Physical properties of the dolomitic chalk can be called standards similarly to other types of chalk such as Hardivilliers chalk (Hjuler et al. 2009) and Haubourdin chalk (Nadah et al. 2013).

  The dolomitic chalk exhibits a brittle type mechanical behaviour. Compare to other types of chalks -such as Lixhe (Shroeder et al. 2003, Homand & Shao 2000) or Estreux chalk (not a (MPa) a (MPa)

Figure

  Figure 2. Stress-strains curves ofuniaxial compression tests as a function o f hydric state o f sample (a) dry, (b) partially saturated (H r= 90%), (c) fully saturated.

Figure 3 .

 3 Figure 3. Compressive strength ratio versus water satu ration degree for different geomaterials such as limestone, sandstone or chalk from Watelet (1996).

Figure 4 .

 4 Figure 4. Ageing uniaxial compressive strength versus deplli (from pillar face) on dry and/or saturated samples of SCI and SC3 core drillings.

Figure 5 .

 5 Figure 5. Densities versus depth from pillar edge on dry and or saturated samples o f SCI and S( 3 core drillings.

Figure 6 .

 6 Figure 6. SEM analysis on m irror polished samples illus-I rating the difference between the edge and centre o f the pillar (SC3 core).

Table 1 . Results o f uniaxial compression testing and physi cal properties o f chalk. E and v are calculated using cycles in the elastic domain and are presented here as average values. Mechanical properties Physical properties Hr %

 1

		Sc	E*	v*	Pi	Ph	PHi
		MPa	GPa	( -)		g.cm -3
	Dry	11	10.0 0.18	1.53 _	_
	90	5	9.9 0.20	1.54 -	1.55
	95	7	7.6 0.24	1.54 -	1.55
	98	5	10.9 0.21	1.69 1.91	1.97
	Saturated	2	11.3 0.19	1.53 1.88

♦average values with low standard deviation.