Didier Henrion

Jean-Bernard Lasserre

Inner approximations for polynomial matrix inequalities and robust stability regions

Keywords: polynomial matrix inequality, linear matrix inequality, robust optimization, robust fixed-order controller design, moments, positive polynomials

Following a polynomial approach, many robust fixed-order controller design problems can be formulated as optimization problems whose set of feasible solutions is modelled by parametrized polynomial matrix inequalities (PMI). These feasibility sets are typically nonconvex. Given a parametrized PMI set, we provide a hierarchy of linear matrix inequality (LMI) problems whose optimal solutions generate inner approximations modelled by a single polynomial superlevel set. Those inner approximations converge in a well-defined analytic sense to the nonconvex original feasible set, with asymptotically vanishing conservatism. One may also impose the hierarchy of inner approximations to be nested or convex. In the latter case they do not converge any more to the feasible set, but they can be used in a convex optimization framework at the price of some conservatism. Finally, we show that the specific geometry of nonconvex polynomial stability regions can be exploited to improve convergence of the hierarchy of inner approximations.

Introduction

Linear system stability can be formulated semialgebraically in the space of coefficients of the characteristic polynomial. The region of stability is generally nonconvex in this space, and this is a major obstacle when solving fixed-order and/or robust controller design problems. Using the Hermite stability criterion, these problems can be formulated as parametrized polynomial matrix inequalities (PMIs) where parameters account for uncertainties and the decision variables are controller coefficients. Recent results on real algebraic geometry and generalized problems of moments can be used to build up a 1 CNRS; LAAS; 7 avenue du colonel Roche, F-31077 Toulouse; France. henrion@laas.fr 2 Université de Toulouse; UPS, INSA, INP, ISAE; UT1, UTM, LAAS; F-31077 Toulouse; France 3 Faculty of Electrical Engineering, Czech Technical University in Prague, Technická 2, CZ-16626 Prague, Czech Republic 4 Institut de Mathématiques de Toulouse, Université de Toulouse; UPS; F-31062 Toulouse, France.

hierarchy of convex linear matrix inequality (LMI) outer approximations of the region of stability, with asymptotic convergence to its convex hull, see e.g. [START_REF] Henrion | Solving nonconvex optimization problems -How Glop-tiPoly is applied to problems in robust and nonlinear control[END_REF] for a software implementation and examples, and see [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF] for an application to PMI problems arising from static output feedback design.

If outer approximations of nonconvex semialgebraic sets can be readily constructed with these LMI relaxations, inner approximations are much harder to obtain. However, for controller design purposes, inner approximations are essential since they correspond to sufficient conditions and hence guarantees of stability or robust stability. In the robust systems control literature, convex inner approximations of the stability region have been proposed in the form of polytopes [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF], ellipsoids [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF] or more general LMI regions [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF][START_REF] Karimi | Robust control of polytopic systems by convex optimization[END_REF] derived from polynomial positivity conditions. Interval analysis can also be used in this context, see e.g. [START_REF] Walter | Guaranteed characterization of stability domains via set inversion[END_REF].

In this paper we provide a numerical scheme for approximating from inside the feasible set P ⊂ R n of a parametrized PMI P (x, u) 0 (for some matrix polynomial P ), that is, the set of points x such that P (x, u) 0 for all values of the parameter u in some specified domain U ⊂ R p (assumed to be a basic compact semialgebraic set1 ). This includes as a special case the approximation of the stability region (and the robust stability region) of linear systems. The particular case where P (x, u) is affine in x covers parametrized LMIs with many applications in robust control, as surveyed e.g. in [START_REF] Scherer | LMI relaxations in robust control[END_REF].

Given a compact set B ⊂ R n containing P, this numerical scheme consists of building up a sequence of inner approximations G d ⊂ P ⊂ B, d ∈ N, which fulfils two essential conditions:

1. The approximation converges in a well-defined analytic sense ;

2. Each set G d is defined in a simple manner, as a superlevel set of a single polynomial.

In our mind, this feature is essential for a successful implementation in practical applications.

More precisely, we provide a hierarchy of inner approximations (G d ) of P, where each G d = {x ∈ B : g d (x) ≥ 0} is a basic semi-algebraic set for some polynomial g d of degree d.

The vector of coefficients of the polynomial g d is an optimal solution of an LMI problem. When d increases, the convergence of (G d ) to P is very strong. Indeed, the Lebesgue volume of G d converges to the Lebesgue volume of P. In fact, on any (a priori fixed) compact set B, the sequence (g d ) converges for the L 1 -norm on B to the function x → λ min (x) = min u∈U λ min (x, u) where λ min (x, u) is the minimum eigenvalue of the matrixpolynomial P (x, u) associated with the PMI. Consequently, g d → λ min in (Lebesgue) measure on B, and g d k → λ min almost everywhere and almost uniformly on B, for a subsequence (g d k ). In addition, if one defines the piecewise polynomial ḡd := max k≤d g k , then ḡd → λ min almost everywhere, almost uniformly and in (Lebesgue) measure on B.

In addition, we can easily enforce that the inner approximations (G d ) are nested and/or convex. Of course, for the latter convex approximations, convergence to P is lost if P is not convex. However, on the other hand, having a convex inner approximation of P may reveal to be very useful, e.g., for optimization purposes.

On the practical and computational sides, the quality of the approximation of P depends heavily on the chosen set B ⊃ P on which to make the approximation of the function λ min . The smaller B, the better the approximation. In particular, it is worth emphasizing that when the set P to approximate is the stability or robust stability region of a linear system, then its particular geometry can be exploited to construct a tight bounding set B. Therefore, a good approximation of P is obtained significantly faster than with an arbitrary set B containing P.

Finally, let us insist that the main goal of the paper is to show that it is possible to provide a tight and explicit inner approximation with no quantifier, of nonconvex feasible sets described with quantifiers. Then this new feasible set can be used for optimization purposes and we are facing two cases:

• the convex case: if f and -g are convex polynomials,

B = {x ∈ R n : x ∞ ≤ 1} and G = {x ∈ B : g(x) ≥ 0} then the optimization problem min x f (x) s.t.
x ∈ G is polynomially solvable. Indeed, functions f (x), g(x), x ∞ are polynomially computable, of polynomial growth, and the feasible set is polynomially bounded. Then polynomial solvability of the problem follows from [3, Theorem 5.3.1].

• the nonconvex case: if -g is not convex then notice that firstly we still have an optimization problem with no quantifier, a nontrivial improvement. Secondly we are now faced with an polynomial optimization problem with a single polynomial constraint and possibly bound constraints x ∈ B. One may then apply the hierarchy of convex LMI relaxations described in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]Chapter 5]. Of course, in general, polynomial optimization is NP-hard. However, if the size of the problem is relatively small and the degree of g is small, practice seems to reveal that the problem is solved exactly with few relaxations in many cases, see [11, §5.3.3]. In addition, if some structured sparsity in the data is present then one may even solve problems of potentially large size by using an appropriate sparse version of these LMI relaxations as described in [START_REF] Waki | Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity[END_REF], see also [11, §4.6].

The outline of the paper is as follows. In Section 2 we formally state the problem to be solved. In Section 3 we describe our hierarchy of inner approximations. In Section 4, we show that the specific geometry of the stability region can be exploited, as illustrated on several standard problems of robust control. The final section collects technical results and the proofs.

Problem statement

Let R[x] denote the ring or real polynomials in the variables x = (x 1 , . . . , x n ), and let R[x] d be the vector space of real polynomials of degree at most d.

Similarly, let Σ[x] ⊂ R[x]
denote the convex cone of real polynomials that are sums of squares (SOS) of polynomials, and Σ[x] d ⊂ Σ[x] its subcone of SOS polynomials of degree at most 2d. Denote by S m the space of m × m real symmetric matrices. For a given matrix A ∈ S m , the notation A 0 means that A is positive semidefinite, i.e., all its eigenvalues are real and nonnegative.

Let P : R[x, u] → S m be a matrix polynomial, i.e. a matrix whose entries are scalar multivariate polynomials of the vector indeterminates x and u. Then

P := {x ∈ R n : ∀u ∈ U, P (x, u) 0} (1)
defines a parametrized polynomial matrix inequality (PMI) set, where x ∈ R n is a vector of decision variables, u ∈ R p is a vector of uncertain parameters belonging to a compact semialgebraic set

U := {u ∈ R p : a i (u) ≥ 0, i = 1, . . . , n a } (2)
described by given polynomials a i (u) ∈ R[u], and P (x, u) is a given symmetric polynomial matrix of size m. As U is compact, without loss of generality we assume that for some

i = i * , a i * (u) = R 2 -u T u,
where R is sufficiently large.

We also assume that P is bounded and that we are given a compact set B ⊃ P with explicitly known moments y = (y α ), α ∈ N n , of the Lebesgue measure on B, i.e.

y α := B x α dx (3)
where x α := n i=1 x α i i . Typical choices for B are a box or a ball. To fix ideas, let

B := {x ∈ R n : b j (x) ≥ 0, j = 1, . . . , n b }
for some polynomials b j ∈ R[x]. Again, with no loss of generality, we may and will assume that for some j = j * , b j * (x) = R 2 -x T x, where R is sufficiently large. Finally, denote by vol A the Lebesgue volume of any Borel set A ⊂ B.

We are now ready to state our polynomial inner approximation problem.

Problem 1 (Inner Approximations) Given set P, build up a sequence of basic closed semialgebraic sets

G d = {x ∈ B : g d (x) ≥ 0}, for some g d ∈ R[x], such that G d ⊆ P, d = 1,
G d = {x ∈ B : g d (x) ≥ 0}, for some g d ∈ R[x], such that G d ⊆ G d+1 ⊆ P, d = 1, 2, . . .

A hierarchy of semialgebraic inner approximations

Given a polynomial matrix P (x, u) which defines the set P in (1), polynomials a i ∈ R[u] which define the uncertain set U in (2), let V = {v ∈ R m : v T v = 1} denote the Euclidean unit sphere of R m and let λ min : B → R be the function:

x → λ min (x) = min u∈U min v∈V v T P (x, u)v (4)
as the robust minimum eigenvalue function of P (x, u). Function λ min is continuous but not necessarily differentiable. It allows to define set P alternatively as the superlevel set

P = {x ∈ R n : λ min (x) ≥ 0}.

Primal SOS SDP problems

Let a 0 ∈ R[u] be the constant polynomial 1. Let 2d 0 ≥ max(2+deg P, max i dega i , max j degb j ), and consider the hierarchy of convex optimization problems indexed by the parameter d ∈ N, d ≥ d 0 :

ρ d = B λ min (x) dx -min g,r,s,t B g(x) dx s.t. v T P (x, u)v -g(x) = r(x, u, v)(1 -v T v) + na i=0 s i (x, u, v)a i (u) + n b j=1 t j (x, u, v)b j (x) ∀(x, u, v) (5)
where decision variables are coefficients of polynomials g ∈ R[x] 2d , r ∈ R[x, u, v] 2dr and coefficients of SOS polynomials s i ∈ Σ[x, u, v] ds i , i = 0, 1, . . . , n a , and t j ∈ Σ[x, u, v] dt j , j = 1, . . . , n b . Note in particular that the degrees of the polynomials should be such that

d r ≥ d -1, d s i ≥ d -⌈(deg a i )/2⌉ and d t j ≥ d -⌈(deg b j )/2⌉
. Since higher degree terms may cancel, the degrees can be chosen strictly greater than these lower bounds. However, in the experiments described later on in the paper, we systematically chose the lowest possible degrees.

For each d ∈ N fixed, the associated optimization problem ( 5) is a semidefinite programming (SDP) problem. Indeed, stating that the two polynomials in both sides of the equation in ( 5) are identical translates into linear equalities between the coefficients of polynomials g, r, (s i ), (t j ) and stating that some of them are SOS translates into semidefiniteness of appropriate symmetric matrices. For more details, the interested reader is referred to e.g. [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]Chapter 2].

Dual moment SDP problems

To define the dual to SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] we must introduce some notations.

With a sequence y = (y α ), α ∈ N n , let L y : R[x] → R be the linear functional

f (= α f α x α ) → L y (f ) = α f α y α , f ∈ R[x].
With d ∈ N, the moment matrix of order d associated with y is the real symmetric matrix M d (y) with rows and columns indexed in N n d , and defined by

M d (y)(α, β) := L y (x α+β ) = y α+β , ∀α, β ∈ N n d . (6)
A sequence y = (y α ) has a representing measure if there exists a finite Borel measure µ on R n , such that y α = x α dµ for every α ∈ N n .

With y as above and h ∈ R[x], the localizing matrix of order d associated with y and h is the real symmetric matrix M d (h y) with rows and columns indexed by N n d , and whose entry (α, β) is given by

M d (y)(h y)(α, β) := L y (h(x) x α+β ) = γ h γ y α+β+γ , ∀α, β ∈ N n d . (7)
With these notations, the dual to SDP problem ( 5) is given by:

ρ * d = B λ min (x)dx -min y L y (v T P (x, u)v) s.t. M d (y) 0, M d-1 ((1 -v T v) y) = 0 M d-da i (a i y) 0, i = 0, 1, . . . , n a M d-d b j (b j y) 0, j = 1, . . . , n b L y (x α ) = B x α dx, ∀α ∈ N n 2d (8)
where y ∈ N n+p+m 2d .

Convergence

Before stating our main results, let us recall some standard notions of functional analysis. Let g : B → R be a function of x, and let (g d ) denote a sequence of functions of x indexed by d ∈ N. Lebesgue space L 1 (B) is the Banach space of integrable functions on B equipped with the norm

g 1 = B |g|dx.
Regarding sequence (g d ), we use the following notions of convergence in B when d → ∞:

• g d → g in L 1 norm means lim d→∞ g -g d 1 = 0; • g d → g in Lebesgue measure means that for every ε > 0, lim d→∞ vol{x : |g(x) -g d (x)| ≥ ε} = 0;
• g d → g almost everywhere means that lim d→∞ g d (x) = g(x) pointwise except possibly for x ∈ A ⊂ B with vol A = 0;

• g d → g almost uniformly means that for every given ε > 0, there is a set A ⊂ B such that vol A < ε and g d → g uniformly on B \ A;

• finally, with the notation g d ↑ g we mean that g d → g while satisfying g d (x) ≤ g d+1 (x) for all d.

For more details on these related notions of convergence, see [1, §2.5].

Lemma 1 For every d ≥ d 0 , SDP problem ( 5) has an optimal solution g d ∈ R[x] 2d and

ρ d = B (λ min (x) -g d (x)) dx = λ min -g d 1 . (9)
A detailed proof of Lemma 1 can be found in §6.1. In particular we prove that there is no duality gap between SOS SDP problem (5) and moment SDP problem (8), i.e. ρ d = ρ * d . For every d ≥ d 0 , let ḡd : B → R be the piecewise polynomial

x → ḡd (x) := max d 0 ≤k≤d g k (x). ( 10
)
We are now in position to prove our main result.

Theorem 1 Let g d ∈ R[x]
2d be an optimal solution of SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] and consider the associated sequence

(g d ) ⊂ L 1 (B) for d ≥ d 0 . Then: (a) g d → λ min in L 1 norm

Polynomial and piecewise polynomial inner approximations

Corollary 1 For every d ≥ d 0 , let g d ∈ R[x] 2d be an optimal solution of SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF], let ḡd be the piecewise polynomial defined in [START_REF] Lasserre | The multi-dimensional version of b a x p dx[END_REF], and let

G d := {x ∈ B : g d (x) ≥ 0}, Ḡd := {x ∈ B : ḡd (x) ≥ 0}. ( 11
)
Then

G d ⊂ P ∀ d ≥ d 0 and lim d→∞ vol(P \ G d ) = 0. ( 12
) Ḡd 0 ⊆ • • • ⊆ Ḡd ⊆ • • • ⊂ P and lim d→∞ vol(P \ Ḡd ) = 0. ( 13
)
That is, sequence (G d ) solves Problem 1 and sequence ( Ḡd ) solves Problem 2 if piecewise polynomials are allowed.

A proof can be found in §6.3.

Nested polynomial inner approximations

We now consider Problem 2 where g d is constrained to be a polynomial instead of a piecewise polynomial. We need to slightly modify SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF]. Suppose that at step d -1 in the hierarchy we have already obtained an optimal solution

g d-1 ∈ R[x] 2d-2 , such that g d-1 ≥ g d 0 on B, for all d 0 ≤ d -1.
At step d we now solve SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] with the additional constraint

g(x) -g d-1 (x) = c 0 (x) + n b j=1 c j (x)b j (x), ∀x (14)
with unknown SOS polynomials

c 0 ∈ Σ[x] d and c j ∈ Σ[x] d-d b j . Corollary 2 Let g d ∈ R[x]
2d be an optimal solution of SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] with the additional constraint [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF] and let G d be as in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]

for d ≥ d 0 . Then the sequence (G d ) solves Problem 2.
For a proof see §6.4.

Convex nested polynomial inner approximations

Finally, for g ∈ R[x] 2d , denote by ∇ 2 g(x) the Hessian matrix of g at x, and consider SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] with the additional constraint

v T ∇ 2 g(x)v = c 0 (x, v) + n b j=1 c j (x, v)b j (x) + c n b +1 (x, v)(1 -v T v), (15)
for some SOS polynomials

c 0 ∈ Σ[x, v] d , c j ∈ Σ[x, v] d-d b j and c n b +1 ∈ Σ[x, v] d-1 .
Corollary 3 Let g ∈ R[x] 2d be an optimal solution of SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] with the additional constraint [START_REF] Scherer | LMI relaxations in robust control[END_REF] and let G d be as in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]

for d ≥ d 0 . Then the sequence (G d ) solves Problem 3.
The proof follows along the same lines as the proof of Corollary 2.

Example

Consider the nonconvex planar PMI set

P = {x ∈ R 2 : P (x) = 1 -16x 1 x 2 x 1 x 1 1 -x 2 1 -x 2 2 0}
which is Example II-E in [START_REF] Henrion | Convergent relaxations of polynomial matrix inequalities and static output feedback[END_REF] scaled to fit within the unit box whose moments (3) are readily given by

B = {x ∈ R 2 : x ∞ ≤ 1} -1 -0.
y α = 4 (α 1 + 1)(α 2 + 1)
.

On Figure 1 we represent the degree two and degree four solutions to SDP problem (5), modelled by YALMIP 3 and solved by SeDuMi 1.3 under a Matlab environment. We see in particular that the degree four polynomial superlevel set G 2 is somewhat smaller than expected. This is due to the fact that the objective function in problem ( 5) is the integral of g(x) over the whole box B, not only over PMI set P. There is a significant role played by the components of the integral on complement set B\P, and this deteriorates the inner approximation.

This issue can be addressed partly by embedding P in a tighter set B, for example here the unit disk

B = {x ∈ R 2 : x 2 ≤ 1}
whose moments (3) are given by

y α = Γ( α 1 +1 2 )Γ( α 2 +1 2 ) Γ(2 + α 1 +α 2 2 )
where Γ is the gamma function such that Γ(k) = (k -1)! for integer k. See [12, Theorem 3.1] for the general expression2 of moments of the unit disk in R n .

On Figure 2 we represent the degree two and degree four solutions to SDP problem [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF]. Comparing with Figure 1, we see that the approximations embedded in the unit disk are much tighter than the approximations embedded in the unit box. Finally, on Figure 3 we represent the tighter degree six and degree eight inner approximations within the unit disk.

Geometry of control problems

As explained in the introduction, inner approximations of the stability regions are essential for fixed-order controller design. The PMI regions arising from parametric stability conditions have a specific geometry that can be exploited to improve the convergence of the hierarchy of inner approximations. In this section, we first recall Hermite's PMI formulation of (discrete-time) stability conditions. Then we recall that the PMI stability region is the image of a unit box through a multi-affine mapping, which allows to derive explicit expressions for the moments of the full-dimensional stability region, as well as tight polytopic outer approximations of low-dimensional affine sections of the stability region. Numerical examples illustrate these techniques for fixed-order nominal and robustly stabilizing controller design.

Hermite's PMI

Derived by the French mathematician Charles Hermite in 1854, the Hermite matrix criterion is a symmetric version of the Routh-Hurwitz criterion for assessing stability of a polynomial. Originally it was derived for locating the roots of a polynomial in the open upper half of the complex plane, but with a fractional transform it can be readily transposed to the open unit disk and discrete-time stability. The criterion says that a polynomial x(z

) = z n + x 1 z n-1 + • • • + x n-1 z +
x n has all its roots in the open unit disk if and only if its Hermite matrix

P (x) = T T 1 (x)T 1 (x) -T T 2 (x)T 2 (x) is positive definite, where T 1 (x) =      1 x 1 x 2 0 1 x 1 0 0 1 . . .      T 2 (x) =      x n x n-1 x n-2 0 x n x n-1 0 0 x n . . .     
are n-by-n upper-right triangular Toeplitz matrices, see e.g. the entrywise formulas of [2, Theorem 3.13] or the construction explained in [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF]. The Hermite matrix is n-by-n, symmetric and quadratic in coefficients x = (x 1 , x 2 , . . . , x n ), so that the interior of the PMI set

P = {x ∈ R n : P (x) 0}
is the parametric stability domain which is bounded, connected but nonconvex for n ≥ 3. Optimal controller design amounts to optimizing over semialgebraic set P.

Multiaffine mapping of the unit box

As explained e.g. in [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF] or [16, §3.5] and references therein, stability domain P can also be constructed as the image of the unit box (in the space of so-called reflection coefficients) through a multiaffine mapping. More explicitly

P = f (K) where K = {k ∈ R n : k ∞ ≤
1} and multiaffine mapping f : R n → R n is defined by

f (k) =   0 1 0 0 0 0 0 1 0 0 0 0 0 1 0         1 0 0 k 3 0 1 k 3 0 0 k 3 1 0 k 3 0 0 1 0 0 0 0           1 0 k 2 0 1 + k 2 0 k 2 0 1 0 0 0       1 k 1 k 1 1 0 0   1 0 =   k 2 k 3 + k 1 (1 + k 2 ) k 2 + k 1 k 3 (1 + k 2 ) k 3  
in the case n = 3. The general expression of f for other values of n is not given here for space reasons, but it follows readily from the construction outlined above.

Using this mapping we can obtain moments (3) of B = P analytically:

y α = P x α dx = K (k 2 k 3 + k 1 (1 + k 2 )) α 1 (k 2 + k 1 k 3 (1 + k 2 )) α 2 k α 3 3 det ∇f (k)dk (16)
where det ∇f

(k) = (1 + k 2 )(1 -k 2 3
) is the determinant of the Jacobian of f , in the case n = 3. For space reasons we do not give here the explicit value of y α as a function of α, but it can be obtained by integration by parts.

Finally, let us mention a well-known geometric property of P: its convex hull is a polytope whose vertices correspond to the n + 1 polynomials with roots equal to -1 or +1. For example, when n = 3, we have conv P = conv{(-3, 3, -1), (-1, -1, 1), (1, -1, -1), (3, 3, 1)}.

(17)

Third degree stability region

Consider the problem of approximating from the inside the nonconvex stability region P of a discrete-time third degree polynomial z → z 3 + x 1 z 2 + x 2 z + x 3 . An ellipsoidal inner approximation was proposed in [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF]. The Hermite polynomial matrix defining P as in ( 1) is given by

P (x) =   1 -x 2 3 x 1 -x 2 x 3 x 2 -x 1 x 3 x 1 -x 2 x 3 1 + x 2 1 -x 2 2 -x 2 3 x 1 -x 2 x 3 x 2 -x 1 x 3 x 1 -x 2 x 3 1 -x 2 3   .
The boundary of P consists of two triangles and a hyperbolic paraboloid. The convex hull of P is the simplex described in [START_REF] Waki | Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity[END_REF]. We have analytic expressions [START_REF] Shcherbakov | On the generation of random stable polynomials[END_REF] for the moments (3) of B = P.

On Figures 4, 5 and 6 we respectively represent the degree two, four and six inner approximations of P, scaled within the unit box for visualization purposes. We observe that the degree six approximation is very tight, thanks to the availability of the moments of the Lebesgue measure on P.

Fixed-order controller design

Consider the linear discrete-time system with characteristic polynomial z → z 4 -(2x 1 + x 2 )z 3 + 2x 1 z + x 2 depending affinely on two real design parameters x 1 and x 2 . It follows from Hermite's stability criterion that this polynomial has its roots in the open unit disk if and only if

P (x) =    1 -x 2 2 -2x 1 -x 2 -2x 1 x 2 0 2x 1 + 2x 1 x 2 + x 2 2 -2x 1 -x 2 -2x 1 x 2 1 + 4x 1 x 2 -2x 1 -x 2 -2x 1 x 2 0 0 -2x 1 -x 2 -2x 1 x 2 1 + 4x 1 x 2 -2x 1 -x 2 -2x 1 x 2 2x 1 + 2x 1 x 2 + x 2 2 0 -2x 1 -x 2 -2x 1 x 2 1 -x 2 2   
is positive definite. As recalled in (4.2), the convex hull of the four-dimensional stability domain of a degree four polynomial is the simplex with vertices (-4, 6, -4, 1), (-2, 0, 2, -1), (0, -2, 0, 1), (2, 0, -2, -1), (4, 6, 4, 1) corresponding to the five polynomials with zeros equal to -1 or +1. Using elementary linear algebra, we find out that the image of this simplex through the affine mapping (-(

2x 1 + x 2 ), 0, 2x 1 , x 2 ) parametrized by x ∈ R 2 is the two-dimensional simplex B = conv{(- 1 4 , 1), ( 7 8
, -

), (-

5 8 , - 1 2
)}.

The (closure of the) stability region P = {x ∈ R 2 : P (x) 0} is therefore included in B, whose moments (3) are readily obtained e.g. by the explicit formulas of [START_REF] Lasserre | The multi-dimensional version of b a x p dx[END_REF]. On Figures 7 and8 we represent the degree two, four, six and eight inner approximations to P, corresponding to stability regions for the linear system. We observe that the approximations become tight rather quickly. This is due to the fact that B is a good outer approximation of P with known moments. Tighter outer approximations B would result in tighter inner approximations of P, but then the moments of B can be hard to compute, see [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF].

Robust controller design

Now consider the uncertain polynomial z →

x 2 + u + 2x 1 z -(2x 1 + x 2 )z 3 + z 4 with u ∈ U = {u ∈ R : u 2 ≤ 1 16
} with uncertain Hermite matrix P (x, u) and the corresponding parametrized PMI stability region P in [START_REF] Ash | Real analysis and probability[END_REF]. Let us use the same bounding set B as in §4.4. On Figure 9 we represent the degree two and degree four inner approximations to P, corresponding to robust stability regions for the linear system. Comparing with Figure 7 we see that the approximations are smaller, and in particular they do not touch the stability boundary to cope with the robustness requirements.

Conclusion

We have constructed a hierarchy of inner approximations of feasible sets defined by parametrized or uncertain polynomial matrix inequalities (PMI). Each inner approximation is computed by solving a convex linear matrix inequality (LMI) problem. The hierarchy converges in a well-defined analytic sense, so that conservatism of the approximation is guaranteed to vanish asymptotically. In addition, the inner approximations are simple polynomial or piecewise-polynomial superlevel sets, so that optimization over these sets is significantly simpler than optimization over the original parametrized PMI set. In particular, we remove the possibly complicated dependence of the problem data on the uncertain parameters.

One may also impose the hierarchy of inner approximations to be nested. Finally, one may also impose the inner approximations to be convex. In this latter case they do not converge any more to the feasible set but, on the other hand, optimization over the parametrized PMI set can be reformulated as a convex polynomial optimization problem (of course at the price of some conservatism). Ideally, beyond convexity, we may also want the inner convex approximation to be semidefinite representable (as an explicit affine projection of an affine section of the SDP cone), and deriving such a representation may be an interesting research direction.

The tradeoff to be found is between tightness of the inner approximation and degree of the defining polynomials. In the context of robust control design, a satisfactory inner approximation can be possibly computed off-line, and then used afterwards on-line in a feedback control setup.

Our methodology is valid for general parametrized PMI problems. However, in the case of parametrized PMI problems coming from fixed-order robust controller design problems, geometric insight can be exploited to improve convergence of the hierarchy. The key information is the knowledge of the moments of the Lebesgue measure on a compact set which tightly contains the parametrized PMI set we want to approximate from the inside. In turns out that for robust control problems this knowledge is available easily, as illustrated in the paper by several examples.

The main limitation of the approach lies in the ability of solving primal moment and dual polynomial sum-of-squares LMI problems. State-of-the-art general-purpose semidefinite programming solvers can currently address problems of relatively moderate dimensions, but problem structure and data sparsity can be exploited for larger problems.

6 Appendix

Proof of Lemma 1

Proof: The dual of polynomial SOS SDP problem ( 5) is moment SDP problem (8). Slater's condition cannot hold for (8) because V has empty interior in R m . However it turns out that Slater's condition holds for an equivalent version of SDP problem (8), i.e., the latter has a strictly feasible solution ŷ. Indeed, let J ⊂ R[v] be the ideal generated by the polynomial v → θ(v) := 1 -v T v so that the real variety V R (J) := {v ∈ R m : θ(v) = 0} associated with J is just the unit sphere V. It turns out that the real radical3 of J is J itself, that is, I(V R (J)) = J (where for S ⊂ R m , I(S) denotes the vanishing ideal). And after embedding J in R[x, u, v], we still have I(V R (J)) = J.

Let H := {(α, β, γ) ∈ N n × N p × N m : γ m ≤ 1}, and let H d := {(α, β, γ) ∈ H : i α i + j β j + ℓ γ ℓ ≤ d}. The monomials (x α u β v γ ), (α, β, γ) ∈ H, form a basis of the quotient space R[x, u, v]/J. Moreover, for every (α, β, γ) ∈ N n+p+m d , x α u β v γ = (a,b,c)∈H d p abc x a u b v c + h(x, u, v) ∈R[x,u,v] d-2 (1 -v T v),
for some real coefficients (p abc ), and some h ∈ R[x, u, v] d-2 . This is because every time one sees a monomial x a u b v c with c m ≥ 2, one uses v 2 m = 1j =m v 2 j to reduce this monomial modulo θ = (1 -v T v). For instance

x a u b v c 1 1 • • • v c m-1 m-1 v 3 m = x a u b v c 1 1 • • • v c m-1 m-1 v m × v 2 m =-θ+(1-j =m v 2 j ) = x a u b v c 1 1 • • • v c m-1 m-1 v m - j =m x a u b v c 1 1 • • • v c j +2 j • • • v c m-1 m-1 v m -x a u b v c 1 1 • • • v c m-1 m-1 v m ∈R[x,u,v] d-2 θ(v),
etc. Therefore, for every (α, β, γ), (α

′ , β ′ , γ ′ ) ∈ H d , x α+α ′ u β+β ′ v γ+γ ′ = (a,b,c)∈H 2d p abc x a u b v c + h(x, u, v) ∈R[x,u,v] 2d-2 (1 -v T v), (18)
for some real coefficients (p abc ), and some

h ∈ R[x, u, v] 2d-2 .
So because of the constraints M d-1 ((1 -v T v) y) = 0, the semidefinite program ( 8) is equivalent to the semidefinite program: Next, let ŷ be the sequence of moments of the (product) measure µ uniformly distributed on B × U × V, and scaled so that for all (α, β, γ)

ρ * d = B λ min (x)dx -min y L y (v T P (x, u)v) s.t. Md (y) 0, M d-1 ((1 -v T v) y) = 0 Md-da i (a i y) 0, i = 0,
∈ N n+p+m 2d ŷαβγ = B×U×V x α u β v γ dµ(x, u, v) = 1 vol U × V B U V x α u β v γ dx du dλ(v) dµ(x,u,v)
(with λ the rotation invariant measure on V). Therefore, for every α

∈ N n 2d , ŷα00 = L y (x α ) = B×U×V x α dµ(x, u, v) = B x α dx.
Moreover, M d-1 ((1-v T v) y) = 0 for every d and importantly, Md (ŷ) ≻ 0, Md-da i (a i ŷ) ≻ 0 and Md-d b j (b j ŷ) ≻ 0. To see why, suppose for instance that h T Md (ŷ)h = 0 for some vector h = 0. This means that for some non trivial polynomial h ∈ R[x, u, v]/J of degree d,

h T Md (ŷ)h = B×U×V h 2 dµ = 0, that is, h(x, u, v) = 0 for µ-almost all (x, u, v) ∈ B × U × V, and so h(x, u, v) = 0 for all (x, u, v) ∈ B × U × V because h is continuous. But as B × U has nonempty interior in R n × R p , then necessarily h ∈ I(V R (J)) (= J) -see Lemma 2 in section 6.5 -which contradicts 0 = h ∈ R[x, u, v]/J.
P (x, u)v -g(x) = r(x, u, v)(1 -v T v) + na i=0 s i (x, u, v)a i (u) + n b j=1 t j (x, u, v)b j (x) ∀(x, u, v) (20
) where now the decision variables are coefficients of polynomials g ∈ R[x] 2d , r ∈ R[x, u, v] 2dr , and coefficients of SOS polynomials s i ∈ Σ[x, u, v] da i , i = 0, 1, . . . , n a , and 20) is guaranteed to have an optimal solution (g * , r * , s * , t * ). But observe that such an optimal solution (g * , r * , s * , t * ) is also feasible in [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF], and so having value ρ ′ d = ρ d = ρ * d , (g * , r * , s * , t * ) is also an optimal solution of (5).

t j ∈ Σ[x, u, v] d b j , j = 1, . . . , n b . That is, ρ ′ d = ρ * d and so ρ d = ρ * d because ρ ′ d ≤ ρ d ≤ ρ * d . If ρ ′ d < ∞ then (
It remains to prove that ρ d is bounded. For any feasible solution y of (5), y 0 ≤ vol B, and

L y (x 2d i ) ≤ R 2d y d 0 ; L y (u 2d j ) ≤ R 2d y d 0 ; L y (v 2d k ) ≤ y d 0 , (21)
for all i = 1, . . . , n, j = 1, . . . , p, k = 1, . . . , m. This follows from , which shows that the feasible set of ( 8) is compact. Hence (8) has an optimal solution and ρ d is finite; therefore its dual (5) also has an optimal solution, the desired result.

M d-da i * (a i * y) 0, M d-d b j * (b j * y) 0 and M d-1 ((1 -v T v) y) = 0, where a i * (x) = R 2 -x T x and b j * (x) = R 2 -u T u;

Proof of Theorem 1

Proof: (a) Let K := B × U × V ⊂ R n+p+m and consider the infinite-dimensional optimization problem

ρ = min µ∈M (K) K v T P (x, u)v dµ(x, u, v) s.t. K x α dµ = B x α dx, α ∈ N n (22)
where M(K) is the space of finite Borel measures on K. Problem ( 22) has an optimal solution µ * ∈ M(K). Indeed, ρ ≥ B λ min (x)dx because for every (x, u, v) ∈ K, v T P (x, u)v ≥ λ min (x); and so for every feasible solution µ ∈ M(K),

K v T P (x, u, v)v dµ(x, u, v) ≥ K λ min (x) dµ(x, u, v) = B λ min (x) dx because K x α dµ = B
x α dx for all α ∈ N and hence the marginal of µ on R n is the Lebesgue measure on B. On the other hand, observe that for every As ε > 0 was arbitrary and G d ⊂ P, we obtain the desired result [START_REF] Lasserre | Solving a class of multivariate integration problems via Laplace techniques[END_REF]. The proof of ( 13) is similar.

x ∈ B, λ min (x) = v T x P (x, u x )v x for some (u x , v x ) ∈ U × V. Therefore, let µ * ∈ M(K) be the Borel measure concentrated on (x, u x , v x ) for all x ∈ B, i.e. µ * (B ′ × U ′ × V ′ ) := B ′ ∩U ′ 1 U ′ ×V ′ (u x , v x ) dx, ∀(B ′ , U ′ , V ′ ) ∈ B(B) × B(U) × B(V)
f (x, u, v) = (v T v -1) h(x 0 , u 0 , v), ∀(x, u, v) ∈ R n × R p × R m ,
i.e., f = (v T v -1)h, which proves the desired result that f ∈ J.

2 Figure 2 :

 22 Figure 2: Degree two (left) and four (right) inner approximations (light gray) of PMI set (dark gray) embedded in unit disk (dashed).

2 Figure 3 :

 23 Figure 3: Degree six (left) and eight (right) inner approximations (light gray) of PMI set (dark gray) embedded in unit disk (dashed).

Figure 4 :

 4 Figure 4: Two views of a degree two inner approximation (red) of nonconvex third-degree stability region (gray).

Figure 5 :

 5 Figure 5: Two views of a degree four inner approximation (red) of nonconvex third-degree stability region (gray).

Figure 6 :

 6 Figure 6: Two views of a degree six inner approximation (red) of nonconvex third-degree stability region (gray).

Figure 7 :

 7 Figure 7: Degree two (left) and four (right) inner approximations (light gray) of PMI stability region (dark gray) embedded in simplex (dashed).

Figure 8 :

 8 Figure 8: Degree six (left) and eight (right) inner approximations (light gray) of PMI stability region (dark gray) embedded in simplex (dashed).

Figure 9 :

 9 Figure 9: Degree two (left) and four (right) inner approximations (light gray) of robust PMI stability region (dark gray) embedded in simplex (dashed).

1 ,

 1 . . . , n a Md-d b j (b j y) 0, j = 1, . . . , n b L y (x α ) = B x α dx, ∀α ∈ N n 2d (19)where the smaller moment matrix Md (y) is the submatrix of M d (y) obtained by looking only at rows and columns indexed in the monomial basis (x α y β v γ ), (α, β, γ) ∈ H d , instead of N n+p+m d . Similarly, the smaller localizing matrix Md-da i (a i y) is the submatrix of M d-da i (a i y) obtained by looking only at rows and columns indexed in the monomial basis (x α y β v γ ), (α, β, γ) ∈ H d-da i , instead of N n+p+m d-da i ; and similarly for Md-d b j (b j y).Indeed, in view of (18) and usingM d ((1 -v T v) y) = 0, every column of M d (y) associated with (α, β, γ) ∈ N n+p+m is a linear combination of columns associated with (α ′ , β ′ , γ ′ ) ∈ H d .And similary for M d-da i (a i y) and M d-d b j (b j y). Hence, M d (y) 0 ⇔ Md (y) 0, and M d-da i (a i y) 0 ⇔ Md-da i (a i y) 0; M d-d b j (b j y) 0 ⇔ Md-d b j (b j y) 0, for all i = 1, . . . , n a , j = 1, . . . , n b .

  see the comments after (2) and (3). Then by [13, Lemma 4.3], one obtains |y α | ≤ R 2d (vol B) d , for all α ∈ N n 2d

where x → 1 B 2 .

 12 (x) denotes the indicator function of set B and B(B) denotes the Borel σ-algebra of subsets of B. Then µ * is feasible for problem (22) with value K v T P (x, u)v dµ * (x, u, v) = B λ min (x) dx which proves that ρ = B λ min (x) dx.Next, λ min being continuous on compact set B, by the Stone-Weierstrass theorem [1, §A7.5], for every ε > 0 there exists a polynomialh ε ∈ R[x] such that sup x∈B |λ min (x) -h ε (x)| < ε Hence the polynomial p ε := h ε -ε satisfies λ min -p ε > 0 on B and so v T P (x, u)v -p ε > 0 on B × U × V.By Putinar's Positivstellensatz, see e.g [11, Section 2.5], there exists SOS polynomials r ε ∈ R[x, u, v], and s iε , t jε ∈ Σ[x, u, v] such that equation (5) is satisfied. Hence for d sufficiently large, say d ≥ d ε , (p ε , r ε , s iε , t jε ) is a feasible solution of (5) with associated value B(λ min (x) -p ε (x)) dx ≤ 3ε 2 B dx.Hence 0 ≤ ρ d ≤ 3ε 2 B dx whenever d ≥ d ε where ρ d is defined in[START_REF] Karimi | Robust control of polytopic systems by convex optimization[END_REF]. As ε > 0 was arbitrary, we obtain the desired result lim d→∞ ρ d = 0. Observe that since g d ≤ λ min for all d, ρ d = B (λ min (x) -g d (x)) dx = B |λ min (x) -g d (x)| dx so that the convergence ρ d → 0 is just the convergence g d → λ min for the L 1 norm on B. Finally the convergence g d → λ min in Lebesgue measure on B follows from [1, Theorem 2.5.1]. (b) For each x ∈ B, fixed and arbitrary, the sequence (ḡ d ) is monotone nondecreasing and bounded above by λ min . Therefore there exists g * : B → R such that for every x ∈ B, ḡd (x) ↑ g * (x) ≤ λ min (x) as d → ∞. Since ḡd ≥ ḡ0 and B ḡ0 dx > -∞, by Lebesgue's Dominated Convergence Theorem [1, §1.6.9] B g * (x)dx = lim d→∞ B ḡd (x)dx = B λ min (x)dx, and so from g * (x) ≤ λ min (x) we deduce that g * (x) = λ min (x) for almost all x ∈ B. Combining the latter with ḡd ↑ g * , we obtain that ḡd → λ min almost everywhere in B. But then since the Lebesgue measure is finite on B, by Egorov's theorem [1, Theorem 2.5.5], ḡd → λ min almost uniformly in B. Finally, convergence in Lebesgue measure on B also follows from [1, Theorem 2.5.2]. 6.3 Proof of Corollary 1 Proof: By Theorem 1, lim d→∞ λ min -g d 1 = 0. Therefore, by [1, Theorem 2.5.1] the sequence (g d ) converges to λ min in Lebesgue measure, i.e. for every ε > 0, lim d→∞ vol{x : |λ min (x) -g d (x)| ≥ ε} = 0. (23) Let ε > 0 be fixed, arbitrary, and let P ε := {x ∈ B : λ min (x) ≥ ε}, so that lim ε→0 vol P ε = vol P. By (23), lim d→∞ vol(P ε ∩ {x ∈ B : g d (x) < 0}) = 0. Next, for all d ∈ N, vol P ε = vol(P ε ∩ {x ∈ B : g d (x) < 0}) + vol(P ε ∩ {x ∈ B : g d (x) ≥ 0}). Therefore, taking the limit as d → ∞ yields vol P ε = lim d→∞ vol(P ε ∩ {x ∈ B : g d (x) < 0}) =0 by (23) + lim d→∞ vol(P ε ∩ {x ∈ B : g d (x) ≥ 0} ε ∩ G d ) ≤ lim d→∞ vol G d .

  etc., to finally obtainv α = p α (v)(1 -v T v) + r α , ∀α ∈ N m d ,for some p α ∈ R[v] |α|-2 and r α ∈ R[v]/J. Therefore, summing up over all α ∈ N m d yields:f (x 0 , u 0 , v) = α∈N m d g α (x 0 , u 0 ) v α = (v T v -1) α∈N m d g α (x 0 , u 0 ) p α (v) h(x 0 ,u 0 ,v) + α∈N m d g α (x 0 , u 0 ) r α (v) = 0 as v → f (x 0 , u 0 , v) ∈ J = (v T v -1) h(x 0 , u 0 , v), ∀v ∈ R m ,(25)for some h ∈ R[x, u, v]. But since (25) holds for every (x, u) ∈ B × U, we obtainf (x, u, v) = (v T v -1) h(x 0 , u 0 , v), ∀(x, u, v) ∈ B × U × R m ,and as B × U × R m has nonempty interior,

  2, . . . and lim

	d→∞	vol G d = vol P.
	In addition, we may want the sequence of inner approximations to satisfy additional
	nesting or convexity conditions.
	Problem 2 (Nested Inner Approximations) Solve Problem 1 with the additional
	constraint
	G d ⊆ G d+1 ⊆ P, d = 1, 2, . . .
	Problem 3 (Convex Inner Approximations) Given set P, build up a sequence of
	nested basic closed convex semialgebraic sets

  Therefore ŷ is a strictly feasible solution of (19) and so Slater's condition holds for (19).

Denote by Σ[x, u, v] d the space of polynomials of degree at most 2d, that are SOS of polynomials in R[x, u, v]/J. As ŷ is a strictly feasible solution of the semidefinite program (19), by a standard result of convex optimization, there is no duality gap between (19) and its dual

ρ ′ d = B λ min (x) dx -min g,

r,s,t B g(x) dx s.t. v T

A basic semialgebraic set is a set defined by intersecting a finite number of polynomial superlevel sets.

Note however that there is an incorrect factor 2 -n in the right handside of equation(3.3) in[START_REF] Lasserre | Solving a class of multivariate integration problems via Laplace techniques[END_REF].

V is Zariski dense in V C (J) (= {v ∈ C n : θ(v) = 0}) so that I(V R (J)) = I(V C (J)). But θ being irreducible, J is a prime ideal and so I(V C (J)) = J.

Acknowledgements

The first author acknowledges support by project No. 103/10/0628 of the Grant Agency of the Czech Republic. We are grateful to Luca Zaccarian for pointing out a mistake in a previous version of this paper.

Proof of Corollary 2

Proof: Let 0 < ε < 1 3 be fixed, arbitrary. As in the proof of Theorem 1, for every k ∈ N there exists a polynomial h k ∈ R[x] such that sup x∈B |λ min (x) -h k (x)| < ε k . Hence for all x ∈ B and all k ≥ 1,

and so the polynomial x → p k (x) := h k (x) -2ε k satisfies p k+1 (x) > p k (x) and λ min (x) > p k (x) for all x ∈ B. Again, by Putinar's Positivstellensatz, see e.g [11, Section 2.5], p k is feasible for [START_REF] Henrion | Positive polynomials and robust stabilization with fixed-order controllers[END_REF] with the additional constraint [START_REF] Nurges | Robust pole assignment via reflection coefficients of polynomials[END_REF], provided that d is sufficiently large, and with associated value

An auxiliary result for the proof of Lemma 1

Remember that J ⊂ R[v] is the ideal generated by 1 -v T v and the real radical

for some polynomial


P Baillieux
email: paul.baillieux@unine.ch

E Schill

L Moresi

Y Abdelfettah

C Dezayes

INVESTIGATION OF NATURAL PERMEABILITY IN GRABEN SYSTEMS: SOULTZ EGS SITE (FRANCE)

Two different approaches have been employed to investigate the localization of fault zones and natural permeability at regional and local levels in extensional tectonic settings with potential for EGS technology.

Regional distribution of potential permeable fault zones has been analyzed dynamically using a new 2D geodynamic model. Numerical simulation of instability of a viscoelastic-plastic lithosphere require a lower crustal viscosity in the order of 5•10 21 to 10 22 Pa s to match the deformation pattern of the Upper Rhine Graben (URG). Models clearly show a localization of deformation at regional level, but do not explain the internal fault distribution on the horst level as it can be seen at the Soultz. This has been investigated further using geophysical means.

The high-resolution 3D structural regional model of the Soultz area is used as a basis for new geophysical interpretation and inversion of gravimetric and magnetic data. Several gravimetric and magnetic anomalies can be linked to thermal anomalies, which in general the Upper Rhine Graben are known to be linked to free convection and thus, linked to porosity and permeability changes.

INTRODUCTION

The URG reveals several heat flow anomalies in the order of up to > 160 mW m -2 . One of those at Soultz-sous-Forêts has been selected to develop the Hot Dry Rock technology in 1987. Given the occurrence of natural fractures and hydrothermal fluid in the subsurface, the project has been developed using Enhanced Geothermal System (EGS) technology and is today the only EGS project under exploitation world-wide. The enhanced surface heat flux in the Upper Rhine valley has been related to the circulation of thermal water along fault zones [START_REF] Illies | Bauplan und Baugeschichte des Oberrheingrabens. Ein Beitrag zum "Upper Mantle Project[END_REF]. Numerical simulation of the neighboring geothermal site of Landau have confirmed this hypothesis [START_REF] Bachler | Impact of grabenparallel faults on hydrothermal convection -Rhine Graben case study[END_REF].

Localization of geothermal highs in the URG can be observed at two different levels. At a region scale for example, these are found in the North Western part of the graben (Fig. 1), in the part which is opposed to the major deformation accommodated by the Eastern master fault [START_REF] Mauthe | Kohlenwasserstoffvorkommen und -potential im deutschen Teil des Oberrheingrabens[END_REF]. At a local scale, the example of the Enhanced Geothermal System (EGS) site Soultz is found to be located in the northern part of a horst structure. Within this horst, the thermal anomaly localizes between the Soultz and Kutzenhausen faults only at its western edge [START_REF] Baillieux | 3D structural regional model of the EGS Soultz site (northern Upper Rhine Graben, France): insights and perspectives[END_REF].

Other recent studies highlighted the link between temperature anomaly and its state of fracturing at Soultz: a thermo-hydraulic simulation around the geothermal wells predicts fluid convection with permeability up to 3x10 -14 m 2 in these strongly fractured zones [START_REF] Kohl | Steps towards a comprehensive thermo-hydraulic analysis of the HDR test site Soultz-sous-Forêts[END_REF]. The low resistivity anomaly at the western part of the horst in the magnetotelluric inversion is interpreted to be the presence of clay minerals, due to hydrothermal alteration of granite or by the presence of the saline brine by itself [START_REF] Geiermann | 2-D Magnetotellurics at the geothermal site at Soultzsous-Forets: Resistivity distribution to about 3000 m depth[END_REF]. A density anomaly of about 250 kg.m -3 has been observed in the granitic basement below the Soultz site on the 3D inversion of existing gravity data in combination with a previous geological model [START_REF] Schill | 2-D Magnetotellurics and gravity at the geothermal site at Soultz-sous-Forêts[END_REF]. It has been attributed to an increase in porosity due to strong fracturation and the presence of leucogranitic bodies in this structure, supported by geological observations in the wells EPS1 and GPK1 [START_REF] Genter | Fracture analysis and reservoir characterization of the granitic basement in the HDR Soultz project (France)[END_REF].

In this study, we aim at understanding the formation and localization of fracturing in graben systems using an integrated approach:

First, we use graben-wide seismic sections [START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF]) (see Fig. 2 and their location in Fig. 1) and geological knowledge of the well-studied URG as an ideal database for benchmarking geodynamic models of crustal extension at a regional level, allowing the visualization of formation of faulted zones through time.

Secondly, our high-resolution 3D structural regional model [START_REF] Baillieux | 3D structural regional model of the EGS Soultz site (northern Upper Rhine Graben, France): insights and perspectives[END_REF]) (location in Fig. 1 ) is used as a basis for new geophysical interpretation and inversion of gravimetric and magnetic data. These are used for re-interpretation of the geophysical anomalies.

KEY TECTONIC INFORMATION ON THE UPPER RHINE GRABEN

The Upper Rhine Graben (URG) is a central segment of the European Cenozoic Rift System (ECRIS) and originates from the Alpine and Pyrenean collisions [START_REF] Ziegler | European Cenozoic rift system[END_REF]. It represents a typical example of synorogenic, intra-continental foreland rifting [START_REF] Schumacher | Upper Rhine Graben: Role of preexisting structures during rift evolution[END_REF][START_REF] Dèzes | Evolution of the European Cenozoic Rift System: interaction of the Alpine and Pyrenean orogens with their foreland lithosphere[END_REF][START_REF] Cloetingh | Lithosphere tectonics and thermo-mechanical properties: An integrated modelling approach for Enhanced Geothermal Systems exploration in Europe[END_REF]. This rifting is related to foreland indenter effects and ensuing escape tectonics: lithospheric over-thickening in orogenic belts, resulting in uplift and extension of their axial part. The main rifting stage of the URG formation was an E-W extension taking place from the end of the Oligocene to that of the Eocene [START_REF] Villemin | L'évolution structurale du fossé rhénan au cours du Cénozoïque : un bilan de la déformation et des effets thermiques de l'extension[END_REF]. Inherited Variscan tectonic structures have a strong impact on the localization of deformation up to recent times [START_REF] Schumacher | Upper Rhine Graben: Role of preexisting structures during rift evolution[END_REF]Edel et al. 2006). [START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF].

Recent thermo-tectono-stratigraphic forward modeling of the URG [START_REF] Hinsken | Thermo-Tectono-Stratigraphic Forward Modeling of the Upper Rhine Graben in reference to geometric balancing: Brittle crustal extension on a highly viscous mantle[END_REF] suggests plane strain deformation perpendicular to the graben during the Middle Eocene to Early Miocene, Late Miocene post-rift stage and renewed rifting during the Pliocene to recent Rifting occurred at very low strain rate of about 1.7x10 -16 s -1 involving brittleelastic deformation of the crust and ductile deformation of the highly viscous, high strength upper mantle and a "pre-rift" necking depth of 29 km coinciding with the Moho. This is in agreement with previous results indicating a factor of extension of on average 1.15-1.2 for the whole rift and a total amount of E-W extension of 6-7 km during approximately 40 Ma [START_REF] Villemin | The Rhinegraben: Extension, subsidence and shoulder uplift[END_REF][START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF].

Today the URG boundary faults are thought to operate in a left lateral strike slip sense under the NW-SE oriented compression [START_REF] Illies | Holocene movements and state of stress in the rhinegraben rift system[END_REF][START_REF] Plenefisch | The stress field in the Rhine Graben area inferred from earthquake focal mechanisms and estimation of frictional parameters[END_REF][START_REF] Lopes Cardozo | Kinematic analysis of the Upper Rhine Graben boundary fault system[END_REF] observed over much of Western Europe.

The URG is characterized by different asymmetries in the distribution of deformation (Fig. 1 and Fig. 2): 1. In the southern part, the sediment deposition center is located on the western side of the Graben (Brun et al. 1991), where the vertical motion (around 3 km) was accommodated by the western border fault [START_REF] Brun | Crustal-scale structure of the southern Rhinegraben from ECORS-DEKORP seismic reflection data[END_REF][START_REF] Cornu | Numerical backward and forward modelling of the southern Upper Rhine Graben (France-Germany border): new insights on tectonic evolution of intracontinental rifts[END_REF] 2. In the northern part, the sediment deposition center is located on the eastern side of the graben where the vertical motion is maximum along the eastern border fault [START_REF] Wenzel | A deep seismic reflection profile across the northern rhine graben[END_REF]. 3. The master fault is interpreted to switch side in the southern transfert zone located at in the central part of the URG [START_REF] Mauthe | Kohlenwasserstoffvorkommen und -potential im deutschen Teil des Oberrheingrabens[END_REF]) (Fig. 1).

Graben wide seismic sections reveal fault zones down to a depth of 10-15 km (Fig. 2).

SHORT

REVIEW OF GEODYNAMIC MODELING OF CONTINENTAL EXTENSION

Recent geodynamic modeling aims at predicting observations in rifting environments such as the geometry and modes of crustal extension [START_REF] Huismans | Symmetric and asymmetric lithospheric extension: Relative effects of frictional-plastic and viscous strain softening[END_REF][START_REF] Huismans | Effect of plastic-viscous layering and strain softening on mode selection during lithospheric extension[END_REF][START_REF] Wijns | Mode of crustal extension determined by rheological layering[END_REF][START_REF] Buiter | Dissipation analysis as a guide to mode selection during crustal extension and implications for the styles of sedimentary basins[END_REF], the fault and strain distribution [START_REF] Lavier | Factors controlling normal fault offset in an ideal brittle layer[END_REF][START_REF] Regenauer-Lieb | Strain localisation and weakening of the lithosphere during extension[END_REF] and the distribution of continental strength in general. Major interests with respect to our study are:

1. What are the effects of extension and compression, and phenomena such as structural inheritance, basin inversion, and far-field stresses on the formation and evolution of sedimentary basins in extensional, compressional, and strike-slip settings? 2. How can we define the structural styles of faulting in the shallow brittle part of the lithosphere, the brittle-ductile transition, and the shear zone in the ductile part of the lithosphere?

On a global geological point of view rifting is firstly controlled by forces controlling the movement and interaction of lithospheric plates such as boundary stresses, slab pull, ridge push, collisional resistance, frictional forces caused mainly by the influence of the convecting mantle on the base of the lithosphere [START_REF] Cloetingh | Tectonic Models for the Evolution of Sedimentary Basins[END_REF]. [START_REF] Cloetingh | Tectonic Models for the Evolution of Sedimentary Basins[END_REF] argue, by looking at the best documented sedimentary basin systems in Europe (including the Upper Rhine graben), that the thermal and mechanical structure of the lithosphere is the main control on the development of sedimentary basins.

By considering direct evidence for fault strength, [START_REF] Scholz | Fault Mechanics[END_REF] conclude that stress in the lithosphere is limited by faulting and is determined, to first-order, by Byerlee"s law with hydrostatic pressure [START_REF] Watts | An Overview[END_REF]). Byerlee"s law describes the stress state in the Earth's upper crust at which fracturing along a geological fault takes place by solving the Mohr-Coulomb criterion:

τ = Co + μ (σ n -Pf)
In which τ is the shear stress and σ n the normal stress. Co is the cohesion or internal strength of the material and the value Pf is the pore fluid pressure inside the rock. The fact that pore pressure is hydrostatic in the crystalline rocks of the crust is a consequence of the high permeability of the crust resulting from the presence of faults [START_REF] Townend | How faulting keeps the crust strong[END_REF][START_REF] Scholz | Fault Mechanics[END_REF].

The geometry of rifts and sedimentary basins and their mode of extension are found to be linked to the amount of strain softening (the relation which describes the decrease of strength of a rock when subject to faulting) in the brittle upper crust, the thickness of the upper crust, the viscosity and strength of the lower crust, the extension rate and other processes linked with the extension such as gravitationally driven deformation, isostasy, magmatism and necking [START_REF] Buiter | Dissipation analysis as a guide to mode selection during crustal extension and implications for the styles of sedimentary basins[END_REF]) (and references therein). They also analyzed the dependence of modes of crustal extension of a brittle crust on the strength of a ductile lower crust, its effective viscosity, the extension rate, and the layer thicknesses. They showed that an asymmetric basin is more likely to develop for a strong brittle layer which has a high amount of strain softening, a weak viscous layer (lower crust) and slow extension.

I-2D CRUSTAL EXTENSION OF THE UPPER RHINE GRABEN AND IMPLICATION FOR GEOTHERMAL PROSPECTION

The aim of the study is to model the crustal extension of the Upper Rhine Graben satisfying the geophysical and geological observations in Fig. 2 and Fig. 3 with a localization of strain at one side of the Graben, in form of concentration of fault zones. We use the geodynamic platform Underworld. It uses a Lagrangian particle-in-cell finite element scheme (Moresi et al. 2007) that enables the accurate tracking of stress and strain-rate history in simulations involving large-scale deformation. This platform includes a toolkit for studying the geodynamics of a viscoelastic-plastic lithosphere in 3D. This methodogy has proven to be successful in simulating plausible scenarios for the creation of rifting environments such as the distributed faulting mode at the North Sea or metamorphic core complexes of the western U.S.A. and the Aegean [START_REF] Wijns | Mode of crustal extension determined by rheological layering[END_REF].

Sensitivity Study

We tested the effect of the initial settings on the faulting style during the extension, adapting input parameters from the model of [START_REF] Buiter | Dissipation analysis as a guide to mode selection during crustal extension and implications for the styles of sedimentary basins[END_REF] to values recorded in the URG (Figure 3 and Table 1). We used a 2D two-layer model with an upper and a lower crust layer with the same uniform density. The upper crust deforms in a rigid-plastic manner using a Drucker-Prager frictional-plastic pressure dependent law which is an attempt to produce a smooth yield surface with otherwise similar characteristics to the Mohr-Coulomb yield surface (Moresi et al. 2007). There is an air layer above the solid material to approximate a free surface and allow the formation of topography. The sidewall boundary condition is freeslip vertically, with velocity conditions applied to the left and right to drive extension. An initial 100% damaged singularity was set at the boundary between the lower and upper crust to have control on the localization of the deformation.

Effect of a variation in viscosity

For this exercise, three models were stretched to reach 10% of deformation (5km on each side). The initial 100% damaged singularity was chosen to be 2x2 km. The lower crust viscosity values are 10 21 , 5x10 21 and 10 22 Pa s are chosen in order to model the strength of an intermediately strong to a strong lower crust [START_REF] Wijns | Mode of crustal extension determined by rheological layering[END_REF]. * The chosen value of extension velocity corresponds to a total amount 7km of graben extension over a 100km model during 40Ma [START_REF] Villemin | The Rhinegraben: Extension, subsidence and shoulder uplift[END_REF][START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF]. This corresponds to a slow extension process [START_REF] Hinsken | Thermo-Tectono-Stratigraphic Forward Modeling of the Upper Rhine Graben in reference to geometric balancing: Brittle crustal extension on a highly viscous mantle[END_REF].

The results are shown in Figure 4, represented in the final Lagrangian particle swarm distribution and cumulative strain rate distribution after 10 km of extension when using a lower crust viscosity of 10 21 Pa s. As expected the singularity leads the deformation. A single fault crossing the entire upper crust appears during the first kilometers of extension when the lower crust strength is intermediate (Figure 4A) to intermediately high (Figure 4B). The deformation then starts to be distributed along a generation of normal and conjugate faults which forms the landscape of an asymmetric basin. The modeling of continental extension with an intermediately strong lower crust (Figure 4B) shows the apparition of an extra border fault on the other side of the graben and a relatively equal distribution of conjugate normal faults between the two border [START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF] but not yet clearly observed. Here, the effect of an initial singularity is probably too strong to see how the faults "naturally" develop. In the case where the lower crust is strong (Figure 4C), a symmetric basin tends to develop.

Effect of the initial geometry

In this part, the initial singularity was removed and replaced by a random distribution of damage within the envelope function φ=sin 2 (k1. x1), where k1=1/L1 and L1 is the length of the model, allowing the damage to be kept away from the boundaries (Figure 5). Three models were stretched to reach 15% of deformation as observed in the Upper Rhine graben (7.5 km on each side).

In the case of the low value of viscosity for the lower crust (i.e. intermediately strong lower crust), numerous faults develop in the early stages of the deformation in all model but they no longer accumulate slip after around 5% of extension and only one main fault take control over the deformation, stretching from the lower crust to the upper crust (Figure 6A). In the case of a stronger lower crust (i.e. viscosity values of 5x10 21 and 10 22 Pa s), in the two other models, this main fault is followed by the generation of conjugate faults at the surface. One of them (after 10% of extension) eventually grows and connects to the first main fault at the border between the two crusts, forming the second border fault of the basin (Figure 6B andC) Faults that are outside the basin tend to no longer accumulate slip, whereas the one that are inside the basin keep their shape and form the landscape of a graben and horst system mainly located on the side of the second border fault.

Another model with a lower crust viscosity value of 10 23 Pa s was run: the resulting model shows the formation of multiple basins as we could expect for this high value [START_REF] Buiter | Dissipation analysis as a guide to mode selection during crustal extension and implications for the styles of sedimentary basins[END_REF]). This result is not comparable to what we observe in the URG.

Comparison with independent geophysical data

The results of our geodynamic modeling find similarities with both the ECORS-DEKORP crustalscale seismic exploration of the URG and the recent 3D geological model of the Soultz area, located at the western side in the central part of the URG:

For an intermediately strong to a strong lower crust, we find that after development of symmetrically conjugated faults, a preferential orientation develops into one main fault crossing the entire upper crust, which accommodates the deformation in this phase of graben formation. The presence of such a fault is suggested in the crustal-scale seismic exploration of the Rhine Graben [START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF]). In the following phase a second boundary fault appears and the minor fault density is increasing in the area next to that boundary fault. This second boundary fault is also observed in the seismic results and provides indication for an asymmetric distribution of fault zones within a larger tectonic unit. This may explain also the fault density and temperature distribution in the area of Soultz.

Discussion and implication for geothermal exploration

Although we used simple assumptions (constant extension rate, no isostatic compensation, no shear heating, no temperature profile, no sediment deposition, no elasticity) to model the crustal extension, we find first order results that can be compared to geophysical data. For lower crust viscosity values of 5•10 21 to 10 22 Pa s, we approached the patterns of deformation observed in the Upper Rhine Graben.

These first order results have an implication for the geothermal prospection of the Cenozoic Europen Rift System. For example the supposed masterfaults could be potential for highly hydraulically conductive rocks as well as the region of horst and graben such as the Soultz area on the other side because of their long period of activity and renewed activity of the URG up to Recent.

II-GEOPHYSICAL RE-INTERPRETATION OF THE AREA OF SOULTZ

As temperature distributions, also other geophysical parameters reveal a number of anomalies throughout the investigation area. The focus of this study is the distribution of different anomalies with relation to the crystalline basement. To analyze the contribution of the basement to the temperature distribution, temperature data from the local 110 boreholes [START_REF] Pribnow | Thermal tracking of upper crustal fluid flow in the Rhine Graben[END_REF] were interpolated in 3D and plotted on the basement top (Fig. 7). Temperature highs are found to be not only located at the north-western part of the horst of Soultz but also about 7 km away in a south-east direction. In this area the top of the basement occurs at significant depth compared to Soultz. The distribution of faults however, shows a similar pattern. In order to investigate this apparent link between the structural pattern and the temperature distribution we used the high-resolution 3D structural regional model of the Soultz area (30x20x6km) (Fig. 8 and location in Fig. 1) as a basis for new geophysical reinterpretation and inversion of gravimetric and magnetic data anomalies. (red). After [START_REF] Baillieux | 3D structural regional model of the EGS Soultz site (northern Upper Rhine Graben, France): insights and perspectives[END_REF].

To avoid physical inconsistency only a selection of the existing data from [START_REF] Rotstein | Insight into the structure of the Upper Rhine Graben and its basement from a new compilation of Bouguer Gravity[END_REF] has been used in this study. The Bouguer anomaly (Fig. 9) shows an important gravity low in the region of Soultz, which has been attributed to a low density granitic and granodioritic basement of Lower Carboniferous age [START_REF] Edel | Structure et évolution du Fossé Rhénan, du Carbonifère à nos jours -apports de la géophysique[END_REF]). It occurs as several types of granite and sub-facies in the geothermal wells [START_REF] Genter | Géothermie roches chaudes sèches: le granite de Soultz-sous-Forêts[END_REF][START_REF] Hooijkaas | Deepseated geology of the granite intrusions at the Soultz EGS site based on data from 5 km-deep boreholes[END_REF]) and reveals a characteristic density of 2630 kg m -3 [START_REF] Rotstein | Insight into the structure of the Upper Rhine Graben and its basement from a new compilation of Bouguer Gravity[END_REF] with local minima of down to 2460 kg m -3 in the hydrothermally altered zones with a high clay content and a porosity of up to 20% [START_REF] Genter | Géothermie roches chaudes sèches: le granite de Soultz-sous-Forêts[END_REF]). Measurement on rock samples from the deeper levels resulted in a density of 2520 kg m -3 for the two-micagranite at 4909-4973 m and 2530 kg m -3 at 5018-5036 m for the biotite-rich granite in GPK3 borehole [START_REF] Grecksch | Thermophysical Study of GPK2 and GPK3 Granite Samples[END_REF]. In the central URG, the magnetic anomalies are attributed to pre-Permian basement bodies with different magnetic susceptibilities [START_REF] Lauer | Propriétés magnétiques des roches au voisinage du contact métamorphique des schistes de Steige et des granites d"Andlau et du Hohwald (Vosges cristallines du Nord)[END_REF][START_REF] Edel | Paleomagnetic investigations and K-AR dating on the variscan plutonic massif of the champ du feu and its volcanic-sedimentary environment, northern vosges, France[END_REF][START_REF] Edel | The Upper Rhenish Shield Basement (Vosges, Upper Rhinegraben and Schwarzwald) -Main Structural Features Deduced from Magnetic, Gravimetric and Geological Data[END_REF]. In the area of Soultz (Fig. 10), at the south-eastern side of the geothermal wells the 15 km-large magnetic structure has been estimated from surface magnetic measurements to be between 1.2 km and 2 km depth using half-slope methods [START_REF] Papillon | Traitements et interpretations des cartes d'anomalies magnétiques et gravimétriques du Fossé Rhénan supérieur[END_REF], which corresponds to the depth of the top basement and to the measurements in the Soultz well. For example, in EPS1 borehole, the large increase in magnetic susceptibility in the granite was measured at 1550 km bsl [START_REF] Rummel | Density, ultrasonic velocities and magnetic susceptibility measurements on the core material from borehole EPS1 at Soultz-sous-Forêts[END_REF] on granitic samples, below the Variscan paleo-surface. In other parts of the graben magnetic structures globally coincide as well with the lower Triassic depth, i.e. the top of the basement [START_REF] Papillon | Traitements et interpretations des cartes d'anomalies magnétiques et gravimétriques du Fossé Rhénan supérieur[END_REF]. Magnetic data of the URG and rock susceptibility measurements from the same areas showed generally that the magnetic bodies in the URG mostly correspond to Carboniferous gabbros, diorites and granodiorites, and their volcanic equivalents (Edel et al. 2006).

In the following, we compared the magnetic anomalies reduced to the pole (Fig. 10) with the tectonic features of our 3D geological model. In the central part of the model, it can be observed that the fault system follows the 15 km-large magnetic structure. Thus, this magnetic structure seems to have generated the localization of deformation during the graben opening.

Inversion of the gravity data

Two 3D forward models with varying densities for the basement were calculated using GraviFor3D (Abdelfettah and Schill subm.) to investigate the gravimetric response of the geological model. Densities were chosen using values estimated in the URG [START_REF] Plaumann | Erweiterungen Kompakter Gruppen Durch Abelsche Gruppen[END_REF][START_REF] Rousset | Structure of the Southern Rhine Graben from Gravity and Reflection Seismic Data (Ecors-Dekorp Program)[END_REF][START_REF] Rotstein | Insight into the structure of the Upper Rhine Graben and its basement from a new compilation of Bouguer Gravity[END_REF]) and measurements in the Soultz boreholes [START_REF] Genter | Géothermie roches chaudes sèches: le granite de Soultz-sous-Forêts[END_REF]) as well as from core samples [START_REF] Rummel | Density, ultrasonic velocities and magnetic susceptibility measurements on the core material from borehole EPS1 at Soultz-sous-Forêts[END_REF]. The first forward model with a basement density of 2600 kg m -3 fits the dynamic of the measurements (22 mgal in the forward model compared to 20 mgal in the Bouguer), but it is not in line with the observed anomalies (Fig. 11). The second forward model with a basement density of 2500 kg m -3 shows a dynamic by 12 mgal lower than the measured Bouguer anomalies (Fig. 12).

To match the dynamic of the forward modeling, we have applied a Butterworth 30 km high-pass filter to the measured Bouguer anomaly (Fig. 13). It results in a good fit in dynamic and density distribution with the results from the forward modeling with a uniform basement density of 2500 kg m -3 . Following this results we used the residual anomaly in Fig. 13 for inversion to determine densities that would fit the measurements. Gravity inversion can be considered as a linear or nonlinear approach. In the first case, the objective is to find (invert) only the density values and geometries are considered as known. In the second case, the geometries and the densities are simultaneously inverted. The challenge of the inversion, especially for the potential method, is the non-uniqueness of the solution caused by the fact that the data number is always less than the parameters number used in the inversion run. For all these considerations, we chose to use the constrained 3D inversion in which the 3D geological model is fixed.

The minimization algorithm used is the routine developed by [START_REF] Beiner | Fortran routine mindef for function minimization[END_REF]. Principle of the inversion process is to minimize the misfit between the observed gravity data and the gravity response of the geological model. At each iteration, the forward modeling is achieved and its response is compared with real data and the discrepancy is minimized. As long as the misfit is great, the inversion parameters, the density values, are modified until that the minimization conditions are satisfied. This approach can be assimilated to a topographical surface in which we search the lowest elevation point [START_REF] Fischer | Topography and Minimization of the Standard-Deviation in One-Dimensional Magnetotelluric Modeling[END_REF]. The misfit function F is defined by

           M i i com i obs i g g = d F 1 2 ) (  (1)
where d is the density or density contrast values, σ is the standard deviation of the observed data, g obs and g com are the observed and computed gravity anomalies, respectively. M is the total gravity measurements. The smoothness term of form

      M i M i j j d i d = d S 1 1 2 )] ( ) ( [ ) (
(2) can be added to the minimization defined in (1). This smoothness term is used to avoid a huge density contrast which may not be natural.

The results of the inversion are the densities determined in Table 3. The density values match the ones found in the literature. The basement density is found to be lower than expected. The gravimetric response of the model computed after inversion (Fig. 14) shows strong similarities with the residual anomaly, but also areas with important differences of up to 7 mgal (Fig. 15).

Interpretation

We find that the 3D geological model can explain most of the residual anomalies, if we use a low density value for the basement. The low density for the basement could be explained by a strong natural fracturing or/and the low density of the two-mica and biotite-rich granites that are found at depth in Soultz boreholes. The residuals also show zones of density lows that cannot be explained by the 3D model as demonstrated by the misfit. Interesting to note is that about 1 mgal of misfit is observed in the area where the temperature distribution reveals a positive anomaly. As observed in Soultz such temperature anomalies go along with porosity anomalies. For example, for the reservoir zones of Soultz, we observe a mean porosity of about 8%.

Within the horst structure of the Soultz site, a density difference can be observed between the eastern and the western part. In the West, the residual anomaly is slightly smaller compared to the East. This, combined with the magnetic anomaly present at that location gives us interesting information: the Kutzenhausen-Soultz fault zone is a major structure in the region and has probably a large influence on the state of fracturing at Soultz.

CONCLUSION

The geodynamic approach we used provides insights in the localization of permeable structures at a regional level, at a first order comparable to the patterns of deformation observed in the Upper Rhine Graben.

The processing and inversion of gravity, the magnetic and temperature measurements and the 3D geological model reveal certain correlations at a local level, such as distribution of low density areas correlated to geothermal positive anomalies. The basement was found to be of a lower density than expected. This is in favor of a naturally porosity or low density basement lithology.

Figure 1 :

 1 Figure 1: Structural interpretation of the Upper Rhine Graben. Red lines denote the spatial distribution of the graben master faults. Blue dashed lines show the location of the ECORS-DEKORP seismic lines. Red stars show the main geothermal anomalies of the upper Rhine Graben.Green rectangle shows the geological model location. Modified after Derer et al. (2005).

Figure 2 (

 2 Figure 2 (previous page): Interpreted cross-sections of the Rhine Graben from the ECORS-DEKORP seismic investigation of the crustal structure of the Upper Rhine Graben. a) Northern profile. b) Southern profile. Red lines denote the interpreted master fault accommodating the maximum offset. Modified after[START_REF] Brun | Deep crustal structure of the Rhine Graben from seismic reflection data: A summary[END_REF].

Figure 3 :

 3 Figure 3: Model setup with plastic upper crust (Drucker-Prager frictional plastic pressure dependent law) and with viscous lower crust and an initial singularity

Figure 4 (

 4 Figure 4(left): Final Lagrangian particle swarm distribution (in km) (upper picture) and cumulative strain rate distribution (s-1) (lower picture) after 10 km of extension using a lower crust viscosity of A) 10 21 Pa s, B) 5•10 21 Pa s, C) 10 22 Pa s

Figure 5 :

 5 Figure 5: Model setup. Grey dots represent the distribution (5%) of initially damaged material.

  Figure 6: Final Lagrangian particle swarm distribution (in km) (upper picture) and cumulative strain rate distribution (s-1) (lower picture) after 15 km of extension using a lower crust viscosity of A) 10 21 Pa s, B) 5•10 21 Pa s, C) 10 22 Pa s.

Figure 7 :

 7 Figure 7: Temperature distribution on basement top. Fault and horizons traces are located on surface. Coordinates are in decimal degrees (UTM32N).

Figure 8 :

 8 Figure 8: 3D structural model of the Soultz area including the basement (grey), the Triassic-Jurassic layers (pink, purple, light and dark blue) and the fault planes(red). After[START_REF] Baillieux | 3D structural regional model of the EGS Soultz site (northern Upper Rhine Graben, France): insights and perspectives[END_REF]).

Figure 9 :

 9 Figure 9: Bouguer Anomalies (in mgal) in the Soultz area. Coordinates are in decimal degrees (UTM32N).

Figure 10 :

 10 Figure 10: Magnetic anomalies reduced to the pole in the Soultz area (in mgal). Coordinates are in decimal degrees (UTM32N).

FormationsFigure 11 :

 11 Figure 11: Forward modeling of the Bouguer anomaly at Soultz area with a uniform basement density of 2600 kg m -3 . Coordinates are in decimal degrees (UTM32N).

Figure 12 :

 12 Figure 12: Forward modeling of the Bouguer anomaly at Soultz area with a uniform basement density of 2500 kg m -3 . Coordinates are in decimal degrees (UTM32N).

Figure 13 :

 13 Figure 13: Residuals of the Bouguer Anomalies after Butterworth 30 km high-pass filter (in mgal) in the Soultz area. Coordinates are in decimal degrees (UTM32N).

Figure 14 :

 14 Figure 14: Gravimetric response of the model computed after inversion (in mgal) in the Soultz area. Coordinates are in decimal degrees (UTM32N).

Figure 15 :

 15 Figure 15: Difference between the gravimetric response of the model computed after inversion and the residuals (in mgal). Coordinates are in decimal degrees (UTM32N).

Table 1 :

 1 Rheological and modeling parameters used during the presented crustal extension modeling.

	Parameter	Natural Value Unit
	Thicknesses
	Crust	35	km
	Brittle crust	17.5	km
	Viscous crust	17.5	km
	Width of the model	100	km
	Extension velocity (per	2,19x10 -2 *	cm.a -1
	side)
	Density
	Crust	2800	kg.m -3
	Angle of internal friction	30 -> 4	°
	Cohesion	20 -> 2	MPa
	Strain softening range	0-1	-
	Linear viscosity (η)	10 21 to 10 22	Pa.s

Table 2 :

 2 Densities of formations chosen for a forward modeling of the Bouguer anomaly at Soultz area.

Table 3 :

 3 Densities determined by inversion.

	Formations	Initial density	After inversion
		(g.cm -3 )	(g.cm -3 )
	Tertiary	2.35	2.43 (+3.4%)
	Jurassic	2.55	2.61 (+2.35%)
	Keuper	2.7	2.58 (-4.4%)
	Muschelkalk	2.7	2.64 (-2.2%)
	Buntsandstein	2.5	2.66 (+6.4%)
	Basement	2.5	2.41 (-3.6%)

ACKNOWLEDGEMENTS

We kindly thank the "Donation Fund Bureau of Neuchâtel University" in for supporting the research project carried out at Monash University Australia, the (geo)dynamic team of the Mathematical Department of Monash University and Dr. Gabriele Morra for their help and support during the geodynamic modeling. We kindly thank Prof. J.-B. Edel for input geophysical data and fruitful discussions, and LIAG for providing input temperature data.


Chrystel Dezayes
email: c.dezayes@brgm.fr

Laurent Beccaletto

Laure Capar

Pauline Chenin

Identification of the fractured basement/cover unconformity on former oil seismic profiles in the Upper Rhine Graben

On the French territory of the Upper Rhine Graben, about 6.800 kilometres of seismic reflection have been acquired since 1954 to the 80's for the petroleum exploration. Raw seismic data which were recovered from oil and seismic companies were retreated by modern processing techniques. However, as the acquisition parameters of these old seismic profiles aimed to better image shallow petroleum reservoir from Tertiary to Jurassic sediments, the deepest levels were not imaged properly.

In the Rhine Graben, the base of the basin is constituted by lower Triassic Buntsandstein and locally by Permian sediments. These both levels are clastic sediments yielded by the degradation of the hercynian massif and are constituted by fluvial deposits in an arid climatic environment (red sandstones) [START_REF] Bourquin | The Permian-Triassic boundary and lower Triassic sedimentation in western European basins: an overview[END_REF]. These levels could not be distinguished from each other on seismic profiles because they have similar origin and composition. They overlain the Palaeozoic basement constituted by different granites, gneiss and metamorphic rocks, volcanic and volcano-sedimentary rocks [START_REF] Edel | The upper Rhenish Shield basement (Vosges, Upper Rhinegaben and Schwarzwald): main structural features deduced from magnetic, gravimetric and geological data[END_REF]. Some seismic profiles, acquired in 70's and 80's, have been recently retreated to improve the seismic images, particularly the continuity of the horizontal reflectors, faults and deep zones [START_REF] Capar | Apport du retraitement sismique à l'étude des bassins sédimentaires: exemple du nord du fossé rhénan[END_REF]. On these seismic profiles, the basement/cover unconformity has been mapped systematically in the northern part of the Rhine Graben, from Wissembourg to Selestat. However, few boreholes are deep enough to identify with certainty the reflector corresponding to the basement/cover unconformity. Moreover, , we tried to find a relation between the signature of the seismic reflector of the basement/cover unconformity and the nature of the underlying basement. It appears low magnetic rocks, such as granites, or Palaeozoic sediments, show a more continuous reflector than high magnetic rocks like gabbro, diorites, … Three recent seismic profiles have been retreated with a different signal process, the post stack time migration or PSTM, in order to improve the image of faults and their dip. The interpretation of these seismic profiles shows two types of faults, which cross-cut the basement /cover unconformity: 1-Major faults, which cross-cut the whole infill basin, except the Quaternary, and extend down-ward in the basement. Two different structural features have been identified: normal faults are bound to the early stage of the graben formation and negative flower structures, which reactive older fault systems in transtension during the recent tectonic stage.

2-Normal faults, which affect the top of the basement, the Buntsandstein and locally the Muschelkalk.

This study shows that the base of the basin is highly affected by large-scale fractures and constitutes a good target for geothermal projects. However, former petroleum seismic profiles, even if they are retreated, are not really appropriate to explore in detail the deepseated basin of the Rhine Graben (detail mapping of the top of basement, trend of the deepdeated fractures…). Thus, in the framework of deep geothermal exploration, dedicated acquisition of 2D or 3D seismic surveys would provide a better image of the underground and would contribute strongly to minimize the risks of future projects.

This study has been done in the framework of an ADEME-BRGM joint project, namely CLASTIQ 2. __________