Jessica Franco

Laurent Carraro

Olivier Roustant

Astrid Jourdan
email: astrid.jourdan@eisti.fr.

A RADAR-SHAPED STATISTIC FOR TESTING AND VISUALIZING UNIFORMITY PROPERTIES IN COMPUTER EXPERIMENTS

Keywords: Computer experiments, Space-Filling Designs, Dimension Reduction, Discrepancy, Kolmogorov-Smirnov Statistic, Cramer-Von Mises Statistic

In the study of computer codes, filling space as uniformly as possible is important to describe the complexity of the investigated phenomenon. However, this property is not conserved by reducing the dimension. Some numeric experiment designs are conceived in this sense as Latin hypercubes or orthogonal arrays, but they consider only the projections onto the axes or the coordinate planes. In this article we introduce a statistic which allows studying the good distribution of points according to all 1-dimensional projections. By angularly scanning the domain, we obtain a radar type representation, allowing the uniformity defects of a design to be identified with respect to its projections onto straight lines. The advantages of this new tool are demonstrated on usual examples of space-filling designs (SFD) and a global statistic independent of the angle of rotation is studied.

INTRODUCTION

For the last 15 years or so, the design of experiment theory initiated by [START_REF] Fisher | The arrangement of field experiments[END_REF] has experienced a revival for the analysis of costly industrial computer codes. This development has led to at least two major changes. First, these codes represent phenomena of an increasing complexity, which implies that the corresponding models are often nonlinear and/or nonparametric. Second, the experiment itself is different. Numerical experiments are simulations and, except for stochastic codes implementing a Monte Carlo-based method, produce the same response for identical conditions (including algorithm and computer-based parameters). Therefore, repeating an experiment under the same conditions does not make sense since no new information is acquired.

In this new paradigm, the experiment planning methods are different. For example when the code is to be analyzed for the first time before any simulation has been made (scanning phase), one often tries to satisfy the following two requirements. Firstly, distribute the points in the space as uniformly as possible to catch non-linearities; this excludes repetitions also.

Secondly, this space coverage should remain well-distributed even when the effective dimension is lowered. The first requirement was the starting point of research work in space filling designs (SFD). The quality of the spatial distribution is measured either by using deterministic criteria like minimax or maximin distances (Johnson, Moore and Ylvisaker, 1990), or by using statistical criteria like discrepancy [START_REF] Niederreiter | Low-Discrepancy and Low-Dispersion Sequences[END_REF][START_REF] Hickernell | A generalized discrepancy and quadrature error bound[END_REF][START_REF] Fang | Design and Modeling for Computer Experiments[END_REF]. The second requirement stems from the observation that codes often depend only on a few influential variables, which may be either direct factors or "principal components" composed of linear combinations of these variables. Therefore dealing only with these influential factors is more efficient. Note that dimension reducing techniques like SIR [START_REF] Li | Sliced inverse regression for dimension reduction (with discussion)[END_REF] or KDR [START_REF] Fukumizu | Dimension Reduction for Supervised Learning with Reproducing Kernel Hilbert Spaces[END_REF] effectively identify the space generated by the principal factors. Hence, it is desirable that the space filling property should be also satisfied in the projection onto subspaces. This is the aim of Latin hypercubes designs (LHD) and orthogonal arrays (OA) in computer experiments. For instance, the space filling LHDs provide a uniform repartition in the projection onto margins, so that there is no loss of information if the code depends on only one variable. In addition, orthogonal arrays with space-filling properties extend this aspect to higher dimensions (see, for example, [START_REF] Koehler | Computer Experiments[END_REF], and more generally, Owen, 1992[START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Nevertheless, considering only the projections onto margins is not sufficient if, for example, the code is a function of a linear combination of 2 variables.

In this article we introduce a statistic based on 1-dimensional projections to test the uniformity for a design of experiment. The choice of 1 dimension is due to the difficulty in computing the theoretical distribution for higher dimensional space. The advantage of restricting to a single dimension is that it offers a simple viewing tool similar to a radar screen. By representing the statistic's value in all directions, one obtains a parameterized curve (or surface), identifying the uniformity defects of a design with respect to its projections onto straight lines. The article is structured as follows. In section 2 we define the new statistic and the associated visualization tool, called uniformity radar, and give some properties of them. In section 3 we show examples of radar applications to space-filling designs. Section 4 is devoted to extend the concept by defining a global statistic which does not depend on a particular axis of rotation. In section 5, a discussion on the uniformity radar allows specifying its scope of application. Proofs are given in the appendix.

UNIFORMITY RADAR

In the analysis of a computer code, let us consider a uniform experiment design on a cubic

domain [ ] 1,1 d Ω = -
. Note 1 ,. μ is uniform. However, this is not realistic when a L is not a coordinate axis.

For example, in the case of the cubic domain [ ] 2

1,1 -, the density of the projected design points is higher in the central part of the axis as can be seen in Figure 1. In the general case, the projection onto a L is a linear combination of independent random variables of uniform distribution, which leads to a traditional problem of probabilities first 0

β α M -M -α
where :

M = |cosθ| + |sinθ| α = | |cosθ|-|sinθ| | β = α + M 1
solved by Lagrange in the 18th century (see discussion in [START_REF] Elias | Convolution of Uniform Distributions and Ruin Probability[END_REF] on this topic).

If we denote

( )

1 ,..., d X X X =
, where 1 ,..., d X X are independent random variables identically distributed following the uniform distribution over [ ] , then the distribution function of Z is given by :

F Z (z) = 1 2a j j=1 d ∏ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ × ε(s) (z + s.a) + d d! s∈{-1,1} d ∑ where ( ) { } 1 ,..., 1,1 d d s s s = ∈- are the corners of the hypercube [ ] 1,1 d Ω = - , 1 ( ) d j j s s ε = = ∏ , .
s a is the scalar product of the vectors s and a, and ( ) max( , 0)

y y + =
the positive part of y. As a result, for a given axis, Z admits a piecewise linear density whose nodes correspond to the projections of the corners of the domain.

Note. The distribution of the projections is known in other situations, as in the case of a spherical domain: if Ω is the unit sphere of R d , a direct calculation shows that a μ admits the density

( ) [ ] ( ) 2 1,1 2 1 1 a f x x x π - = - , the distribution function being equal to 2 1 1 ( ) (Arcsin 1 ) 2 a F x x x x π = + + - for [ ] 1,1 x ∈ - .
In sum, for a uniform experiment design to have good distribution properties on the 1dimensional projections, it will be necessary that in all the directions a, the empirical distribution of the projections onto a L is close to their theoretical distribution under the hypothesis ( ) 0 H . There exist many distribution adequacy statistics (see D'agostino and [START_REF] Stephens | Goodness-of-fit Techniques[END_REF], which allow for a large number of choices to define a criterion adapted for this purpose. However, possibilities are limited by special requirements. To start with, it is preferable that the statistic's distribution be known to avoid the approximate calculation of its distribution. Furthermore, one would also like the statistic to be distribution free, that is, its distribution doesn't depend on the projection direction to have a unique rejection threshold for all the angles. Also, for the sake of consistency, it would be desirable for the retained statistic to be interpretable in terms of discrepancy when projections are made onto a coordinate axis.

Finally, two famous statistics (at least) correspond to these requirements: the Kolmogorov-Smirnov statistic

D N (a) = sup ˆ F N ,a (z) -F a (z) (1)
and the Cramér-Von Mises statistic μ is the uniform distribution on [0,1] and these statistics correspond, respectively, to the discrepancies L ∞ and L 2 [START_REF] Niederreiter | Low-Discrepancy and Low-Dispersion Sequences[END_REF]. In what follows, we decided to work on the first because the conclusions seem equivalent, while the corresponding graphics are a little more readable (see section 5). By analogy with the case of coordinate axes, we will talk about a discrepancy of projections to designate the Kolmogorov-Smirnov statistic of the formula (1).

Nω N 2 (a) = ( ˆ F N ,a (z) -F a (z)) 2 dz ∫ (
The discrepancy of projections provides a tool for visualizing the defects in uniformity based on the 1-dimensional projections. Since this tool looks like a radar screen, we propose to call it uniformity radar. Its utilization depends on the dimension of Ω . In 2 dimensions, the discrepancy of projections is calculated in 360° by continuously projecting onto a rotating axis. Thus, one obtains a parameterized curve in polar coordinates directions, the 2D radar provides decision support on design uniformity. In 3 dimensions, one calculates the discrepancy of projections onto an axis , L θ ϕ , for all directions, pivoting around the center of the domain. This axis is defined in spherical coordinates by an angle θ in longitude and ϕ in latitude. This time a parameterized surface is obtained, called 3D radar,

( ) ( ) , , N D θ ϕ θ ϕ a defined over [ ] 0, 2 , 2 2 π π π ⎡ ⎤ × - ⎢ ⎥ ⎣ ⎦ .
In higher dimensions, it seems unrealistic to make an angular scan of the space Ω . In addition, it becomes almost unfeasible to represent the result graphically (although the calculation is still possible). However, the hypothesis (H 0 ) remains valid on 2-and 3-dimensional coordinate spaces. Therefore, the uniformity radar may be applied to all pairs and/or triplets of possible dimensions.

In practice, the quality of the representation can be degraded by discretizing. Here the 2D and 3D radars are continuous applications. But N D is not differentiable with respect to all the axes a L so that at least two points of Ω are projected onto the same point, which explains why the parameterized curves of the next section are not smooth and contain many singularities.

APPLICATIONS OF UNIFORMITY RADAR

In this section we present a few examples to show the interest in using uniformity radar to test the uniformity of the distribution of experimental points by looking at their projection.

We consider cases where the hypothesis ( )

0
H of a uniform distribution in the experimental domain is plausible. For each representation of the uniformity radar we have added the circle (or the sphere) of radius ks equal to the statistic of the Kolmogorov-Smirnov test associated with a 95% confidence level. Recall that since the statistic is distribution free, ks does not depend on a. This provides a decision-making support or, at least, a means of comparison with the random designs obtained by a uniform sampling. Should the studied design be stochastic (pseudorandom, Latin hypercubes or randomized orthogonal arrays), the graph displays the directions a along which the hypothesis ( ) 0 H is rejected. If the design is deterministic, we are outside the usual scope of application of the test. If one of the values of ( ) N D a is greater than ks, then we can only say that this design is worse than a random design since the probability that a random design will have a better discrepancy exceeds 95%.

The following examples apply essentially to the latter cases because we preferred to use known SFD designs without transforming them. Nonetheless, in practice, it would suffice to apply randomization or scrambling (see, for example, [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF] to obtain a stochastic design, and thus be under the usual assumptions of statistical tests.

Example 1. Analysis of a 15-dimensional Halton sequence using 2D radar. We consider the first 250 points of a 15-dimensional Halton sequence of low discrepancy (1960). Since the design is high dimensional, we apply the radar to all pairs of possible dimensions. Among the rejected pairs we have, for example, the pair (14, 15), represented on Figure 3. In this example, since there exist values of ( ) N D a outside the circle of radius ks, the uniformity radar detects a non-uniform distribution for the 2-dimensional projections onto ( )

14 15
, X X . The largest deviation in uniformity is observed in direction a associated with an angle of approximately 135°, corresponding to the direction that is perpendicular to the visible diagonal alignments on the figure to the left. Here, we find ourselves faced with the wellknown defect of high-dimensional Halton sequences, which do not preserve a low discrepancy in projection [START_REF] Thiémard | Sur le calcul et la majoration de la discrépance à l'origine[END_REF], [START_REF] Morokoff | Quasi-random sequences and their discrepancies[END_REF]). Note, however, that the radar is not designed to systematically detect directions of alignment as we will see in the next example.

Example 2. Analysis of orthogonal arrays using 3D radar. Let us consider a 49 points linear orthogonal array of strength two in 3 dimensions [START_REF] Owen | Orthogonal arrays for computer experiments, integration and visualization[END_REF] and [START_REF]Analyse statistique et échantillonnage d'expériences simulées[END_REF]).

Like Latin hypercubes, these designs are often recommended for numeric experimental designs because of their appreciable properties in projection. Projected onto a surface, an orthogonal array of strength two always yields a regular grid of points. However, the nonredundancy of the 2-dimensional projections does not imply a good distribution of the points neither on the axes of the domain (here 7 packets of 7 points) nor in the 3D space, as we will see. The considered design is a linear orthogonal array of strength two. The way it is built implies that the points satisfy the relation x 1 +3x 2 +x 3 =0 (mod 7). Therefore, these points are located on 5 parallel planes. As a result, the distribution of the projections onto the axis perpendicular to these planes will not be satisfactory. This problem is never mentioned in computer experiments. However, it is well-known in the literature of experimental designs.

We applied the uniformity radar to the studied orthogonal array and represented the sphere of radius ks (the Kolmogorov-Smirnov statistic with a 95% confidence level), the points Note that all the points of the uniformity radar are inside the circle of radius ks and, as a result, the radar accepts the design although we can see that the points of the plane ( )

3 6 , X X
are not uniformly distributed. However, the discrepancy values are rather scattered with a low value for angle θ = 0°, and rather high values, for example, for angle 29 θ = °, which seems to correspond to the orthogonal direction of alignments. The idea to reject this type of design amounts, therefore, to defining a new statistic which introduces minimum and maximum discrepancies. In order to avoid scale problems, we suggest taking the ratio of these quantities.

We have:

[ ] [ ] 0,2 0,2 sup ( ) inf ( ) N N N D G D θ π θ π θ θ ∈ ∈ = (3)
This statistic should filter out designs which have a relative poor distribution in one direction.

This statistic has the advantage of being global. This increase the power of the corresponding statistical test. For a fixed value of N, the distribution of N G seems difficult to obtain other than by simulation. In the appendix, the table corresponding to N=1…100 is given. According to this table, the rejection threshold at level 95% for a design of 100 points is equal to 4.23. In 6 X 3 X example 3, the observed value of the statistic N G is equal to 6.07, which very clearly eliminates this design.

In the statistic table, we observe that the values stabilize as N increases, which leads us to believe that it is possible to specify its asymptotic behavior. This is in fact the case. The following demonstration is based on the theory of empirical processes.

Asymptotic distribution of the global statistic (3). Let

( )

1 ,..., N X X X =
, where X 1 ,..., X N are independent random variables from uniform distribution on [ ] 2

1,1

- . Denote ( ) R X θ as
the projection of X on axis L θ . Let ( )

, t t N A A Y θ be the following process: ] ] ( ) ( ) ( ) ( ) , , 1 1 1 1 1 ( ) ( ) 1 t t N N N i A i t t A i i Y R X PR X t X PX A N N θ θ θ -∞ = = = - < = - ∈ ∑ ∑ where A t = x ∈ R d ,R θ (x) < t { } .
Using these notations, the discrepancy of projections may be written as ( )

, sup t N N A t D Y θ θ ∈ = .
Therefore, , , sup ( ) sup sup inf ( ) inf sup

t t N N A t N N N A t D Y θ θ = = .
The limit distribution of the process ( )

, t t N A A
Y θ is obtained using the theory of empirical processes. With the central limit theorem, one obtains the convergence of the finite dimensional distributions toward the Gaussian process of same expectation and covariances.

Furthermore, the family (A t ) t is a VC-class (see …). Therefore, ( ) ( )

, t t t t loi N A A A A Y Y θ θ ⎯⎯→
where ( )

t t A A
Y θ is a centered Gaussian process with covariance:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 , 1 cov , ( ) , ( ) 1 ( ) . ( ) 1 ( ) . ( )
( ) . ( ) cov , ( ) , ( ) ( ) . ( ) cov , ( ) ( ). ( ). s t s t s t N A A i j N i N j A A s t s t A A Y Y P R X s R X t N P R X s P R X t N P R X t P R X s N P R X s P R X t Y Y P R X s R X t P R X s P R X t = < < - < < - < < + < < = < < - < < = ∩ - ∑ ∑ ∑ (4)
and

P(G N < y) N →+∞ ⎯ → ⎯ ⎯ P sup θ sup t Y A t θ inf θ sup t Y A t θ < y ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟
The probability ( )

s t P A A ∩
is interpreted as the surface of a polygon delimited by the domain and the straight lines s t A et A (see Figure 8). This surface can be analytically calculated, by noting (for example) that the polygon is a partition of triangles with a common corner to be selected from among the vertices of the polygon. The asymptotic distribution of N G can then be obtained by simulating a centered Gaussian field with covariance matrix defined by (4).

Figure 8. Interpretation of the probability ( )

P A A s t ∩ .
The procedure is semi-analytic because simulations are needed to perform the calculations. In theory, the global statistic may be extended to 3 dimensions. For a finite sample, the statistic 1 O=(0,0)

L μ t s t A s A -1 1 -1 L θ
table can be calculated with simulations. However, the asymptotic distribution is much more difficult to obtain.

DISCUSSION

Due to the complexity of the phenomena simulated by computer codes, which implies nonlinearities, distributing numerical experiments as uniformly as possible in the domain is preferable. In addition, this distribution should continue to be satisfactory when projecting onto subspaces, especially when the code depends only on a small number of factors or principal components. In this article we introduce a statistic based on 1-dimensional projections to test the uniformity for a design of experiment. In 2 and 3 dimensions, it provides a simple viewing tool called uniformity radar, which graphically tests the uniformity hypothesis by omni-directional scanning. Moreover, we introduced a global statistic in 2 dimensions to further identify unsatisfactory designs that had been accepted by the uniformity radar.

These tools were applied on usual SFD designs. On these case studies, some of the designs have very poor properties when projecting to subspaces, as the 15-dimensional sequence with low discrepancy in example 1 or the 3-dimensional orthogonal array in example 2. The uniformity radar was able to detect these defects. It succeeds when there is a rectangular shaped empty area in the domain, as in the aforementioned examples. In such cases, the distribution is unsatisfactory when projecting onto the rectangle's width. The uniformity radar can also detect the alignments of points, but may not identify them if the directions of alignments are well distributed such as with a factorial design (see example 2). This underscores the lack of power of the Kolmogorov-Smirnov test when the sample is generated from a continuous distribution supported by the union of small intervals regularly distributed.

In practice, this situation is not very detrimental because the SFD obtained by a deterministic process are often randomized or scrambled (see, for example, [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF].

Our uniformity radar may be adapted to other goodness-of-fit statistics, such as Cramér-Von Mises (see section 2), which corresponds to the discrepancy L 2 for a projection onto a coordinate axis. For instance, we repeated the examples 1 to 3 with the corresponding radar.

As expected, the conclusions are the same because the Kolmogorov-Smirnov and Cramér-Von Mises tests do not present any clear-cut difference in terms of power. Interestingly, the main difference is graphic. The curve of the radar defined with the Cramér-Von Mises statistic is smoother, which is due to the norm L 2 , and introduces sometimes large scale variations from one design to another, while these differences are attenuated by the norm L ∞

in the examples given here. For this reason, the L ∞ radar may be preferred because the conclusions are more apparent.

APPENDICES

Distribution of a linear combination for uniformly distributed independent variables.

The final result demonstrated by [START_REF] Ostrowski | Two Explicit Formulae for the Distribution Function of the Sums of n Uniformly Distributed Independent Variables[END_REF] often uses very technical methods.

However, [START_REF] Elias | Convolution of Uniform Distributions and Ruin Probability[END_REF] propose a simpler operational computational method that we apply to our special case.

Let

1 1 . ... d d Z X a X a X a = = + +
. This is, therefore, a sum of independent random variables of

uniform distribution on [ ] , , 1,..., i i a a i d - = .
As a result, it admits the density given by

[ ] [ ] ( ) 1 1 , , 1 1 1 1 2 d d d Z a a a a j j f a - - = ⎛ ⎞ = × * * ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ∏ L
, where * designates the convolution product. To calculate this product of convolutions, [START_REF] Elias | Convolution of Uniform Distributions and Ruin Probability[END_REF] suggest using translations. Let

( )( ) ( ) a T f x f x a =
+ , the translation of the function f from -a. Then we can write (with the convention 0 0, 0 (  )

α α = ∀ ≥ ): ( ) ( ) ( ) ( ) ( )( ) 0 0 , 0 0 0 1 ( ) ( ) ( ) (
1 1 0 0 1 1 ( ) 2 d d d Z a a a a j j f x T T x T T x a - + - + = ⎛ ⎞ = × - * * - ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ∏ L .
The proof follows from two basic lemmas.

Lemma 1. The translation commutes with the convolution.

( ) ( ) ( ) 1 2 1 2 1 2 1 2 ( ) ( ) , a x x x a x a x x T f f f T f T f f x x * = * = * ∀ ∈ Lemma 2. ( ) 1 , , ! ! 1 ! m n m n x x x n m x m n m n + + + + + * = ∀ ∈ ∀ ∈ + +
By applying lemma 1, we obtain:

( ) (
)( ) ( )

1 1 0 0 1 1 ( ) 2 d d d Z a a a a j j d fois f x T T T T x x a - - + + = ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ = × - - * * ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ∏ L L 14 4244 3 .
Then by applying lemma 2:

( )( ) ( )

1 1 2 2 1 1 1 ( ) ... 2 (1 ) ! d d d d Z a a a a a a j j x f x T T T T T T a d - + - - - = ⎛ ⎞ ⎛ ⎞ = × - - - ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ - ⎝ ⎠ ⎝ ⎠ ∏ . and { } ( ) 1 1 1 1 1 1,1 1 ( ) 2 ( 1 ) ! d d d d d Z d s a s a j j s x f x s s T T a d - + = ∈ - ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ = × ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ - ⎝ ⎠ ⎝ ⎠ ∑ ∏ L L .
Finally,

1 1 ( . ) 1 ( ) ( ) 2 ( 1 ) ! d d Z j s j x s a f x s a d ε - + = ⎛ ⎞ ⎛ ⎞ + = × ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ - ⎝ ⎠ ⎝ ⎠ ∑ ∏ where ( ) ( ) 1 1 ,..., [0,1] , , sup N U N z z z D z z F z z ∈ = - K
is the Kolmogorov-Smirnov statistic for the uniform distribution over [ ] For differentiability, the problem can be reduced to the existence of partial derivatives of U D .

• Let us assume that the supremum is attained for only one j z where { } 1,..., j N ∈ .

In this case, the partial derivatives related to i z , i j ≠ , exist and are equal to 0 because any local perturbation in i z would not change the supremum. Let us now prove the existence of

U j D z ∂ ∂ .
Let us assume, for example, that:

[0,1] ˆŝup ( ) ( ) 0 N N j j z F z z F z z ∈ -= -> . We choose ε such that ˆ, ( ) ( ) N j j N i i i j F z z F z z ε ∀ ≠ - > -+ et ˆ( ) N j j F z z ε -> . Then : ( ) ( ) 1 1
,..., ,..., ,...,

U U j N j D z z z D z z ε ε + = -et ( ) ( ) 1 1
,..., ,..., ,...,

U U j N j D z z z D z z ε ε - = + Therefore, ( ) 1 ,..., U N j D z z z ∂ ∂ exists in ( ) 1 ,..., N z z and is equal to -1. Similarly, if ( ) [0,1] ˆŝup ( ) N j N j z F z z z F z ∈ -= -
, we prove that Let us assume that the supremum is attained at j z and at least at another k z for

{ } , 1,..., j k N ∈
. Now, we are going to prove that the partial derivative in j z does not exist.

Assume, for example, that ˆ( ) 0 ,..., ,..., ,...,

U U j N N D z z z D z z ε ε + = -, but ( ) ( ) 1 1
,..., ,..., ,...,

U U j N N D z z z D z z ε - =
. Therefore, the right-hand derivative is different from the left-hand derivative.

• Otherwise, for all k z , where k j ≠ at which the supremum is attained,

j k z z ≠ . Thus ( ) ( ) 1 1
,..., ,..., ,...,

U U j N N D z z z D z z ε ε - = + but ( ) ( ) 1

Figure 1 .

 1 Figure 1. Left : Projections of points onto an axis L a . Right : The histogram of the projections.

Figure 2 .

 2 Figure 2. The distribution of the projections for a 3 dimensional cubic domain.

  2π], called 2D radar. By displaying the quality of distributions in all

Figure 3 .

 3 Figure 3. Left: Design projections of a Halton sequence onto the plane (X 14 ,X 15 ). Right: The 2D radar curve.

Figure 4 .

 4 Figure 4. Left. A 3-dimensional linear orthogonal array of strength two with 49 points. Right. Projections onto the plane (X 1 ,X 2 ).

Figure 5 .

 5 Figure 5. Left. The p-values in -log10 of the Kolmogorov-Smirnov test. Right. The values of the Kolmogorov Smirnov test with the radar's representation in pins.

Figure 6 .Figure 7 .

 67 Figure 6. Left to right. the distribution functions (after transformation) of the projected points onto L1, L2 and L3.

Figure 9 .

 9 Figure 9. Uniformity radar with the Cramér-Von Mises statistic for examples 1 and 3.

  value of the supremum, which is still attained at k z . , the right-hand derivative is therefore also different from the left-hand derivative.

ACKNOWLEDGMENTS

This work was conducted within the frame of the DICE (Deep Inside Computer Experiments) Consortium between ARMINES, Renault, EDF, IRSN, ONERA and TOTAL S.A. We also thank Chris Yukna for his help on editing.

Continuity and differentiability of the radar.

For the sake of simplicity, we give only the proofs in 2 dimensions. However, the result is also valid in 3 dimensions. Let

, the Kolmogorov-Smirnov statistic for the uniform distribution over [ ] 0,1 .

Proposition.

(i) U D is continuous.

(ii) U D admits partial derivatives with respect to all the , 1,...,


Severine Furst
email: severine.furst@univ-smb.fr

J Chéry

M Peyret

B Mohammadi

Tiltmeter data inversion to characterize a strain tensor source at depth: application to reservoir monitoring

Keywords: Strain tensor source, Reservoirs monitoring, Tiltmeters, Anisotropic deformation, Long-term drift, Volcano geodesy

Surface deformation measured by geodetic data is the sum of single strain sources deforming at depth. A combination of volume changes from several analytical models (e.g. a point source or dislocation along a plane) can be used to model the different sources. However, solving for the best fit of volume variations, dislocations, position and orientation parameters of all sources is a non-linear problem, and its solution is generally non-unique. This problem can be converted into a linear one by assimilating the sum of sources to a simplified model formed by three orthogonal planes of dislocations at fixed position and orientation. This strain source model is equivalent to having all neighbouring deformation sources contained in a small size volume. The determination of the strain tensor components can be performed by inverting geodetic data. Because of their high resolution, tiltmeters are well adapted to survey shallow deformation of volcanoes and geological reservoirs. However, they are known to display unknown long-term drift. We propose an approach to jointly estimate the temporal evolution of the strain source and time-dependent instrumental parameters. We verify the approach using synthetic data, giving confidence intervals for each component of the strain tensor. Finally, we link geological information to the internal deformation by interpreting the strain tensor as principal directions of deformation. This approach seems promising for the identification of fracture onset and fault reactivation in geothermal, hydrocarbon exploitations or volcanic systems.

Introduction

Surface ground deformation associated with resource extraction in subsurface or volcanic eruptions is commonly interpreted as change in pressure in deep reservoirs (e.g. oil, gas, salt, magma bodies). Analytical models including pressure sources [START_REF] Davis | Surface Deformation Due to Inflation of an Arbitrarily Oriented Triaxial Ellipsoidal Cavity in an Elastic Half-Space, With Reference to Kilauea Volcano, Hawaii[END_REF][START_REF] Mctigue | Elastic stress and deformation near a finite spherical magma body: Resolution of the point source paradox[END_REF] and planar dislocations [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF] are often used to model and characterize the behavior and dynamics of such reservoirs in purely elastic medium. For such models, the constitutive equations linearly link the surface deformation to the source at depth: for point source models described by [START_REF] Mctigue | Elastic stress and deformation near a finite spherical magma body: Resolution of the point source paradox[END_REF], the linearity lies in the source volume variation, while for Okada's models, it resides in the variation in fracture slip or opening. Others parameters, such as the position of the source, and for planar dislocation, azimuth, dip and geometry of the fracture, are non-linearly linked to the geodetic deformation signal.

These analytical models are suitable to represent the deep strain source at first approximation. They are widely used because of their low computational cost due to simple source geometries and homogeneous rheologies (e.g. [START_REF] Dvorak | Volcano Geodesy: the search for magma reservoirs and the formation of eruptive vents[END_REF][START_REF] Masterlark | Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions[END_REF]. Models of magma plumbing systems often include a magma reservoir (pressure source), a conduit (cylinder or ellipsoidal sources) or a dike (planar dislocation), which may or may not reach the surface [START_REF] Bato | Assimilation of Deformation Data for Eruption Forecasting: Potentiality Assessment Based on Synthetic Cases[END_REF][START_REF] Bonaccorso | Stromboli 2007 eruption: Deflation modeling to infer shallow-intermediate plumbing system[END_REF][START_REF] Montgomery-Brown | Geodetic evidence for en echelon dike emplacement and concurrent slow slip during the June 2007 intrusion and eruption at Klauea volcano, Hawaii[END_REF][START_REF] Palano | Ground deformation patterns at Mt. Etna from 1993 to 2000 from joint use of InSAR and GPS techniques[END_REF][START_REF] Segall | Volcano deformation and eruption forecasting[END_REF], while models of stimulated reservoir volume in hydraulic-fracture treatment consider the superposition of several nearby fractures assimilated to planar dislocations [START_REF] Astakhov | SPE 151017 A New Method of Characterizing the Stimulated Reservoir Volume Using Tiltmeter-Based Surface Microdeformation Measurements[END_REF][START_REF] Warpinski | Surface Tiltmeters : A Proven Tool for New Answers in Unconventional Reservoir Stimulations[END_REF][START_REF] Zhou | Tiltmeter Hydraulic Fracturing Mapping on a Cluster of Horizontal Wells in a Tight Gas Reservoir[END_REF].

The deformation of a complex fractured medium can be modelled as the sum of multiple fractures whose azimuth, dip, slip, opening and geometry are free parameters. It is fair to assume that the inversion of geodetic data considering many sources instead of a single one would improve the fit between observed and modelled ground deformation either at a specific time or over a long time period. However, this is made at the expense of the uniqueness of the solution of the selected mathematical formulation. For such ill-posed and non-linear formulation (when position and orientation is also sought), multiple combinations of parameters can explain the surface observations, and assumptions or a priori information are necessary to discriminate between solutions. In this paper, we propose to assimilate a given volume of a complex fractured medium (e.g. a fracture network of a reservoir) at a fixed position, to a uniformly deforming effective medium. We assume that the size of this unit effective medium is small relative to its depth, so that the surface deformation of the unit can be explained by a single strain tensor involving 6 independent components. In order to simulate the surface deformation associated to this anisotropic source, we use three orthogonal Okada's planes of fixed position and geometry. Similar study from [START_REF] Nikkhoo | Compound dislocation models (CDMs) for volcano deformation analyses[END_REF] has already combined three orthogonal Okada's planes to build a strain model, capturing all possible degrees of freedom and free of artefact singularities along the edges and below and above the vertices of a plane. [START_REF] Nikkhoo | Compound dislocation models (CDMs) for volcano deformation analyses[END_REF] considered arbitrary sizes and orientations of the dislocation planes in space, but only for a tensile dislocation. In our model, we consider 6 degrees of freedom: the 6 independent components of the strain tensor and their relation to the dislocation parameters along three orthogonal planes.

Because analytical solutions of ground displacement and tilt given by [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF] depend linearly on the dislocation values, we enforce the linearity of the problem. Once the components of the strain tensor are inferred by inversion of geodetic data, we can use geological, geophysical or geomechanical information to interpret the strain tensor in term of natural or induced fractures, dikes and other deformation features.

The geodetic techniques, such as Global Navigation Satellite System (GNSS), Interferometric Synthetic Aperture Radar (InSAR), tilt and levelling surveys, provide a time-space record of the surface deformation above various geologically active areas (volcanoes: e.g, [START_REF] Dzurisin | Volcano Deformation[END_REF][START_REF] Poland | Insights into shallow magmatic processes at Kīlauea Volcano, Hawaiʻi, from a multiyear continuous gravity time series[END_REF]; geological reservoirs: e.g, [START_REF] Maisons | Monitoring of slow ground deformation by satellite differential radar-interferometry . A reference case study[END_REF][START_REF] Vasco | Reservoir monitoring and characterization using satellite geodetic data: Interferometric synthetic aperture radar observations from the Krechba field[END_REF][START_REF] Verdon | Comparison of geomechanical deformation induced by megatonne-scale CO2 storage at Sleipner, Weyburn, and In Salah[END_REF]. Tiltmeters (long-base and borehole) are sensitive instruments commonly installed on volcanic systems to monitor structure deformation in real time along with geodetic network [START_REF] Montgomery-Brown | Geodetic evidence for en echelon dike emplacement and concurrent slow slip during the June 2007 intrusion and eruption at Klauea volcano, Hawaii[END_REF][START_REF] Peltier | Early detection of large eruptions at Piton de La Fournaise volcano (La Réunion Island): Contribution of a distant tiltmeter station[END_REF][START_REF] Gambino | Volcanic processes detected by tiltmeters: A review of experience on Sicilian volcanoes[END_REF][START_REF] Alpala | Monitoring remote volcanoes: The 2010-2012 unrest at Sotará volcano (Colombia)[END_REF][START_REF] Medina | Twenty years (1990-2010) of geodetic monitoring of Galeras volcano (Colombia) from continuous tilt measurements[END_REF]. Due to their high resolution (up to 1 nrad, [START_REF] Jahr | Monitoring fluid induced deformation of the earth's crust: A large scale experiment at the KTB location/Germany[END_REF], tilt sensors record a wide range of geological deformation signals at large distances from the source that are generally outside the resolution range of space geodetic techniques. For example, [START_REF] Kamo | A Preliminary Experiment on Automated Judgement of the Stages of Eruptive Activity Using Tiltmeter Records at Sakurajima, Japan[END_REF] and [START_REF] Nishi | Micro-tilt changes preceding summit explosions at Semeru volcano[END_REF] detected deformation of the order of a nanorad in volcano monitoring, with a maximum precursory tilt change of 10-200 nrad for periods ranging from 10 minutes to 7 hours prior to eruptions. Hydraulic fracturing of low permeable reservoirs [START_REF] Astakhov | SPE 151017 A New Method of Characterizing the Stimulated Reservoir Volume Using Tiltmeter-Based Surface Microdeformation Measurements[END_REF]Castillo et al., 1997;[START_REF] Zhou | Tiltmeter Hydraulic Fracturing Mapping on a Cluster of Horizontal Wells in a Tight Gas Reservoir[END_REF], phreatic eruption [START_REF] Honda | Precursory tilt changes associated with a phreatic eruption of the Hakone volcano and the corresponding source model[END_REF], long term injection of water into borehole [START_REF] Jahr | Monitoring fluid induced deformation of the earth's crust: A large scale experiment at the KTB location/Germany[END_REF][START_REF] Jahr | Deformation, seismicity, and fluids: Results of the 2004/2005 water injection experiment at the KTB/Germany[END_REF] or changes in pump rates in the vicinity of fluid-producing wells [START_REF] Chen | Field layout of a tiltmeter array to monitor microdeformation induced by pumping through a horizontal collector well[END_REF] produce small deformation signals detectable by tiltmeters. In addition to the source signal, tilt signal is altered by other strain sources due to tides, hydraulic loading, temperature effect or pressure gradient. In general, installing tiltmeters in deep boreholes leads to a decrease in noise amplitude [START_REF] Jahr | Monitoring fluid induced deformation of the earth's crust: A large scale experiment at the KTB location/Germany[END_REF]. Semi-empirical tidal models permit the time series to be corrected from tidal effects [START_REF] Van Camp | Tsoft: Graphical and interactive software for the analysis of time series and Earth tides[END_REF]. Nevertheless, the remaining noise in the signal needs to be filtered. Similarly to colored noise in GNSS data or Brownian noise in levelling data, the remaining noise in tilt time-series can be considered as a white noise accumulation over time, and thus can be considered as Brownian noise [START_REF] Boudin | A silica long base tiltmeter with high stability and resolution[END_REF]. Finally, a long-term drift plagues the tilt signal, making it difficult to use tiltmeters for determining slow deformation processes. Usually, the drift is removed from the tilt data before the inversion by applying a linear regression on the time series. However, the source can have a linear term in time. [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF] developed a methodology to invert time series of tilt data induced by linearly time-dependent strain variations of a source at a given position and depth. This method permits to automatically separate the instrumental drift from the source signal through a two-step approach: first, tilt data are inverted to determine a set of admissible parameters, and second, uniqueness of the solution, with regards to the mathematical problem setting, is enforced by assuming the independence between instrumental drift and source parameters. This provides a consistent estimate of instrumental and source parameters. Because this approach only requires the linearity of the forward model between the data and the strain source parameter, it is easily generalizable to the strain tensor model described above. In this paper, we extend the methodology to this model to jointly retrieve the components of the strain tensor and drift rates of tiltmeters for a given set of fixed model parameters. The paper opens with the description and advantages of the strain tensor model. Although our approach may be applied to any geodetic techniques, we illustrate its application to the continuous monitoring of small surface deformation using tiltmeters. We then describe the methodology including

(1) the tilt data parametrization, (2) a two-step approach (optimization and uniqueness enforcement) used to jointly retrieve source and instrumental parameters and (3) a resolution analysis performed to discuss the reliability of the optimal parameters. This two-step approach is then verified against a synthetic data set: the deformation induced by a 3D strain tensor of a single source at depth, lasting for 11 months. Finally, we discuss the potential of the developed methodology in fractures mapping in hydrocarbon exploitations and geothermal reservoirs or dike opening in volcanic systems.

Strain tensor model

Model description

A deforming fracture at depth can be mathematically assimilated to a planar dislocation. [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF] presents complete analytical expressions for the displacements and tilts caused by a strike-slip, dip-slip or tensile rectangular dislocation in a homogeneous, isotropic flat and elastic half-space. These expressions depend linearly on the dislocation (slip and tensile) parameters of the source and non-linearly on the position ( 𝑥, 𝑦, 𝑧 ) and orientation (azimuth and dip) of the source. Now, let us consider a complex fractured medium embedded in an elementary volume at a certain depth (Fig. 1a).

The surface deformation is the resultant of the slip (or opening) of all single fracture forming the complex unit.

Looking for the best fit of slip, opening, position and orientation parameters of all fractures by modelling the surface deformation is a non-linear problem and its solution is generally non-unique. Therefore, a linear form of the strain nucleus is the global strain tensor that sums up the contribution of all individual internal strain tensors. This model can be setup by the use of three orthogonal planes of fictitious dislocations at fixed position (the centroid of the considered volume) and orientation. This model is equivalent to a complex source model of the deformation at the surface (Figure 1b). A similar relation, linking the mean of the amount of displacement over the area of the fracture to the sum of seismic moments [START_REF] Aki | Generation and Propagation of G Waves from the Niigata Earthquake of June 16, 1964. Estimation of earthquake moment, released energy, and stress-strain drop from the G wave spectrum[END_REF] occurring during a period of time, was described by [START_REF] Brune | Seismic Source, Fault-plane Studies, and Tectonics[END_REF] and later used by [START_REF] Kostrov | Seismic moment and energy of earthquakes, and seismic flow of rock[END_REF].

To be representative of a complex medium, the size of the equivalent model must be small with respect to its depth.

The surface deformation can be inverted using this equivalent model, allowing to retrieve a unique solution for the effective deformation in the three orthogonal planes.

(𝑈 1 1 , 𝑈 2 1 , 𝑈 3 1 ), (𝑈 1 2 , 𝑈 2 2 , 𝑈 3
2 ) and (𝑈 1 3 , 𝑈 2 3 , 𝑈 3 3 ) correspond to strike, dip and tensile dislocations for plane 1, 2 and 3 respectively.

In the far field, this model is parametrized to produce the same deformation as the model in a).

Surface displacement and tilt result from changes of the strain tensor at depth, under the hypothesis of small deformations. The equivalent model of the fractured medium is the combination of three orthogonal planar dislocations with fixed position and orientation to take advantage of the linear dependency with slip and tensile parameters. The approximation of a single strain tensor stands as long as the source is relatively small compared to its depth. We estimate that the ratio of the cube edge (ℎ) over the depth of the source (𝑧 𝑠 ) should be smaller than 0.2. Above this value, the surface deformation pattern produced by the system of fractures significantly differs from the one given by the equivalent strain tensor. The analytical solution for the surface and internal displacements, strains and tilts due to shear and tensile on rectangular faults in a semi-infinite elastic half space has been described by [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. The source is a planar facet characterized by its length (𝐿), width (𝑊), location (𝑥 𝑠 , 𝑦 𝑠 , 𝑧 𝑠 ), orientation (dip, azimuth) and dislocation type (strike 𝑈 1 , dip 𝑈 2 or tensile 𝑈 3 ). The induced deformation is linear with dislocation 𝑈 𝑗 and non-linear with all other parameters.

Geodetic techniques such as GNSS, InSAR, levelling and tiltmeters measure the surface deformation induced by such a source. We propose analytical equations only for a tiltmeter network (gradient of vertical displacement).

Nevertheless, similar equations can be easily derived for other geodetic measurements. We will focus on the longterm monitoring of a reservoir using a network of tiltmeters. Using solutions given by [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF], the tilt signal 𝑑 𝑠 ⃗⃗⃗⃗ measured by 𝑁 instruments in both X and Y horizontal directions, can be expressed as the linear combinations of the dislocation values 𝑈 𝑖 with their associated deformation model vector 𝛼 𝑖 :

𝑑 𝑠 ⃗⃗⃗⃗ = ∑ 𝑈 𝑖 3 𝑖=1 𝛼 𝑖 (1)
where 𝑑 𝑠 ⃗⃗⃗⃗ and 𝛼 𝑖 are vectors of dimensions (2𝑁) and 𝛼 𝑖 correspond to the contribution of dislocation parameter i to the signal recorded by each instrument. In this study, we consider three orthogonal Okada's planes of fixed dimension (we consider the edge ℎ of the plane so that ℎ = 𝑊 = 𝐿) and orientation arbitrarily chosen in the 𝑥𝑦𝑧 reference system to build the elementary volume (Figure 1b). Okada's plane can be dipping at an arbitrary angle, and be orientated at any arbitrary azimuth angle but no plunge (angle between the upper edge of the plane and the free surface) is considered. Our model has a fixed orientation with only potential azimuth variations. In [START_REF] Nikkhoo | Compound dislocation models (CDMs) for volcano deformation analyses[END_REF], they address this lack of freedom by considering an extended rectangular dislocation with a full rotational degrees of freedom (i.e. characterizing the plunge).

The frame associated with each plane is represented in Figure 1, permitting the correspondence between dislocation types and the components of the strain tensor (Eq. 2). We associate to each plane a slip (strike or dip)

in one direction and a tensile deformation (𝑈 𝑖 𝑗 ) representing the strain deformation (𝜺). Due to the static equilibrium of the source within the medium, the strain tensor is symmetrical:

( 𝑈 3 1 𝑈 1 1 𝑈 1 3 𝑈 2 3 𝑈 3 2 𝑈 2 2 𝑈 1 2 𝑈 2 1 𝑈 3 3 ) = 1 𝑎 ( 𝜀 𝑥𝑥 𝜀 𝑥𝑦 𝜀 𝑥𝑧 𝜀 𝑦𝑥 𝜀 𝑦𝑦 𝜀 𝑦𝑧 𝜀 𝑧𝑥 𝜀 𝑧𝑦 𝜀 𝑧𝑧 ) = ( 𝑝 1 𝑝 2 𝑝 3 𝑝 2 𝑝 4 𝑝 5 𝑝 3 𝑝 5 𝑝 6 ) (2)
where 𝑈 𝑖 𝑗 are the dislocation types for all three planes: the subscript 1 to 3 stands for strike, dip and tensile for plane 1 to 3 (Figure 1). Hence, the tilt signal produced by a strain tensor source made of three Okada's planes is linear with respect to the 6 components of the strain tensor 𝑝 𝑖 (the strain tensor being symmetrical):

𝑑 𝑠 ⃗⃗⃗⃗ = ∑ 𝑝 𝑖 𝛼 𝑖 6 𝑖=1
(3)

The strain tensor source is hence described by 6 characteristic patterns of deformation 𝛼 𝑖 associated to each source component whose amplitude depends on the strength of the deformation, represented by the amount of dislocation 𝑝 𝑖 . In the next section, we describe these patterns of deformation as being the sensitivity associated to each component of the strain tensor.

Sensitivity map

The source tilt signal produced by a strain tensor given by Eq. 3 is the sum of the product between each strain parameter 𝑝 𝑖 and the associated deformation vector 𝛼 𝑖 . These vectors represent the known spatial variability of the forward model. They depend on the network distribution relative to the source location, dimension and orientation.

Thus, for a unit deformation of a strain source at a fixed depth and location, the values of ‖𝛼 ‖ show, for each component, the sensitivity of the tiltmeter to the deformation.

On Figure 2, we perform this general sensitivity analysis to the domain used for the synthetic application. To ease the interpretation, we choose to represent the sensitivity of each strain component

𝑠 𝑖 = ‖𝛼 ⃗⃗ 𝑖 ‖ ‖𝛼 ⃗⃗ ‖ 𝑚𝑎𝑥 ( 𝑖 = 1,6 )
normalized by the maximum value of ‖𝛼 ‖ for all 𝑖. The scale of sensitivity to the deformation ranges between 0 to 1, from insensitive to totally sensitive areas (Figure 2). Mean values of sensitivities are estimated for each component, over the 10x10 km domain. The 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 components present higher sensitivities with a mean value of 𝑠 𝜀 𝑥𝑧 ̅̅̅̅̅ = 𝑠 𝜀 𝑦𝑧 ̅̅̅̅̅ = 0.17 and 𝑠 𝜀 𝑧𝑧 ̅̅̅̅̅ = 0.32 compared to 𝜀 𝑥𝑥 , 𝜀 𝑥𝑦 and 𝜀 𝑦𝑦 components (𝑠 𝜀 𝑥𝑥 ̅̅̅̅̅ = 𝑠 𝜀 𝑦𝑦 ̅̅̅̅̅ = 0.11 and 𝑠 𝜀 𝑥𝑦 ̅̅̅̅̅ = 0.10). Indeed, tiltmeters measure the horizontal variations of the vertical changes of the ground (𝑡 = -∇ ⃗ ⃗ 𝑢 𝑧 ). This suggests that 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 components should always be better resolved through the optimization process regardless of the network design.

Parameters identification

We estimate both source and instrument parameters from tilt signals using a two-step approach, employing 1) a least-square inversion and 2) a uniqueness enforcement.

Model Parametrization

The input data is measured using an array of 𝑁 tiltmeters composed by 𝑇 monthly-discrete measurements associated to the variation of the strain tensor components at depth. Following [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF], we assume that the tilt signal 𝑑 𝑜 ⃗⃗⃗⃗ (𝑡) (2𝑁×𝑇 ) is the sum of the signal produced by the source 𝑑 𝑠 ⃗⃗⃗⃗ (𝑡) (2𝑁×𝑇 ), a time-linear instrumental drift 𝑑 𝑑 ⃗⃗⃗⃗ (𝑡) = 𝑎 𝑡 (2𝑁×𝑇) (𝑎 is the drift rate vector of dimension 2𝑁) and coloured noise 𝑐𝑛 ⃗⃗⃗⃗ (𝑡) (2𝑁×𝑇) associated to the tiltmeters (Eq. 3). We assume that the tidal signal has been removed from the tilt signal.

This synthetic signal is then used in the inversion process as the observation input. Seeking for source (𝑝 𝑖 ) and instrumental (𝑎 ) parameters, we parametrize the modelled tilt as the sum of a source signal evolving in time and a time-linear component standing for the drift component.

Step 1: Global optimization framework

Although the problem is presented with linear governing equations for the strain tensor model, we present the inverse problem in its nonlinear form and use a global optimization framework. This enables the inversion platform to easily account for complex forward models, like poroelasticity and elastoplasticity, or extend the identification strategy to source locations.

This study considers 2𝑁 instrumental parameters and 6𝑇 source parameters to be retrieved using time series of tilt data including 2𝑁 • 𝑇 observations. The optimization parameters are initialized with no a priori information, providing a first model prediction to be compared with the observations. During this stage of initialization, we also estimate the deformation model vector 𝛼 𝑖 of a unit source for both components of each tiltmeter (using [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. Then, the modelled signal estimated at each iteration is simply the product of the parameter 𝑝 𝑖 (𝑡) by the deformation vector 𝛼 𝑖 , allowing for a significant gain in computational time. The comparison between observed and modelled data is made using the weighted squared error as cost function, integrated over time following the trapezoidal rule (Eq. 5 from [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF]. The functional needs to converge below one (ideally to 0) for the optimization to be complete, i.e. reaching a set of admissible parameters explaining the observations within the data uncertainties. To reach this minimum, our inversion process is based on a recursive algorithm defining a multi-criteria global optimization, varying not only the parameters of the model, but also the initial guesses [START_REF] Ivorra | Design of code division multiple access filters based on sampled fiber Bragg grating by using global optimization algorithms[END_REF][START_REF] Mohammadi | Applied Shape Optimization for fluids (2nd Editio)[END_REF]. This ensures that a given set of optimal model parameters achieves the global minimum of the cost function.

At the end of the optimization, the parametrization of tilt data leads to a system of linear equations whose solutions are infinite combinations of these parameters producing the same signal [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF]:

𝑑 𝑎 ⃗⃗⃗⃗ (𝑡) = ∑ 𝛼 𝑖 𝑝 𝑖𝑎 (𝑡) 6 𝑖=1 + 𝑎 𝑎 ⃗⃗⃗⃗ 𝑡 = ∑ 𝛼 𝑖 𝑝 𝑖 * (𝑡) 6 𝑖=1 + 𝑎 * ⃗⃗⃗⃗ 𝑡 (4)
The subscript a represents any of all admissible sets including the best solution we reach with respect to the data uncertainties (expressed by the exponent * ), hereafter referred to as the target solution.

Step 2: Uniqueness enforcement

This set of admissible source and instrumental parameters (𝑝 𝑖𝑎 (𝑡) and 𝑎 𝑎 ⃗⃗⃗⃗ ) is not unique since an infinite number of combinations of these parameters leads to the same value of the functional [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF]. This implies that the solution of the optimization process can notably diverge from the target solution. The latter can be recovered from the admissible solution (Eq. 4) as,

𝑎 * ⃗⃗⃗⃗ = 𝑎 𝑎 ⃗⃗⃗⃗ -∑ 𝛼 𝑖 𝑅 𝑖 6 𝑖=1 , (5a)
with,

𝑅 𝑖 = p 𝑖 * -𝑝 𝑖𝑎 𝑡 (5b)
where 𝑅 𝑖 are the correction coefficients associated with each component of the strain tensor. Following Furst et al.

(2019), we assume that the position of the instrument relatively to the source is not correlated to the drift rate components. Hence, the distributions of α ⃗ ⃗ 𝑖 and 𝑎 must be statistically independent for each component of the strain tensor, that is 𝑐𝑜𝑣(𝑎 * ⃗⃗⃗⃗ , 𝛼 𝑖 ) = 0. Introducing Eq. 5a in the latter equation leads to

𝑐𝑜𝑣(𝑎 𝑎 ⃗⃗⃗⃗ , α ⃗ ⃗ 𝑖 ) -∑ 𝑅 𝑘 6 𝑘=1 𝑐𝑜𝑣(𝛼 𝑘 ⃗⃗⃗⃗ , α ⃗ ⃗ i ) = 0 (6a)
that we rewrite in matrix form,

𝐶 -𝑅 ⃗ 𝑩 = 0 (6b)
The coefficients 𝑅 ⃗ can be inferred from Eq. 6b if the square matrix 𝑩 = 𝑐𝑜𝑣(𝛼 𝑘 ⃗⃗⃗⃗ , 𝛼 𝑖 ⃗⃗⃗ ) is invertible (det (𝑩) must be different from 0). We normalized matrix 𝑩 by its Frobenius norm to evaluate the determinant of the matrix. In our study, 𝑩 is invertible but this property should be checked for new instrumental distribution. This allows us to determine the coefficients 𝑅 ⃗ such as:

𝑅 ⃗ = 𝑩 -1 𝐶 (7)
The uniqueness enforcement produces six independent values 𝑅 𝑘 which depend on the covariance matrix between admissible drift and deformation model vector, but also involve the covariance between the components of the deformation model vector. In [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF], 𝑅 coefficients corresponds to the linear regression between the drift rates and the components of the deformation model vector. Here, 𝑅 𝑘 also involves the covariance matrix (𝑩)

between the components of the deformation model vector, and therefore they depend on all 𝛼 𝑖 associated with each component of the strain tensor. This implies that the distribution of 𝑎 𝑎 ⃗⃗⃗⃗ as a function of 𝛼 𝑖 , would not present a clear linear trend.

Incorporating 𝑅 𝑘 into Eq. 5a and 5b provides us with a new admissible set of parameters (target solution) which should be closer to the exact solution without altering the residuals between observed and modeled data.

Nevertheless, the target solution does not necessarily match the exact solution due to data uncertainties that are not captured by the model.

Parameters resolution analysis

We account for the data uncertainties in the inversion process by using the weighted Euclidian norm in 𝐽. Data uncertainties propagate through the inversion process to the parameters so that our target solution is the best guess with respect to the resolution of the data. For linear problems, we can directly estimate the resolution matrix associated with the least-square inversion process. Nevertheless, we propose a generic approach to estimate the parameter resolution, regardless of the linearity of the problem. We attempt to estimate ranges of confidence for each component of the strain tensor by adding small perturbations on the data. The mathematical development to obtain this parameter resolution analysis is detailed in Appendix A. This leads to two confidence intervals of the parameters: an instantaneous resolution 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ that only depends on data residuals, and a long-term resolution 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ linked to a priori range of perturbations on the drift rates (see Appendix A for description of drift rates perturbations choice).

Algorithm

The following guideline summarizes the implementation of the method:

1. Tilt time series. The observed data are corrected from the tidal effect or any other known effects that may obscure the effect we are interested in.

2. Time intervals. Data are discretized in a user-defined number of epochs. This also defines the number of different source parameters that are sought.

Initialization.

User defines the number of strain sources. For each one, a certain number of model parameters are fixed (𝑥 𝑠 , 𝑦 𝑠 , 𝑧 𝑠 , 𝑊, 𝐿). Research intervals are given for free parameters (𝑝 𝑖 (𝑡)) of each source and for instrumental parameters (𝑎 ). The source model vector 𝛼 𝑖 is calculated during the initialization.

4.

Step 1: Inversion.

5.

Step 2: Uniqueness enforcement.

Correction coefficients.

Estimation of the coefficient 𝑅 𝑘 defined by Eq. 7.

5.2. Target solution. Admissible solution is corrected using Eq. 5a and 5b producing the target solution.

In the next section, we apply this algorithm to a synthetic case of a strain source deforming at depth.

Application to synthetic data

Model parametrization

The two-step optimization approach generalized from [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF] provides unique and constrained parameters in linearly dependent forward models. We create the synthetic dataset from a single strain source (three Okada's planes) deforming at depth. We assume a continuous deformation during 11 months split into 12 1-month intervals (initial step of the model is considered at 𝑡 = 0, final step at 𝑡 = 11 months), so we have 12 steps of the strain model to be characterized. The strain source is expressed in the reference frame made of three orthogonal planes of 100 m side embedded in an elastic medium at 1500 m depth, at the center of a 10x10 km domain. We choose random and independent time steps for each component of the strain tensor (black curves in Figure 3a-b), producing a complex surface signal. In Figure 3c, we represent one state (𝑡 = 7 months) of the vertical deformation (color scale) and synthetic tilts (black arrows) associated with the values of strain tensor components. The total volume variation at depth (Figure 3b) is estimated using the trace of the strain tensor,

Δ𝑉 𝑖𝑛𝑗𝑒𝑐𝑡𝑖𝑜𝑛 = (𝑝 1 + 𝑝 4 + 𝑝 6 )ℎ 2 (8)
Besides the signal of the source, a time-dependent linear drift and a Brownian noise are included in the synthetic tilt signal. The drift components are randomly chosen using a uniform probability density function. The range of probability for drift rates is set to ±2.4 µrad/yr, corresponding to low drifting tiltmeters [START_REF] Chawah | A simple pendulum borehole tiltmeter based on a triaxial opticalfibre displacement sensor[END_REF][START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF]. The Brownian noise is arbitrarily generated by accumulating random values within a standard deviation, depending on the time of the experiment. The initial standard deviation fixed using the short-term resolution of tiltmeters 𝜎 𝑠ℎ𝑜𝑟𝑡 =5 nrad grows towards a maximum standard deviation 𝜎 𝑚𝑎𝑥 =180 nrad at the end of the experiment. In our synthetic measurements, the covariance matrix is diagonal and uniformly set to the final value of standard deviation 𝜎 𝑚𝑎𝑥 for all data. The noise signal is significantly smaller than the drift signal but it can affect both the instantaneous and long term trend in the tilt signal. Indeed, the linear component of Brownian noise can be superimposed to the linear instrumental drift. As a result, the optimization estimates the sum of both linear signals and interprets the value as the linear drift only.

The synthetic deformation induced by the strain source is monitored using an array of 50 tiltmeters randomly distributed (Figure 3c, [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF]). From Figure 2, we can see that this network is not optimized to capture at best all 6 components of the strain tensor since less than half the instruments are located in sensitive areas.

Continuous tilt data are downsampled to monthly time-interval measurements to decrease the number of parameters to infer, leading to 72 strain parameters and 100 drift parameters for 1200 tilt observations. In Figure 3c, we represent the vertical displacement of the surface at 𝑡 =7 months. The components of the strain tensor are indicated by red dots in Figure 3a, and the associated total volume variation is of 66 667 cubic meters (red dot in

Results

To retrieve the 6 parameters of the strain tensor, we invert the synthetic tilt dataset using the two-step optimization.

Once the inversion of tilt data has converged, we obtain a set of optimal parameters including components of strain tensor and drift rates given by the tilt residuals. This admissible solution needs to be corrected using the correction coefficient to restore uniqueness. We hereafter describe the results of the two-step optimization in terms of tilt residuals, drift rates parameters and strain tensor components. We complement the strain tensor components with instantaneous resolutions and long-term ranges of confidence as previously described.

Residual tilt signal

The best fit model produces a dataset giving the lowest residual between synthetic and modelled data. The data residuals are determined using the norm of the difference between modelled and observed tilt vector for each instrument and can be estimated at 3 different scales: 1) at data point scale, 2) at time-step scale and 3) at experiment scale. Figure 4a plots these residuals (crosses) as a function of time-steps. The inversion process provides a fairly homogeneous tilt residual over time for the whole set of tiltmeters (Figure 4a). Figure 4b displays the spatial distribution of data residuals at time 𝑡 = 7 months (blue crosses on Figure 4a). By averaging the residuals of all data points at a given time, we obtain a global residual for the considered time interval (red crosses on Figure 4a). Finally, we integrate the time-interval residuals over the entire experiment to get a mean global residual (green line on Figure 4a). Time-interval residuals range from 0.007 to 0.024 µrad with an average experiment residual of 0.010 µrad (green line on Figure 4a). Compared to the tilt uncertainty of 0.18 µrad, the residuals are significantly smaller. This is probably due to the parametrization of the tilt data, as being the sum of source and time-linear signals. Hence, the optimization tends to estimate all linear trends which may be present in the tilt signal, including drift rates and any noise that would be linear in time.

Drift rates

This admissible solution given by the first step of the optimization process belongs to a family of solutions described by Eq. 4. In Figure 5a and b, we plot the elements of 𝑎 (𝑥 and 𝑦 components) as a function of 𝛼 𝑧𝑧 ⃗⃗⃗⃗⃗⃗ associated to the 𝜀 𝑧𝑧 -component of the strain tensor (Figure 5b being an enlargement of the grey area of Figure 5a). We choose to represent and describe only one component of the strain tensor where blue dots stand for admissible drift rates after the inversion, red inverted triangles for the corrected values and black crosses for the true synthetic values.

Before the uniqueness enforcement, the distribution of the drift rates was between bounds (given by tiltmeter manufacturers) resulting from the first step of the optimization. One can identify a trend between elements of 𝛼 and elements of 𝑎 inducing some correlation between the drift rates and the position of the instrument relatively to the source (blue line). In [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF], the correction is based on removing this correlation. In this study, the same process (that is, without taking into account the cross correlation between drift rates and model deformation vector) produces the distribution of drift rates (green diamonds), which is far from the true distribution of drift rates (black crosses). Here we use a coupled coefficient described by Eq. 7 to correct the optimal drift rates leading to values of drift rates (red inverted triangles) almost perfectly retrieved (black crosses).

Strain tensor components and volume variations

Along with the drift rates, the inversion produces admissible variations of strain tensor components 𝑝 𝑖𝑎 (blue line in Figure 6) depending on the tilt residuals. The shape of each component variation is preserved compared with the target history (black line in Figure 6), but a strong linear trend persists due to non-uniqueness of the solution.

We interpret this linear trend as the consequence of incorrect estimation of the instrumental drifts. Admissible solutions 𝑝 𝑖𝑎 are corrected using the coupled coefficients 𝑅 𝑖 in Eq. 5b to obtain the target solutions 𝑝 * (red line in Figure 6). For all components, the uniqueness enforcement largely improves the solution. For instance, a slight linear trend remains in 𝜀 𝑦𝑦 while 𝜀 𝑥𝑥 and 𝜀 𝑥𝑦 are almost perfectly restored.

While we demonstrate the efficiency of our strategy on synthetic data, the target solution is not available when using real case datasets. As a result, we complement these target solutions 𝑝 * with two resolution analyses: one considering each time-step separately 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ and a second one involving the entire time series 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ . The instantaneous resolution is determined using the tilt residuals for each component of a tiltmeter at a given time. The residuals associated to each tilt component of the experiment range between 0 and 0.08 µrad (Figure 4a), which is largely lower than the Brownian noise generated in the synthetic. We assume that the level of resolution of the tilt in the inversion depends on the parametrization of the tilt data. Therefore, the instantaneous resolution (Eq. A.3, blue dots in graphs below each parameter variations on Figure 6) is 2 orders of magnitude lower than the parameters.

In Figure 6, we only represent the positive instantaneous resolution (the negative bound being the symmetric). One can notice that the instantaneous resolution is higher for 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 than for 𝜀 𝑥𝑥 , 𝜀 𝑥𝑦 and 𝜀 𝑦𝑦 . This difference lies in the matrix 𝑨, which contains the deformation model vectors associated to each instrument for all of the strain tensor components. 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 components show high sensitivity compared to the three other components (Figure 2), resulting in almost constant instantaneous resolution for 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 and more variable instantaneous resolution for 𝜀 𝑥𝑥 , 𝜀 𝑥𝑦 and 𝜀 𝑦𝑦 components.

The long-term resolution is built using a constraint on the drift rates (user defined, see Appendix A) and increases with time. As part of the optimization parameters, the drift rate distribution is defined by its dispersion around a mean value of 0 (i.e. its standard deviation) and by the gradient of the functional for each drift parameter (Eq. A.4).

Hence, 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ depends on inaccuracy in the minimization of the functional and on the quality of the tiltmeters. Indeed, highly drifting instruments would present a high standard deviation, increasing the range of sensitivity.

Perturbations of the drift rates impact all components of the strain tensor in the same way but 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ also depends on the tilt network through matrix 𝑨 similarly to 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ . Therefore, the range of sensitivity is different for each component as shown in Figure 6 (green area). In our study, the resolution matrix is close to the identity, meaning that the tilt network is not blind to any component of the deformation at depth. All confidence intervals include the target solution but as expected from the sensitivity map. 𝜀 𝑥𝑥 , 𝜀 𝑥𝑦 and 𝜀 𝑦𝑦 components present a wide range of resolution contrarily to that of 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 . The values of at the end of the experiment, 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ for the 6 components of the strain tensor are given in Table 1. Optimizing the spatial distribution of the tilt network improves the target solution and the range of resolution, but 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 components will always present a better range of resolution than the other components.

𝑃 𝐿𝑇 (𝑡 = 𝑡 11 )

(m) 𝑝 1 𝑝 2 𝑝 3 𝑝 4 𝑝 5 𝑝 6
±9.3 ±5.9 ±1.7 ±4.4 ±1.5 ±0.7

Tab. 1 Long-term resolution 𝑃 𝐿𝑇 associated with each component of the strain tensor at the end of the experiment (𝑡 = 𝑡 11 ) From the components of the strain tensor, we can derive the history of volume variations defined by Eq. 8 and represented in Figure 7a. Instantaneous and long-term resolutions can also be estimated using the same relation as

Eq. 8 and involving 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ and 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ associated to 𝜀 𝑥𝑥 , 𝜀 𝑦𝑦 and 𝜀 𝑧𝑧 components. Figure 7a shows the target volume as a black line, the admissible volume as a blue line and the corrected volume as a red line within a range of sensitivity in green. Because the long-term resolution of each diagonal component is added, the range of sensitivity is wide.

Discussion

In this study, we present a concept, based on standard Okada's model, for describing deep complex fracture sources. Under the assumption that the source size is small compared to its depth, we state that a strain tensor model captures the sum of all tensile, dip and strike-slip fracture motion. Therefore, the surface motion associated to the strain tensor is equivalent to the one produced by the fracture network. This strain tensor is an effective image, characterized by six parameters, which intrinsically reflects the fracture network constituting the considered volume. Associated to the strain components, patterns of deformation describe the sensitivity of tiltmeters to the deformation. The sensitivity analysis highlights that all components of the strain tensor cannot be solved equally well. This is probably due to the limited capability of tilt data (2-D) to characterize all components of the strain tensor (3-D entity). As tiltmeters are more sensitive to the vertical displacement of the source, the strain components involving horizontal motion are poorly determined. To overcome the problem, the use of other geodetic measurements sensitive to horizontal displacement of the surface, such as 3-D GNSS data or strainmeters may improve the estimation and resolution of the 𝜺 𝒙𝒙 , 𝜺 𝒙𝒚 and 𝜺 𝒚𝒚 components. Such a network design is conceivable for oil and gas exploitations but also for some volcanoes where the monitoring involves a wide range of geodetic network.

A major advantage of our approach is to use a well-posed inverse problem. Indeed, the two-step methodology developed in the paper allows to retrieve the unique solution of the strain tensor using surface tilt data (or any other geodetic measurements). However, the issue remains to associate the parameters of the global strain tensor with the individual tensile and slip motions of a fracture set. Without any complementary information, the strain tensor can only be analyzed in terms of principal values and directions of the volume. To determine orientation and internal slip motions associated with fracture families constituting the volume, the use of external information is necessary like geological studies, seismic events (in volcanic systems) and microseismicity (in geothermal, mining and oil and gas exploitations). This a priori information permits to characterize active fracture or fault planes in terms of position and geometry, and the question is whether or not their dynamic can be obtained using the strain tensor.

In oil and gas exploitations, the injection of fluid in a tight reservoir produces a nanometric scale deformation along with a microseismic signal typically aligned to form a fracture. With fracture positions and orientations determined through microseismicity (e.g. [START_REF] Wessels | Identifying faults and fractures in unconventional reservoirs through microseismic monitoring[END_REF], we can infer the type (tensile, strike or dip) and amplitude of slip motion. Figure 8 illustrates schematic examples of 2-D fractured domains characterized by microseismicity (blue and red dots) with strain tensors (black and with arrows) resulting from the optimization approach. For each case, the strain tensor can be interpreted using a priori knowledge about position and orientation of the fractures, leading to slip motions of the fractures (tensile for Figure 8a and strike slip for Figure 8b). To do so, we solve a system of linear equations. In 2-D the model considers 1 Okada's plane having an infinite width.

This system includes 3 equations (number of independent strain parameters) and two unknowns defined by the types of slip associated to each fracture family. When only one fracture family is considered (Figure 8a-b), the solution is overdetermined. However, this overdetermined system should have solutions because of fracture configuration: tensile motion can only occur in the orthogonal direction of the fracture plane. Hence, the associated system should contain only two equations (𝜺 𝒙𝒙 or 𝜺 𝒚𝒚 = 𝟎) for two unknowns. On the other hand, if the system truly includes three independent equations, this implies that several fracture families are involved, leading to an under-determined problem. In such a case, more assumptions (e. g. fixing type of slip motion) need to be made to discriminate between the solutions. Similarly, in 3-D (3 Okada's planes) the system of equations is defined by the six independent parameters of the strain tensor, while the equivalent fracture model has three unknowns (tensile, strike and dip-slip). Hence, types and amplitudes of slip motions can be fully determined if two fracture families have been identified. The general case involving all slips of a fracture family need to be defined as an inverse problem with additional constraints such as the strength and friction of the fractures, in order to reduce its intrinsic non-uniqueness. This approach seems to be feasible by using the methods developed by [START_REF] Angelier | Tectonic analysis of fault slip data sets[END_REF][START_REF] Angelier | The tectonic development of the hellenic arc and the sea of crete: A synthesis[END_REF] and by [START_REF] Etchecopar | Etude des états de contrainte en tectonique cassante et simulations de déformations plastiques: approche mathématique[END_REF] for fault slip and stress inversion. This model of strain tensor can be used to characterize complex fracture networks occurring in hydraulic fracturing in mining, geothermal or oil and gas extraction. Indeed, the injection of material at depth along kilometric scale horizontal wells, generates fractures in addition to natural fracking, creating networks with various fracture orientations. Therefore, during the exploitation of the reservoir, these fractures may behave differently depending on the applied stress. Thanks to the full strain tensor model, we can estimate from the surface data the state of the deformation at depth. The strain tensor would allow the identification of preferential deep movements, leading to families of fractures or faults that have played to produce the surface signal. In the specific case of oil exploitation, this would improve monitoring of exploitations by identifying potential pre-existing fractures or faults at depth, that are reactivated during injection or extraction. Furthermore, modelling the deformation of a superficial source requires to consider a collection of strain tensors to respect the far field equivalence. Hence, an efficient discretization of the strain sources representing stimulated volumes along wells would allow to estimate variations of production but also temporal and spatial strain orientations. Finally, the use of the strain tensor model in volcanic systems would provide information on deformation of the magma chamber or magma movements in dikes or sills, and potentially identify anisotropic deformation of the various sources.

Conclusion

In this study, we address the issue of monitoring and characterizing complex fracture reservoirs. Such reservoirs are common in mining, oil or gas exploitations and need a careful monitoring scheme to optimize the production and to follow the induced deformation. Even if the individual characteristics of such complex sources cannot be uniquely determined from surface deformation using a classical analytical source, the deformation of a complex fracture network still results in a unique surface deformation. When the source considered is significantly small with respect to its depth, we can express the deep deformation in a new basis to obtain the global effective deformation of the medium. The choice made here for this basis is three mutually orthogonal Okada's dislocations.

The resulting strain tensor of the deep source produces the same far field signal as the true complex medium.

We develop a source model of fixed position, dimensions and orientation that can represent various amplitude of shear and tensile deformations at depth. By doing so, we reduce a complex non-linear problem to a linear problem with a unique solution between the surface deformation and the source signal. We adapt the two-step methodology developed in [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF] to infer the source and instrumental parameters from tilt signal at the surface. The linearity of the strain tensor model significantly reduces the computational cost when inverting tilt data (or any other geodetic measurements).

The synthetic study shows the strong dependence of the instrument location with respect to the source and the sensitivity of the instrument. In particular, tiltmeters measure the gradient of the vertical displacement and are therefore more sensitive to the z-component of the strain tensor than to purely horizontal components. The confidence intervals illustrate this dependency with a better resolution for 𝜀 𝑥𝑧 , 𝜀 𝑦𝑧 and 𝜀 𝑧𝑧 components than for 𝜀 𝑥𝑥 -, 𝜀 𝑥𝑦 -and 𝜀 𝑦𝑦 -components. Complementing the tiltmeter network with additional geodetic data, such as GNSS or strainmeters network, could bring a better resolution of the parameters in the horizontal components.

Implemented on real datasets, the unique solution from the tilt data inversion can be interpreted using the principal directions of the tensor, leading to main deformations of the complex fracture network. Using available geology of the subsurface structures, we could identify families of fractures producing the surface deformation and report any reactivation of faults or fractures.

Both resolutions depend on matrix 𝑨 which means that they are influenced by the position of the tiltmeter network relatively to the source and by the physics of the source itself (i.e. the type of model used, here the strain tensor).

Instruments placed in sensitive areas should improve determination of the parameter and consequently narrow the confidence interval. The difference between the two resolutions lies in the perturbation component: the short-term resolution 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ depends on the perturbation of tilt data 𝑆 while the long-term resolution 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ is influenced by perturbations on the drift rates 𝐷 ⃗ ⃗ and increases with time.

In this study, we set the perturbations on the data 𝑆 using the residual value for each component of the tilt signal, defined by the difference between optimal and observed signal. This residual is linked to both the coloured noise given by the inversion (at the optimum). Small values of gradients mean better resolved parameter 𝑎 𝑖 . We can estimate the dispersion associated to the optimal drift rate distribution through its associated standard deviation 𝜎 𝑑 . We choose to estimate the 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ using a 2-𝜎 𝑑 uncertainty. However, this value being homogeneous for all drift parameters, we weight this standard deviation according to the gradient of the functional for the drift parameters such as , 𝜺 𝟏 and 𝜺 𝟐 being user defined. By doing so, we want to increase the uncertainty on less resolved parameters. Returning to the interpretation of parameter resolutions, the short-term resolution 𝑷 𝑰𝑻 ⃗⃗⃗⃗⃗⃗ , depends on the capacity of the inversion to converge towards an optimal and unique solution. Because tiltmeters present an excellent instantaneous resolution (up to 5 nrad), 𝑷 𝑰𝑻 ⃗⃗⃗⃗⃗⃗ at a given time would be well retrieved. On the contrary, the long-term resolution 𝑷 𝑳𝑻 ⃗⃗⃗⃗⃗⃗⃗ depends on errors on the drift rate estimation and increases with time. For highly drifting tiltmeters, the high value of variance [START_REF] Furst | Joint estimation of tiltmeters drift and volume variation during reservoir monitoring[END_REF] implies large ranges of confidence.

Fig. 1

 1 Fig. 1 a) Sketch of a complex fractured medium (e.g. fracture network in reservoirs characterization) with principal directions of deformation. b) Model of the effective unit built with three planar dislocations of edge ℎ and associated frame such that

Fig. 2

 2 Fig. 2 Normalized sensitivity distributions for the 6 components of the strain tensor. The sensitivity is estimated using the norm of the vector 𝛼 for each component normalized by the maximum value of the norm. The color code indicates low sensitive areas in white and highly sensitive areas in black. The locations of the tiltmeters in the synthetic network (see section 4) are represented by black dots.

4. 1 . 4 . 4 .

 144 Forward model. A set of modelled data is obtained using the parameters associated to the strain tensor model. 4.2. Functional. Computation of the cost function 𝐽. 4.3. Optimization. While 𝐽 ≤ 𝜖 (𝜖 is user defined), the solver defines new values of free parameters for the forward model (4.1). Admissible solution. Inversion gives a set of admissible solutions for 𝑝 𝑖 * (𝑡) and 𝑎 * ⃗⃗⃗⃗ .

Fig

  Fig 3b).

Fig. 3 a

 3 Fig. 3 a) Evolution of the six components of the strain tensor (𝜀 𝑥𝑥 , 𝜀 𝑥𝑦 , 𝜀 𝑥𝑧 , 𝜀 𝑦𝑦 , 𝜀 𝑦𝑧 , 𝜀 𝑧𝑧 ) for an anisotropic (black curve) strain source. b) Evolution of the total source volume variations (Eq. 8) for the 11 months of experiment. c) Vertical deformation (color scale) and synthetic tilts (black arrows) produced by an anisotropic strain source observed at the surface at 𝑡 =7 months (red dot on a and b).

Fig. 4 a

 4 Fig. 4 a) Time distribution of the norm of tilt residuals for all data points (crosses). Red crosses indicate the averaged residual for each time interval, blue crosses correspond to residuals displayed in b) and the green line displays the mean experiment residual. b) Spatial distribution of the tilt residual vectors at 𝑡 = 7 month (black arrows). The source is located by the orange square at the center of the domain.

Fig. 5 a

 5 Fig. 5 a) Relation between all components (x and y components) of drift rates 𝑎 𝑎 ⃗⃗⃗⃗ and model coefficient 𝛼 𝑧𝑧 ⃗⃗⃗⃗⃗⃗ associated to the 𝜀 𝑧𝑧 component of the strain tensor. Black crosses represent the true synthetic solution. After the first step of optimization, the admissible solution given by the blue dots shows a correlation (blue line). Correcting only this correlation leads to a distribution (green diamonds) largely dependent on the components of 𝛼 𝑧𝑧 ⃗⃗⃗⃗⃗⃗ . In addition to this correlation coefficient, the correction involves the covariance of the components of 𝛼 as provided by Eq. 7. The resulting distribution of 𝑎 * ⃗⃗⃗⃗ (red inverted triangles) is totally uncorrelated with elements of 𝛼 associated to all six components of the strain tensor. b) Enlargement of the grey rectangle in a).

Fig. 6

 6 Fig. 6 Evolution of the 6-strain tensor components over time. For each component, we represent in the upper graphs, the longterm evolutions and resolutions, and in the lower graph using blue dots, the upper bound of the instantaneous resolution 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ . The target values are indicated by the black line, the admissible solution 𝑝 𝑖𝑎 from the inversion by the blue line and the target solution 𝑝 𝑖 * from the correction by the red line. Ranges of long-term sensitivity, 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ (green areas) bound the corrected solution for each component.

Figure

  Figure 7b represent the instantaneous resolution of the volume variation reaching a maximum of ±490 m 3 at 𝑡 = 11.

Fig. 8

 8 Fig. 8 Schematic representation of two fractured domains in 2-D. Blue and red circles represent microseismic events defining fracking planes (blue and red line). a) The strain tensor given by the inversion is represented by black arrows around the elementary volume. Knowing the geometry and orientation of the fracture (blue line), the strain tensor is interpreted as an opening along the fracture plane (blue arrows). b) Considering fracture orientation from a priori knowledge (red line), the strain tensor represented by black and white arrows around the elementary volume can be interpreted as a strike slip (red arrows).

𝑐𝑛

  ⃗⃗⃗⃗ and the probable lack of convergence of the inversion due to data uncertainties. Regarding drift rates perturbations 𝐷 ⃗ ⃗ , we only have drift values 𝑎 𝑖 * (𝑖 = 𝑥, 𝑦) from the optimization process and the gradient of the functional for each drift parameters 𝑑𝑗 𝑖 = 𝑑𝐽 𝑑𝑎 𝑖

Acknowledgments

The post-doctoral research project of S. Furst was fully supported by the Total Company. We thank Maurizio Battaglia, Gábor Papp and an anonymous reviewer for their advice and suggestions, which improved the quality, value and clarity of paper. We also thank Megan Cook for her careful proofreading of this paper.

Data availability statement

The synthetic datasets generated and analysed during the current study are available on request from the corresponding author.

Author contribution statement J.C. devised the main conceptual ideas. S.F. and M.P. developed the theory and B.M. worked out the global algorithm. S.F. designed and performed the numerical calculations for the experiments. J.C., B.M. and M.P.

supervised the findings of this work. S.F. wrote the manuscript and designed the figures with support from M.P., J.C. and B.M. All authors provided critical feedback and discussed the results.

Appendix A: parameters resolution analysis

The approach leading to parameter resolution analysis [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]) can be extended to a generic formulation.

Considering small perturbations, Eq. 4 gives:

Eq. A.1 leads to 2𝑁 equations for both 𝑥 and 𝑦 components of tilt data and instrumental drift and can be rewritten as:

where 𝑆 (2𝑁) and 𝑃 ⃗ (6) are vectors of source signal and strain parameter perturbations at a given time, 𝐷 ⃗ ⃗ (2𝑁) are perturbations on the drift rates parameters and 𝑨 (2𝑁, 6) includes all components of the deformation model vector.

To assess the parameter resolution 𝑃 ⃗ , we invert the matrix 𝑨 using a least-square method. Using Eq. A.2, we define the confidence intervals of the parameters 𝑃 ⃗ as the sum of an instantaneous resolution, 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ = (𝑨 𝑡 𝑨) -1 𝑨 𝑡 𝑆 and a long-term one 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ = (𝑨 𝑡 𝑨) -1 𝑨 𝑡 (-𝐷 ⃗ ⃗ 𝑡):

𝑃 ⃗ = 𝑃 𝐼𝑇 ⃗⃗⃗⃗⃗ + 𝑃 𝐿𝑇 ⃗⃗⃗⃗⃗⃗ (A.3)