Sylvain Pion
email: sylvain.pion@sophia.inria.fr

Andreas Fabri
email: andreas.fabri@geometryfactory.com

A Generic Lazy Evaluation Scheme for Exact Geometric Computations

Keywords: computational geometry, exact geometric computation, numerical robustness, interval arithmetic, lazy evaluation, generic programming, C++, Cgal

We present a generic C++ design to perform exact geometric computations efficiently using lazy evaluations. Exact geometric computations are critical for the robustness of geometric algorithms. Their efficiency is also important for many applications, hence the need for delaying the costly exact computations at run time until they are actually needed, if at all. Our approach is generic and extensible in the sense that it is possible to make it a library that users can apply to their own geometric objects and primitives. It involves techniques such as generic functoradaptors, static and dynamic polymorphism, reference counting for the management of directed acyclic graphs, and exception handling for triggering exact computations when needed. It also relies on multi-precision arithmetic as well as interval arithmetic. We apply our approach to the whole geometry kernel of Cgal.

Introduction

Computational geometry can be regarded as the art of designing algorithms and data structures that organize geometric data and extract structure out of it. The convex hull of a set of points, its Delaunay triangulation (see Figure 1) and Voronoi diagram are basic examples of such algorithms and data structures. They are usually described in the literature using the simple Real-RAM model, which assumes in particular that operations on real numbers are performed exactly and in constant time. Geometric algorithms are unique in the sense that they mix combinatorial techniques with numerical computations. Indeed, they build structures akin to graphs, and they handle point coordinates. This dual nature is at the root of many practical non-robustness issues, because approximate numerical computations violate key geometric properties, which in turn interferes with the discrete logic of combinatorial algorithms. Indeed, the crucial exactness assumption is not satisfied when floating-point arithmetic is used to implement operations on real numbers. With the increasing use of digital sensors like scanners or digital cameras in numerous areas, the amount of geometric data that needs to be processed is growing tremendously, making geometric algorithms more relevant than ever. This also pushes for a careful, efficient and general treatment of the numerical issues. In particular, the development of the Cgal library [1], a large collection of geometric algorithms implemented in C++, strongly expresses such a need.

The exact geometric-computation paradigm investigated by Yap [START_REF] Yap | Robust geometric computation[END_REF] offers a solution with a wide practical applicability. This paradigm states that, in order to achieve robustness of an algorithm, it is enough to guarantee the exactness of the results of the geometric predicates that it evaluates. These predicates are the functions that take as input the numerical data and compute a discrete result, such as a Boolean value, that is used in the control flow of the algorithm to build its combinatorial output. A simple example of predicate is the lexicographic comparison of the coordinates of two points and a more interesting, but more involved example is the orientation predicate described below.

Guaranteeing the exactness of the predicates can be done using an exact arithmetic, which is implemented, for example, by multi-precision integers or rational numbers (also known as arbitrary precision arithmetic). Moreover, for efficiency reasons, tools for certified approximate computations based on floating-point arithmetic, like interval arithmetic, are also used as a first step in the evaluation of the predicates, since they allow amortizing the cost of exact arithmetic operations. Such techniques are known as arithmetic filtering, since the interval arithmetic step filters out the easy instances. Naturally, it is common for geometric algorithms to call predicates over intermediate geometric objects that have been built from the input numerical data. A simple example of geometric construction is the computation of the intersection point of two segments specified by their end points. In this case, the exact geometric-computation paradigm implies that the intermediate geometric objects also need to be known exactly.

An efficient approach to implement such exact geometric constructions is to perform exact computations in a lazy way at the level of geometric objects, rather than at the underlying arithmetic level. It consists in computing a certified approximation of the objects using interval arithmetic and memorizing the necessary information such that a representation using exact arithmetic can be computed on demand, that is, when the exact computation of a predicate will later require it due to insufficient precision of the initial approximation. This approach is mentioned in [START_REF] Pion | De la géométrie algorithmique au calcul géométrique[END_REF] and an implementation is described in [START_REF] Funke | Look -a lazy object-oriented kernel for geometric computation[END_REF]. However, this implementation does not use the generic-programming paradigm, although the approach is fundamentally general. This missing part is the contribution of our paper.

In this paper we describe a generic design using C++ that allows these filtering and lazy exact computation techniques to apply to a large number of geometric primitives. In particular, we consider the complete geometry-kernel of Cgal.

The paper is organized as follows. Section 2 gives some background information about arithmetic tools and specific C++ features that are necessary for the detailed understanding of the paper. Section 3 describes the context and motivation in geometric computing in more detail, as well as the basics of a generic geometry-kernel allowing the use of different arithmetics and what can be done at this level. Section 4 presents our main design, namely, how geometric predicates, constructions and objects are adapted, while Section 5 illustrates how our scheme can be applied to the user own geometric objects and primitives. We also provide some benchmark data in Section 6 that confirms the benefit of our design and implementation. Finally, we list a few open questions related to our design in Section 7 and conclude with ideas for future work.

Prerequisites

This section introduces some auxiliary notions that help understanding the core of the paper.

Variadic C++ Templates

We use the C++ programming language in the examples of code throughout this paper. However, in order to improve clarity, we consider two small deviations. First, given that we use templates as a means to implement generic programming, we improve readability by omitting some uses of the typename keyword that are normally formally required when referring to a nested type in dependant contexts.

Moreover and more importantly, for the sake of genericity, we use functions of arbitrary arity. This is implementable in C++ by providing overloads up to a fixed arity, but for the purpose of simplifying the exposition, we use variadic templates [START_REF] Gregor | Variadic templates for C++[END_REF]. This feature, which will most likely become part of C++0x, the next revision of the C++ standard, allows a template to take an unspecified number of parameters and to apply expressions to them repeatedly using an ellipsis "..." as the syntax. In order to avoid confusion, we use additional spaces ". . ." for normal (non-variadic) ellipses. The following example: template < typename... Args > int fog(Args ... args) { return f(g(args)...); } is equivalent to the following set of overloads for any arity:

int fog() { return f(); } template < typename A1 > int fog(A1 a1) { return f(g(a1)); } template < typename A1, typename A2 > int fog(A1 a1, A2 a2) { return f(g(a1), g(a2)); } . . .
Let us also introduce some related tools. First, the variadic template class tuple stores any number of objects of any types. Then, the apply function calls a functor to each element of a tuple, and returns the results in a tuple of the appropriate type. We do not specify here how to obtain this type, it is named X... in the following example. Here is an example putting all these tools together, where apply is called with a functor that doubles its argument, and bind tuple is applied to a functor performing additions.

Arithmetic Tools

We describe now the computer arithmetics that we use: floating-point arithmetic, interval arithmetic and multi-precision arithmetic.

Floating-Point Arithmetic

Floating-point arithmetic is the usual way to compute with approximations of real numbers on a computer. The Ieee-754 standard defines all the details of its behavior, and its correct use is described in a well known article by Goldberg [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF]. We only recall the properties that matter here. The operations considered are the addition, subtraction, multiplication, division and square root. Floating-point values are mostly finite real numbers, but +∞ and -∞ are also representable.

First, arithmetic operations are necessarily approximate, which means that exactness is generally lost. This process is called directed rounding, as the result of operations can be seen as if they were performed exactly, but then rounded to one of the nearest floating-point values enclosing the exact value. By default, the nearest value is selected, this is rounding to nearest. The other relevant modes are rounding towards plus and minus infinity, which respectively select the representable floating-point value that is immediately smaller or larger than the exact result.

The rounding error, that is, the distance between the exact result of an operation and its floating-point equivalent, can be bounded, since it depends on the distance between consecutive floating-point values, which is specified by the standard. Its value is 0.5 × ulp(x) where x is the result, and ulp is the unit in the last place function, which is the difference between x and the next floating-point value. In double precision, ulp(1) equals 2 -53 . This property is interesting, but it exhibits a major problem when it comes to computation, which means composing operations. Indeed, considering the expression a + b + c, it is not possible to bound the rounding error as easily in terms of the floating-point result. Simply consider what happens when a = 1, b = 2 -60 and c = -1 to observe what is referred to as a catastrophic cancellation.

As a consequence, this model of approximate computation with real functions is hard to deal with when it comes to controlling the rounding errors. Instead, we present interval arithmetic, a tool adapted to solve this problem.

Interval Arithmetic

Interval arithmetic can be seen as an alternative to the floating-point model, but it can also be implemented easily on top of it. This concept has been introduced by Moore [START_REF] Moore | Interval Analysis[END_REF] as the basis for the field of interval analysis. The idea is simple; an interval is used to represent the rounding error attached to a variable by guaranteeing that it contains its exact value. This key property is named the containment property and the arithmetic operations are carefully defined to propagate it. In practice, two floating-point numbers are used to encode the lower and upper bounds of the interval, and therefore an approximation of a real number x is represented by an interval [x.lower, x.upper] which contains x. Note that the lower bound can equate -∞ and the upper bound +∞.

The arithmetic operations are typically implemented by making use of the rounding modes towards +∞ and -∞. For example, the addition of two intervals x + y is computed as [x.lower+y.lower, x.upper+y.upper], where + and + respectively denote the floating-point addition with rounding towards -∞ and +∞. Similarly, the subtraction xy is computed as [x.lower-y.upper, x.upper-y.lower]. The multiplication follows the same principle except that its definition requires a few case distinctions and is therefore less straightforward. Other functions can be defined and implemented similarly.

Due to this outward rounding, we can intuitively say that the width of the intervals will tend to grow as operations are accumulated. This can lead to a complete ineffectivity of the method for some algorithms which perform lots of such accumulations, but fortunately this is not a problem for many practical other cases. In summary, it is possible to build an interval arithmetic on top of the floatingpoint one relatively simply and at a running time that is a few times slower.

Rounding Mode Changes

We now consider a remark on the implementation. Current hardware is advantageous to the floating-point model, since it implements the basic operations directly, while interval arithmetic requires both twice the number of registers (for storing the two bounds) and frequent changes of the rounding mode in general. Most CPUs provide a way to change the rounding direction as a special register that can be modified by specific instructions. These rounding mode changes slow down the interval operations substantially. Fortunately, it is possible to remove some of them by making a few observations. The first is that it is sometimes possible, like in the common case of the interval addition, to replace an operation rounded towards -∞ by one rounded towards +∞. This gives [-((-x.lower)-y.lower), x.upper+y.upper]. This allows the replacement of some rounding mode changes by some faster equivalent computations that can even sometimes be folded by the compiler. But the bulk of the potential improvements follows from the observation that consecutive interval operations can keep the rounding mode set towards +∞, so that a block of code using only +, -, × and ÷ can perform a single change of rounding mode at the beginning and one restoration at the end.

For the purpose of code illustration, we show a C++ class that performs this initialization in its constructor and the restoration in its destructor. We also show how to use it for defining the interval addition as an example. A final note on this issue: besides the parameterization by the types of the bounds, it is useful to consider parameterizing an interval class by a Boolean flag specifying if the rounding mode has been set to +∞ before executing the function. Indeed, this allows to skip the rounding mode change at each call by moving this responsibility higher up in the call tree.

Interval Comparisons

For reasons of generic programming, it is also useful to consider the definition of comparison operators over intervals, corresponding to the comparison over real numbers. The semantic of a comparison between intervals being true (resp. false) is that it is true (resp. false) for all possible real values chosen in both interval arguments, otherwise it is indeterminate. The latter case cannot be returned through a Boolean value, but can be returned through a more complex type holding this uncertainty information, or more simply, by throwing an exception.

Multi-Precision Arithmetic

Geometric data naturally has a numeric aspect since it models entities using approximations of real numbers. However, geometric computations also have a strong algebraic structure, and in order to compute algebraic expressions exactly, it is necessary to be able to compute without approximations. For this, we can use multi-precision arithmetic, which performs exact computations over the integers Z and the rational numbers Q (up to the available memory).

It is worth noting that the time complexity and memory size used by this model are typically here at least one to two orders of magnitude worse than for the approximate models described above. Libraries providing such functionality exist, for example the C library Gmp [START_REF] Gmp | The GNU Multiple Precision Arithmetic Library[END_REF].

3 Exact Geometric Computations and the Cgal Kernel

Exact Geometric Computations

As mentioned in the introduction, when programming geometric algorithms, we usually consider some low-level functions -called geometric predicatesthat take some numerical data as input in the form of geometric objects such as points, and return a Boolean or an enumerated value. This discrete return value is then used to direct the flow of the program to construct some higherlevel data structure such as a graph. Internally, these predicates typically perform comparisons of numerical values computed from the input. A classic example is the orientation predicate of three points in the plane, which returns an enumerated value that indicates whether the three points form a right turn, are collinear, or form a left turn (see Figure 2). Using Cartesian coordinates for the points, the orientation is computed as the sign (as a threevalued function: -1, 0, 1) of the following 3-dimensional determinant, which reduces to a 2-dimensional one:

1 1 1
p.x() q.x() r.x() p.y() q.y() r.y() = q.x() -p.x() r.x() -p.x() q.y() -p.y() r.y() -p.y()

Most predicates can be expressed as signs of polynomial expressions over the coordinates of the input points. Evaluating such a function with floatingpoint arithmetic is going to introduce rounding errors, which can lead to an approximate value having a sign that differs from the exact one. The impact of wrong signs on the geometric algorithms can be disastrous, as for example, it can break some invariants like planarity of a graph, or make the program loop forever or simply crash. Didactic examples of harmful consequences can be found in [START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF]. r q q q p p p r r Fig. 2. The orientation predicate of 3 points p, q, r in the plane: right turn, collinear and left turn.

Operations computing new geometric objects, like the point at the intersection of two non-parallel lines, the circumcenter of three non-collinear points, or the midpoint of two points, are called geometric constructions. We use the term geometric primitives when referring to either predicates or constructions.

In order to tackle the non-robustness issues, many solutions have been proposed. We focus here on the exact geometric-computation paradigm [START_REF] Yap | Robust geometric computation[END_REF], as it is a general solution. This paradigm states that in order to ensure the correct execution of the algorithms, it is enough that all predicates return exact results. Basically, this means that all comparisons of numerical values need to be performed exactly.

A natural way to perform the exact evaluation of predicates is to evaluate the numerical expressions using exact arithmetic. For example, since most computations are signs of polynomials, it is enough to use multi-precision rational arithmetic provided by libraries such as Gmp [START_REF] Gmp | The GNU Multiple Precision Arithmetic Library[END_REF]. Note that exact arithmetic is also available for all algebraic computations using libraries such as Core [START_REF] Karamcheti | The CORE Library Project[END_REF] or Leda [START_REF] Burnikel | The LEDA class real number[END_REF], which is useful when doing geometry over curved objects. This solution works well, but it tends to be very slow.

We next describe the implementation of the geometry kernel of Cgal, which provides a generic implementation of these methods using C++ templates. We use it later to describe the design of our lazy exact evaluation framework.

The Geometry Kernel of Cgal

The geometric algorithms of Cgal are parameterized by the geometry they apply to. The geometry takes the form of a kernel [START_REF] Hert | An adaptable and extensible geometry kernel[END_REF][START_REF] Brönnimann | 2D and 3D geometry kernel[END_REF] regrouping the types of the geometric objects such as points, segments and lines, as well as the geometric primitives operating on them, under the form of functors. The Cgal kernel provides over 100 predicates and 150 constructions. Hence, uniformity and genericity are crucial for both the interface and the maintenance.

Cgal provides several models of kernels. The basic families are the class templates Cartesian and Homogeneous, which are parameterized by the number type used to store the coordinates of the points and to compute arithmetic operations inside the primitives. They respectively differ in the Cartesian (such as (x, y) for a 2D point) or homogeneous (such as (hx, hy, hw), equivalent to the Cartesian (hx/hw, hy/hw)) representation of the points. Their interfaces are otherwise similar and look like the following: Variants named Simple cartesian and Simple homogeneous are also provided, and they differ from the non-simple variants only in that they do not offer reference counting for the geometric objects.

These basic template models allow to use floating-point, multi-precision rational or interval arithmetic for example. More precisely, Cgal provides a hierarchy of concepts for the number types that describe the requirements on types to be valid number types for instantiating some kernel primitives. For example, some primitive might require a division or a square root, while most only need additions, subtractions, multiplications and comparisons. We refer the reader to the Cgal documentation for more details on this particular topic.

The functors are implemented in the following way (here the return type of the predicate is a three-valued enumerated type): result_type operator()(Point p, Point q) const { return Point( (p.x() + q.x()) / 2, (p.y() + q.y()) / 2 ); } };

Conversions Between Objects of Different Kernels

As much as conversions between number types are useful, Cgal also provides tools to convert geometric objects between different kernels. We briefly present these below, as they are used later. A kernel converter is a functor whose function operator is overloaded for each object of the source kernel and returns the corresponding object of the target kernel. Such conversions may depend on the details of representation of the geometric objects, such as homogeneous versus Cartesian representation. Cgal provides such converters parameterized by converters between number types. Here is for example the converter between kernels of the Cartesian family: Cgal also provides a way to find out the type of a geometric object (say a 3D segment) in a given kernel, given its type in another kernel and this second kernel. This is in practice the return type of the function operator of the kernel converter described above. The current implementation works by specializing on all known kernel object types like K1::Point 2 or K1::Segment 3. A more extensible approach could be sought, but this is not the main point of this paper. Now that we have described all the tools, we can present the heart of our paper, the generic design of the lazy exact computation framework.

Design of the Lazy Exact Computation Framework

In this section, we describe the filtering and lazy exact computations techniques. We first devise the filtered predicates, then we propose a lazy exact number-type, which we then generalize to all geometric constructions. We also show how to put things together in geometry kernels.

Filtered Exact Predicates

Evaluating a filtered predicate means first calling the predicate instantiated with interval arithmetic [START_REF] Brönnimann | Interval arithmetic yields efficient dynamic filters for computational geometry[END_REF]. If this fails because some interval is not precise enough, the predicate is reevaluated, this time instantiated with an exact number-type. As all predicates of a Cgal kernel are functors, we can use the following Filtered predicate adaptor, which is itself parameterized by four functors: (i) the exact version of a predicate, which is using exact arithmetic such as rational numbers, (ii) the approximate version of the predicate, which is using interval arithmetic, (iii) a functor that converts the input geometric objects to their exact counterpart (say, points with rational coordinates), and (iv) a functor that converts the input geometric objects to their enclosing interval equivalent. The class Filtered kernel is then obtained from a kernel K by repeating the adaptation for each predicate of K. The geometric objects as well as the constructions remain unchanged.

Lazy Exact Numbers

In order to extend the filtering technique to exact constructions, we can implement a number type that evaluates computations using interval arithmetic. Moreover, it must remember the way it has been constructed by storing the history of operations as a directed acyclic graph (Dag) [START_REF] Benouamer | A lazy solution to imprecision in computational geometry[END_REF]. Figure 3 illustrates the history Dag of the expression

√ x + √ y -x + y + 2 √ xy.
When a comparison is performed on variables of this number type and the intervals overlap, the Dag is used to recompute the values of each node recursively with an exact multi-precision type, hence giving the exact result. Cgal provides such a lazy number-type called Lazy exact nt<NT> that is parameterized by the exact type used to perform the exact computations when needed (such as a rational number-type). This can be seen as a wrapper on top of its template parameter that delays the computations until they are needed, as they might not be needed at all.

We can compare the performance of this approach against a straightforward use of rational arithmetic. Such a test is described in detail in Section 6. In particular, a comparison of the first two rows of Table 1 indicates that the lazy exact number-type approach gains a factor of 8 to 10 in running time, while it causes a loss of a factor of 7 in memory use caused by the storage of the Dag. We next describe a way to recover from this loss and optimize the running time even further.

Optimizing by Regrouping Operations

There are mostly two aspects that can be optimized. First, a node of the Dag is created and the corresponding heap-memory allocated for each arithmetic operation. Secondly, as mentioned in Section 2.2.2.1, interval arithmetic benefits from the elimination of redundant rounding mode changes, which is made possible for consecutive operations. In both cases, regrouping operations in more complex functions would help. For example, if we consider the computation of the determinant in the orientation predicate, we see that seven operations are needed, but we could just as well use one function of arity six for the whole determinant. This would introduce only one Dag node and allow to perform the optimizations mentioned by amortizing the costs over more arithmetic operations.

These remarks have lead to a new scheme mentioned in [START_REF] Pion | De la géométrie algorithmique au calcul géométrique[END_REF], and the description of an implementation has also been proposed in [START_REF] Funke | Look -a lazy object-oriented kernel for geometric computation[END_REF]. The idea is to introduce a Dag at the geometric level, by considering geometric objects and primitives for the nodes, and not only numbers and arithmetic operations. The next section describes such an optimized setup. Our design differs from the one in [START_REF] Funke | Look -a lazy object-oriented kernel for geometric computation[END_REF] in that we followed the generic-programming paradigm and made extensive use of templates, to make it as easily extensible as possible.

Lazy Exact Geometric Objects

Performing lazy exact constructions means performing constructions with interval approximations, and storing the sequence of construction steps. When later a predicate applied to these approximations cannot return a result that is guaranteed to be correct, the sequence of construction steps is performed again, this time with an exact arithmetic. Now the predicate can be evaluated correctly.

The sequence of construction steps is stored in a Dag. Each node of the Dag stores (i) an approximation (at), (ii) a pointer to an exact result (et), (iii) a functor able to perform the same construction exactly (ec), that is, to construct *et, and (iv) reference-counted pointers to the lazy objects that were arguments to the function (args...). The out-degree of a node is therefore the arity of the function. The example shown in Figure 4 illustrates that lazy objects can be of the same type, without being the result of the same sequence of constructions. a, b, i, m, and p are all point-ish. Therefore, we use a handle class template, Lazy exact, storing a reference-counted pointer to a node of the Dag. The types of the nodes of the Dag are organized in a hierarchy using Construction base as abstract base class. The reference-counting mechanism follows a design similar to the one described in [START_REF] Kettner | Reference counting in library design -optionally and with unionfind optimization[END_REF].

We now explain some of the classes in Figure 5 in more detail.

Lazy exact has Lazy exact nt as subclass, which provides arithmetic operations. Note that this framework handles arithmetic and geometric objects in a unified way. For example, a distance between geometric objects yields a lazy exact-number, and a lazy exact-number can become the coordinate of a The subclass Construction is used for generic functions. This class template is parameterized by the types of (i) the approximate construction functor (AC), (ii) the exact construction functor (EC), (iii) the Lazy kernel class (LK) described in more detail in Section 4.6, which provides the required conversion functors, and (iv) the arguments of the function (Args...). The arguments, which may be of arbitrary types, are stored in the class together with the exact version of the function. In the case of lazy exact geometric-objects or lazy exact numbers, the arguments are handles as described above. The constructor stores the arguments. It then takes their approximations and calls the approximate version of the functor.

In case the exact version of the construction is needed, this gets computed in the update exact method. The exact versions of the arguments are fetched, which in turn may trigger their exact computation if they are not already computed and cached. From the exact object one computes an approximate object directly again, as this may improve the precision of the approximation by avoiding any accumulation of past rounding errors and this can help in further uses of this lazy object.

Finally, the Dag is pruned. As the nodes of the Dag are reference counted, some of them may be deallocated by the pruning. Most often, the types forming Args are lazy exact objects. For performance reasons their default constructors generate handles to shared static nodes of the Dag.

The other derived classes store the leaves of the Dag. For efficiency considerations, there is a general purpose leaf class, but also more specialized ones, for example for creating a lazy exact number from an int.

Lazy Exact Constructions

So far we have only explained how lazy constructions are stored, but not how new nodes of the Dag are generated. The following functor adaptor maps the basic generic functors encoding the geometric constructions to their corresponding lazy versions. The functor first tries to construct a new node of the Dag. If inside the approximate version of the construction an exception is thrown because of an interval arithmetic comparison failure, the exact version of the construction is called, and only a leaf node for the Dag is created.

Special Case Handling

The generic functor adaptor works out of the box for all functors that return geometric objects or numbers.

Functors returning objects that are not made lazy are easily handled exceptions. An example in Cgal is the functor that computes a bounding box with double coordinates around geometric objects, whose width is not required to be tight. As the intervals corresponding to the coordinates of the approximate geometric object are already 1-dimensional bounding boxes, we never have to resort to the exact geometric object. The functor adaptor is trivial.

Other functors return a polymorphic object. For example, the intersection of two segments may either be empty, a point or a segment. In order to store such a polymorphic result, Cgal offers a class Object capable of storing one object of any type. The problem we have to solve is that the lazy exact functor must not return a lazy exact Object, but instead must return an Object holding a lazy geometric object. This is solved by iterating over all Cgal kernel types, to try to cast, and if it works to construct the lazy geometric object and put it in an Object again. This functor appears as Lazy construct returning object in the next code example.

Less trivial cases are functors that pass results of a computation back to reference parameters, or that write into output iterators. They need a special functor as well as special Construction classes. It is not hard to write them, but the problem is that they must be dispatched by hand, as we have no means of introspection.

Lazy Exact Geometry Kernel

We are ready to put all pieces together, by defining a new kernel that takes an approximate and an exact kernels as template arguments. It can then be instantiated, for example, as:

typedef Lazy_kernel<Cartesian<Interval_nt>, Cartesian<Gmpq> > Lazy_cartesian_kernel;

typedef Lazy_kernel<Homogeneous<Interval_nt>, Homogeneous<Gmpz> > Lazy_homogeneous_kernel;

In the current implementation, we use the preprocessor to generate the typedefs from a list of types, and template specializations for dispatching the special cases. Approx converter simply fetches the stored approximate object, while Exact converter fetches the exact object, possibly triggering its computation.

Extensibility

We have to distinguish between different levels of extensibility.

When Cgal kernels get extended by geometric objects and constructions, changes are needed in the lazy construction framework if the new constructions have "new" interfaces, e.g., two output iterators, followed by two reference parameters to return a result. This would need a new node type for the Dag, a new functor, and hard wired dispatching in the lazy kernel. Otherwise, there is nothing special to do.

When the user wants to extend the lazy kernel with his own geometric objects and constructions, he first has to add them to the kernel that then gets into the lazy computation machinery, as described in [START_REF] Hert | An adaptable and extensible geometry kernel[END_REF]. Then, our statement in the previous paragraph applies.

The Circular kernel 2 class of Cgal [START_REF] De Castro | 2D circular geometry kernel[END_REF][START_REF] Emiris | Towards an open curved kernel[END_REF] is an example for both. It extends the basic Cgal kernels by adding geometric types such as circular arcs, and it provides additional predicates and constructions around them, like circle instersections. A lazy exact variant of it is also provided by applying the generic filtered tools described in this paper.

Benchmarks

We present the results of two benchmarks, which show the benefit of our techniques. Firstly, we use a simple artificial algorithm and then a more complex and useful algorithm, the computation of the 3D flow complex. The benchmarks compare the running time and the memory consumption of various kernel choices.

Here is the description of our artificial algorithm:

• generate 2000 pairs of 2D points with random coordinates (using drand48()). • construct 2000 segments out of these points.

• intersect all pairs of segments among these, and store the resulting intersection points. • shuffle the resulting points • iterate over consecutive triplets of these points, and compute the orientation predicate of these.

Table 1 provides the resulting data for a choice of four different kernels:

• SC<Gmpq> stands for the Simple cartesian representation kernel instantiated with Gmpq, which is a C++ wrapper around the multi-precision rational number type provided by Gmp, • SC<Lazy exact nt<Gmpq>> uses the lazy exact evaluation mechanism at the arithmetic level, • Lazy kernel<SC<Gmpq>> is our approach for performing lazy exact evaluations at the geometric object level, The results show that our approach wins almost a factor of 10 on memory over the basic lazy-evaluation scheme. It is also between 2 and 3 times faster. While it remains 4 times slower than the approximate floating-point evaluation, it is guaranteed for all cases.

The benchmarks also show the gain obtained thanks to the elimination of rounding mode changes, which is now allowed by the regrouping of operations on intervals.

Another data point illustrating the improvements is that we measured the number of Dag nodes allocated. For SC<Lazy exact nt<Gmpq>>, 29 million nodes were allocated, while for Lazy kernel<SC<Gmpq>> only 2.5 million nodes were needed. So we have won a factor of more than 10, due to the regrouping allowed by our design.

Note that the algorithm uses random data, hence it does not produce many filter failures, so almost no exact evaluation is performed. Another thing worth noticing is that it uses relatively simple 2D primitives. More complex primitives, especially in higher dimensions, should render the method even more advantageous. Finally, real-world geometric applications tend to produce more combinatorial output, hence the relative runtime cost of primitives is smaller, so the slow down factor is lower in those cases.

We now show experiments on the computation of the flow complex, a 3D graph built on top of the 3D Delaunay triangulation, which is used, for example, in a 3D surface reconstruction algorithm. The flow complex algorithm is described in more detail in [START_REF] Cazals | Robust construction of the threedimensional flow complex[END_REF]. The data sets are composed of 12593 and 54811 points, respectively. They represent scanned objects and are referred to as Mecanic and Fish. The implementation is based on the 3D Delaunay triangulation of Cgal. The latter exclusively uses geometric predicates, and no construction, whereas the flow complex computation itself uses cascaded constructions that demonstrate our point (repeated projections of points on the segments representing the Voronoi edges). This benchmark has been run in the same conditions as the previous one, except that the system had 2GB of memory, a 3GHz processor and was running Linux Fedora Core 5. Table 2 illustrates a gain of a factor of 4 in memory, and about 30% in time. Note that the reference kernel SC<double> lead to crashes of the algorithm on both data sets due to numerical issues. 2 Benchmarks comparing different kernels for the construction of the flow complex.

Open Design Questions

Here is a list of open questions related to our framework.

The first question concerns the regrouping of operations as described in Section 4.3. Our framework asks the user to pass functors specifying the level at which the regrouping of operations is made. In Cgal this is not a problem since the primary interface of the kernel towards the geometric algorithms is a list of functors. However it has the drawback of not being automatic. We can think of approaches based on expression templates [START_REF] Veldhuizen | Expression templates[END_REF] that would automatically detect sequences of operations and regroup them. Unfortunately, expression templates are limited to single statements and they tend to increase compilation times considerably. Could there be a way to extend the automatic regrouping to more than single statements? Maybe the auto keyword of C++0x will allow to propagate this through several statements? Or maybe the axiom feature of the concept extension of C++ [START_REF] Gregor | Concepts: linguistic support for generic programming in C++[END_REF] could be used to specify this kind of transformation.

Another question is whether similarly delayed computations are used in other areas, and if yes, whether it is possible to find out a common design, more general than the one we propose.

Conclusion and Future Work

In this paper we have presented a generic framework that implements lazy exact geometric computations, motivated by the needs for robustness and efficiency of geometric algorithms. This framework enables the delay of costly exact evaluations using multi-precision arithmetic when the faster interval arithmetic suffices.

The proposed design is easily extensible to new geometric primitives (predicates and constructions), as well as new geometric objects. It is based on a template family for representing lazy objects, as well as generic functor adaptors that produce them.

Our implementation is now part of Cgal public releases, and its entire geometry kernel already benefits from it.

Future work in this area will consist of various added special-case optimizations as well as generalizations. It is for example possible to refine the filtering scheme by growing the precision incrementally instead of switching directly to full multi-precision computation in case of insufficiency of precision of the intervals. We would also like to study the possibilities of merging the Filtered predicate and Lazy construct functor adaptors. Possible optimizations for specific cases also can be done, using faster schemes than interval arithmetic (so-called static filters). Moreover, the current way of providing a full kernel is by a list of types for the objects and functors, which is generated using of the preprocessor. We will therefore try to provide a better design on this particular point.

The authors wish to thank the anonymous reviewers for their careful work.

A Benchmark Code

#include <CGAL/Simple_cartesian.h> #include <CGAL/Lazy_kernel.h> #include <CGAL/Gmpq.h> #include <CGAL/Lazy_exact_nt.h> #include <CGAL/intersections.h> #include <CGAL/Timer.h> #include <CGAL/Memory_sizer.h> using namespace CGAL; // Choose a kernel: } std::cout << "orientation results : (-) = " << negative_ort << " (+) = " << positive_ort << " (0) = " << collinear_ort << std::endl; t.stop(); std::cout << "Total time = " << t.time() << std::endl; std::cout << "Total memory = " << ((Memory_sizer().virtual_size() -init_mem) >>10) << " KB" << std::endl; }

Fig. 1 .

 1 Fig. 1. Convex hull and Delaunay triangulation of a set of points in the plane.

  template < typename F, typename... Args > tuple<X...> apply(F f, tuple<Args...> args);Finally, the bind tuple function binds the first arguments of a functor f to values passed as a tuple. It returns a functor of type G, which simply forwards its arguments to f after those of the tuple.template < typename F, typename... Args > G bind_tuple(F f, tuple<Args...> args);

  mode = FE_UPWARD) : bak(std::fegetround()) { std::fesetround(mode); } ~Protect_FPU_rounding() { std::fesetround(bak); } }; Interval operator+(Interval a, Interval b) { Protect_FPU_rounding P; return Interval(-((-a.lower()) -b.lower()), a.upper() + b.upper()); }

  class Unsafe_comparison {}; bool operator<(Interval a, Interval b) { if (a.upper() < b.lower()) return true; if (a.lower() >= b.upper()) return false; throw Unsafe_comparison(); }

  operator()(K1::Point_2 p) const { return K2::Point_2( cv(p.x()), cv(p.y()) ); } . . . };

  template < typename O1, typename K1, typename K2 > struct Type_mapper { typedef . . . type; };

Fig. 3 .

 3 Fig. 3. Example Dag: √ x + √ yx + y + 2 √ xy.

Fig. 4 .

 4 Fig. 4. The Dag represents the midpoint m of an intersection point i and the orthogonal projection p of a point b onto a line l. Testing whether a, m, and b are collinear may trigger an exact construction.

Fig. 5 .

 5 Fig. 5. UML diagram of the classes involved in the nodes of the Dag.

•

  Lazy kernel<SC<Gmpq>> (2) is similar to the previous one, but it does not include the additional optimization that consists in eliminating rounding mode changes between consecutive interval computations, • finally, SC<double> is the Simple cartesian representation kernel instantiated with double. It is given for reference as it does not provide any guar-antee of robustness, since it does not use any interval nor multi-precision arithmetic. It only shows what the optimal performance could be.These benchmarks have been compiled using the Gnu g++ compiler version 4.1 with the -O2 optimization option. The code was executed on a Pentium-M laptop clocked at 1.7 GHz, equipped with 1 GB of RAM and 1 MB of cache, running the Fedora Core 3 Linux distribution. Timings are reported in seconds and memory usage in megabytes.

	Kernel	Time Memory
	SC<Gmpq>	70	70
	SC<Lazy exact nt<Gmpq>>	7.4	501
	Lazy kernel<SC<Gmpq>> (2)	3.6	64
	Lazy kernel<SC<Gmpq>>	2.8	64
	SC<double>	0.72	8.3
	Table 1
	Benchmarks comparing different kernels on an artificial algorithm.

This research is supported in part by the IST Programme of the EU as a Sharedcost RTD (FET Open) Project under Contract No IST-006413 (ACS -Algorithms for Complex Shapes).


Jean Sulem

Coupled thermo-hydro-mechanical processes in fault zones during rapid slip

Coupled thermo-hydro-mechanical processes in fault zones during rapid slip

INTRODUCTION

During the rupture of a fault, an earthquake occurs because the frictional resistance to slip on the fault walls decreases with increasing slip, causing an acceleration of sliding. To quantify the energy dissipated by an earthquake and assess the hazard of future ruptures, it is critical to understand the mechanics of slip weakening, i.e. how and how much fault friction drops in due course of the rupture.

The physical processes which occur during an earthquake exhibit many coupled phenomena as large variations of stress, pore pressure and temperature take place in the material. Thermo-poromechanical couplings due to shear heating [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]) can be associated to phase transition such as vaporization of the pore fluid at shallow depth [START_REF] Sulem | Thermo-Poro-Mechanical Properties of Clayey Gouge and Application to Rapid Fault Shearing[END_REF], melting of the rock at seismogenic depth [START_REF] Boullier | Fluid inclusions in pseudotachylites from the Nojima fault, Japan[END_REF][START_REF] Otsuki | Fluidization and melting of fault gouge during seismic slip: Identification in the Nojima fault zone and implications for focal earthquake mechanisms[END_REF], and to chemical effects such as dehydration of minerals [START_REF] Brantut | High-velocity frictional properties of a clay-bearing fault gouge and implications for earthquake mechanics[END_REF][START_REF] Brantut | Thermo-chemical pressurization of faults during coseismic slip[END_REF] or decomposition of carbonates (Sulem and Famin, 2009). Different competing effects may counter balance one another depending on the kinetics of the various physical processes.

MECHANISMS OF FAULT WEAKENING

Shear localization in faults

Field observations of mature faults, i.e. faults that have experienced a large slip, show a generally broad zone of damaged rock, but nevertheless suggest that shear in individual earthquakes occur in very narrow localized zones of few millimeters thick or even less. For example, during the drilling of Aigion Fault in the southern shore of the Gulf of Corinth, Greece, (Cornet et al, 2004a,b), the fault was intercepted at depth 760m, dipping at an angle of about 60°. The heart of the fault was found to be a 1m-thick band made of clay-rich material. Inside the clayey core, a clear shearing surface with marked slip lines is visi-ble on a plane that makes a 68° angle with respect to the core axis [START_REF] Sulem | Stress orientation evaluated from strain localization analysis in Aigion fault[END_REF]. This failure surface was not induced by the coring and extraction process but is indeed a slip plane as clear striation is observed at the interface (Fig. 1). Such an image clearly shows that the so-called failure plane is in reality a thin band of rapidly deforming material whereas the surrounding rock is moving, at least in the first stage of the process, as a rigid body. As the fault slip zones are very thin, thermal effects during rapid slip may be of primary importance. Fault zones exhibit commonly the presence of fluid interacting with the rock and thus coupled effects including shear heating, pore fluid pressurization, thermal mineral decomposition, should be taken into account in the understanding of fault weakening process.

The heat flow paradox

When an earthquake occurs, elastic strain energy is converted into radiated energy, i.e. the energy radiated by seismic waves, fracture energy, i.e. the energy associated with expanding the rupture area over the fault zone and frictional energy, i.e. energy required to overcome fault friction. It is generally estimated that more than 90% of the frictional work is converted into heat [START_REF] Cocco | On the mechanical work absorbed on faults during earthquake ruptures[END_REF].

Considering the extreme thinness of the slip zone (<1-5mm) as mentioned above and assuming that the shear strength along a fault may be represented by a Coulomb friction law with friction coefficients of rocks relatively high (~ 0.6-0.85), one should observe measurable heat outflow near major faults and extensive melting along exhumed faults. However these effects are generally not found. The actual strength of San Andrea fault has been debated for more than 35 years because heat flow measurements did not show a narrow heat flow anomaly centered on the fault that would be expected from a conductive model of generated heat [START_REF] Lachenbruch | Frictional heating, fluid pressure and the resistance to fault motion[END_REF]Sass, 1988, Scholz, 2006). Two main lines of explanation, as summarized by [START_REF] Rice | Thermo-and hydro-mechanical processes along faults during rapid slip[END_REF], have been explored to explain this paradox. The first one assumes that mature faults are weak either because the fault core material is different from the host rock and has a much lower friction coefficient due to the presence of clay or talc or because the pore pressure is high and nearly lithostatic. The other alternative is that major faults are statically strong but dynamically weaken during seismic slip. What physical process governs this dynamic weakening is still an open question.

Dynamic weakening mechanisms

As emphasized above, the slip zones are of extreme thinness and consequently it is expected that thermal weakening mechanisms are of major importance.

Thermal weakening of the friction coefficient

Flash heating and shear weakening of frictional micro-asperity contacts is a mechanism that has been first proposed for dry metal friction and recently applied to rock friction by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Beeler | Constitutive relationships and physical basis of fault strength due to flash heating[END_REF], Noda et al. 2009. From the microscopic point of view, the contact between two rough surfaces is effective on micro-asperities with an area which is much smaller than the nominal contact area. The stresses supported by these asperities are much higher than the macroscopic applied stresses so that the local rate of heat production and consequently the local temperature rise during sliding is large and maybe sufficient to activate thermal weakening of the contacts.

Other mechanisms are related to the lubrication of the contact zone with the formation of silica gels [START_REF] Toro | Friction falls toward zero in quartz rock as slip velocity approaches seismic rates[END_REF]). This mechanism is only possible at large slip in wet silica-rich fault zones. The formation of a molten layer along the sliding interface also lubricates the contact zone and reduces the friction resistance. Field observations reveal the presence of veins of pseudotachylytes (solidified friction-induced melts) (e.g. [START_REF] Boullier | Fluid inclusions in pseudotachylites from the Nojima fault, Japan[END_REF][START_REF] Otsuki | Fluidization and melting of fault gouge during seismic slip: Identification in the Nojima fault zone and implications for focal earthquake mechanisms[END_REF][START_REF] Toro | Earthquake rupture dynamics frozen in exhumed ancient faults[END_REF], 2006). However, onset of macroscopic melting can only occur deep in the fault zone at seismogeneic depth and pseudotachylytes are relatively rare. Whether this scarcity of melt products is due to the fact that conditions for the onset of macroscopic melting do not generally occur or to the fact that they are not preserved is still an open debate. As discussed below, other mechanisms such as thermal pore fluid pressurization and thermal decomposition of minerals can limit the temperature rise and prevent melting to occur.

Thermal pore fluid pressurization

Fluid pressurisation under shear heating is a mechanism of thermal softening. The principle of slip weakening by thermal pressurization is based on the fact that pore fluids trapped inside the fault zone are put under pressure by shear heating, thus inducing a reduction of the effective mean stress, and of the shearing resistance of the fault plane [START_REF] Lachenbruch | Frictional heating, fluid pressure and the resistance to fault motion[END_REF], Mase and Smith, 1985[START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF][START_REF] Sulem | Thermo-poro-mechanical properties of the Aigion fault clayey gouge -application to the analysis of shear heating and fluid pressurization[END_REF][START_REF] Sulem | Stress orientation evaluated from strain localization analysis in Aigion fault[END_REF][START_REF] Ghabezloo | Stress dependent thermal pressurization of a fluid-saturated rock[END_REF]. This mechanism has also been suggested for weakening in catastrophic landslides [START_REF] Habib | Sur un mode de glissement des massifs rocheux[END_REF], 1975[START_REF] Anderson | An earthquake induced heat mechanism to explain the loss of strength of large rock and earth slides[END_REF], Voigt and Faust, 1982[START_REF] Vardoulakis | Dynamic thermo-poro-mechanical analysis of catastrophic landslides[END_REF][START_REF] Veveakis | Thermoporomechanics of creeping landslides: The 1963 Vaiont slide, northern Italy[END_REF]. The presence of clay material in fault zones also affects thermal pressurization as possible collapse of the clay under thermal loading may activate fluid pressurization [START_REF] Sulem | Experimental characterization of the thermo-poromechanical properties of the Aegion fault gouge[END_REF][START_REF] Sulem | Thermo-poro-mechanical properties of the Aigion fault clayey gouge -application to the analysis of shear heating and fluid pressurization[END_REF][START_REF] Sulem | Thermo-Poro-Mechanical Properties of Clayey Gouge and Application to Rapid Fault Shearing[END_REF].

The values of the undrained thermal pressurization coefficient Λ , defined as the pore pressure increase due to a unit temperature increase in undrained condition, is largely dependent upon the nature of the material, the state of stress, the range of temperature change and the induced damage. In the literature we can find values that differ of two orders of magnitude: In Campanella and Mitchell (1968), different values are found from 0.01MPa/°C for clay to 0.05MPa/°C for sandstone. Palciauskas and Domenico (1982) estimate a value of 0.59MPa/°C for Kayenta sandstone. On the basis of Sultan (1997) experimental data on Boom clay, [START_REF] Vardoulakis | Dynamic thermo-poro-mechanical analysis of catastrophic landslides[END_REF] estimates this coefficient as 0.06MPa/°C. For a clayey fault gouge extracted at a depth of 760m in Aigion fault in the Gulf of Corinth (Greece), the value obtained by [START_REF] Sulem | Experimental characterization of the thermo-poromechanical properties of the Aegion fault gouge[END_REF] is 0.1MPa/°C and for intact rock at great depth, the value given by [START_REF] Lachenbruch | Frictional heating, fluid pressure and the resistance to fault motion[END_REF] is 1.5MPa/°C. For a mature fault at mid-seismogenic depth (7km) at normal stress of 196MPa, ambient pore pressure of 70MPa, and ambient temperature of 210°C, [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF] estimates this coefficient as 0.92MPa/°C in case of intact fault walls and 0.31MPa/°C in case of damage fault wall. On the basis of an experimental study on a sandstone, [START_REF] Ghabezloo | Stress dependent thermal pressurization of a fluid-saturated rock[END_REF] have emphasized the importance of pressure-dependent compressibility of rocks and pressure and temperature dependent compressibility of pore water on thermal pressurization and have found values between 0.02 and 0.7MPa/°C.

If for simplicity one assumes that the pressurization coefficient Λ does not change with temperature and pore pressure, an analytical solution can be obtained for the evolution of the temperature and of the pore pressure in an infinite long shear band under uniform shear for no heat and fluid transfer [START_REF] Sulem | Thermo-poro-mechanical properties of the Aigion fault clayey gouge -application to the analysis of shear heating and fluid pressurization[END_REF][START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]. This undrained adiabatic limit is applicable as soon as the slip event is sufficiently rapid considering the low permeability and thermal diffusivity of fault materials. This solution shows that the pore-pressure increases with the slip displacement from its initial hydrostatic value 0 f P towards its geostatic limit n σ (normal stress acting on the shear band), which corresponds to full fluidization. The corresponding temperature rise is simply obtained by dividing the pore pressure rise

0 n f P σ -
by the undrained pressurization coefficient Λ . For example, if one considers a fault at 7km depth with an initial hydrostatic pore pressure at 70MPa, a geostatic stress at 180MPa, an initial temperature at 210°C and an undrained pressurization coefficient comprised between 0.1 and 1MPa/°C, the temperature rise obtained is between 110 and 1100°C. This temperature rise must be seen as a lower bound because, as discussed by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF], the adiabatic and undrained solution appears to be unstable for localization so that the deformation actually localizes in very narrow zones with thickness between few millimeters and few tens of millimeters as observed in mature fault zones. Another idealized situation is met if one assumes that the strain localization zone is infinitely small. This situation corresponds to the model of slip on a plane for which an analytical so-lution has been proposed by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF]. The solution of slip on a plane gives an upper bound for the temperature field and the maximum temperature rise is obtained as ( ) ( , and taking for the undrained pressurization coefficient Λ an average value of 0.5MPa/°C, the temperature rise obtained is between 692 and 1440°C. These simple analytical solutions show that the temperature rise which occurs during a rapid slip event may easily reach the temperature limit for the onset of mineral decomposition. This will be discussed in the next section.

Thermal decomposition of minerals

Frictional heat released in a fault zone during seismic slip may induce chemical reactions. Thermal decomposition of minerals for large increase of temperature is potentially a major source of fluid release in the system which may enhance the generation of additional pore pressure excess.

Field evidences that such reactions may take place during an earthquake have been reported in several recent studies. The Taiwan Chelungpu-fault Drilling Project (TCDP) undertaken in 2002 to investigate the faulting mechanism of the 1999 Taiwan Chi-Chi earthquake allowed to recover core samples inside the active fault-zones within the Chelungpu fault. The analysis of these cores showed that the fault core was depleted in carbon relative to the damage zone, and the depletion was attributed to a decomposition of carbonates induced by frictional heat (Hirono et al., 2006a;[START_REF] Hirono | A chemical kinetic approach to estimate dynamic shear stress during the 1999 Taiwan Chi-Chi earthquake[END_REF][START_REF] Famin | Earthquake produce carbon dioxide in crustal faults[END_REF]. Hirono et al (2006b) have also measured a higher magnetic susceptibility and a lower density inside the fault core than outside. This high magnetic susceptibility which is degree to which a material can be magnetised by an external magnetic field can be attributed to the decomposition of siderite (FeCO 3 ) into magnetite (Fe 3 O 4 ) with degassing of CO 2 . Dehydration of hydrous clay minerals during during seismic slip have also been reported by [START_REF] Hirono | Clay mineral reactions caused byfrictional heating during an earthquake: An example from the Taiwan Chelungpu fault[END_REF]. These authors have observed that the fault core is relatively low in kaolinite and smectite contents as compared to the surrounding zones. This is attributed to the dehydration of interlayer (adsorbed) water (above 200°C) and the dehydroxylation (dehydration of structural water) of kaolinite (above 550°C) and smectite (illitization of smectite above 400°C).

Such reactions of decarbonation of calcite and dehydration of clay minerals have been reproduced experimentally during high velocity friction experiments [START_REF] Han | Ultralow friction of carbonate faults caused by thermal decomposition[END_REF][START_REF] Brantut | High-velocity frictional properties of a clay-bearing fault gouge and implications for earthquake mechanics[END_REF]).

CO 2 degassing during seismic slip

There is evidence of CO 2 release in several active crustal faults. In addition to the relatively low content in inorganic carbon (mainly calcite) in fault core as mentioned above, other field evidences have been reported. In the Corinth rift (Greece) for example, chemical analyzes of water springs near the seismogenic Heliki and Aegion faults revealed an anomalously high content of dissolved CO 2 compared with the regional values [START_REF] Pizzino | Fluid geochemistry along the Eliki and Aigion seismogenic segments Gulf of Corinth, Greece[END_REF]). The surface trace of the San Andreas fault also displays a positive anomaly of CO 2 fluxes [START_REF] Lewicki | CO 2 degassing along the San Andreas fault, Parkfield, California[END_REF], and this CO 2 comes from a shallow source, not from the mantle [START_REF] Lewicki | Shallow soil CO2 flow along the San Andreas and Calaveras Faults, California[END_REF]. Moreover, there is growing evidence that CO 2 release coincides with seismic slip in crustal faults, active and/or exhumed. [START_REF] Sato | Geochemical changes in anomalously discharged groundwater in Awaji Island -after the 1995 Kobe earthquake[END_REF] also reported that the HCO 3 -concentration increased by 30 wt% in springs located near the Nojima fault (Japan) immediately after the 1995 Kobe earthquake. This carbon discharge, together with other coseismic geochemical anomalies, decreased gradually to normal values in the following ten months. A microinfrared analysis of exhumed pseudotachylites (i.e. friction induced melts produced by seismic slip) from the Nojima fault revealed a carbon supersaturation in the melts, and the quantity of CO 2 released by friction melting during the 1995 Kobe earthquake was evaluated to 1.8 to 3.4 10 3 tons [START_REF] Famin | Earthquake produce carbon dioxide in crustal faults[END_REF]. These authors emphasize that the carbon had a biogenic isotopic signature both in springs and in the fault rocks, and did not correspond to carbon dioxide degassing from the mantle (Ueda et al., 1999;[START_REF] Arai | Behavior of gases in the Nojima Fault Zone revealed from the chemical composition and carbon isotope ratio of gases extracted from DPRI 1800 m drill core[END_REF][START_REF] Arai | Sealing process with calcite in the Nojima active fault zone revealed from isotope analysis of calcite[END_REF]Lin et al., 2003). The CO 2 of the Nojima fault might thus be attributed to the decomposition of biogenic carbonates. In the Central Apennines, [START_REF] Italiano | CO 2 degassing over seismic areas: The role of mechanochemical production at the study case of central Apennines[END_REF] also reported enhanced fluxes of crustal CO 2 (i.e. not mantellic) during the 1997-1998 seismic crisis of major faults, and proposed that coseismic decarbonation was responsible for the CO 2 emission.

This production of fluid which enhances the increase of pore fluid pressure is in competition with the change of porosity related to the reduction of solid volume induced by the mineral decomposition. Obviously the competing effects of temperature rise (due to shear heating), thermal fluid pressurization, mineral decomposition and porosity/permeability increase induce strong non-linearities in the problem of pore pressure and temperature evolution of a rapidly sheared fault zone. The kinetics of the chemical reaction of mineral decomposition is also a major factor that competes with heat and fluid diffusion processes. The endothermic character of the chemical reaction plays also a major role in the energy balance of the system as a large part of the frictional heat is consumed into the mineral decomposition [START_REF] Hamada | Energy taken up by co-seismic chemical reactions during a large earthquake: An example from the 1999 Taiwan Chi-Chi earthquake[END_REF]). As shown by Sulem and Famin (2009) and [START_REF] Brantut | Thermo-chemical pressurization of faults during coseismic slip[END_REF], these chemi-cal reactions which are thermally activated limit the temperature increase which tends to stabilize to a value slightly higher than the temperature of onset of the reaction.

THERMO-HYDRO-MECHANICAL MODEL-ING OF RAPID SLIP IN CARBONATE ROCKS

Kinetics of thermal decomposition of calcite

We focus here on the kinetics of chemical decomposition of calcite (decarbonation) because carbonates are present in every fault zones from the ductilebrittle transition (~15 km) to the subsurface, and because, as mentioned above, there is evidence of CO 2 release in several active crustal faults.

The rate of emitted CO 2 is evaluated using the kinetics of the chemical reaction of calcite thermal decomposition:

3 2 CaCO CaO+CO → (1)
The reacted fraction, α (0 ≤ α ≤ 1; α = 1 if the entire substance is reacted) is expressed by the Arrhenius equation (L'vov 2002):

( ) exp a E f A t RT α α ∂   = -   ∂   (2)
where ( ) f α is a kinetic function determined by the reaction mechanism, A is a constant (pre-exponential term of the Arrhenius law), E a is the activation energy of the reaction, R is the gas constant (8.31447JK-1mol-1). In the following we will take the values corresponding to CaCO3 mixed with silica:

-1 15 -1 319000Jmol , 2.95 10 s more et al., 1996). For thermal decomposition of calcite the kinetic function given by [START_REF] Criado | The effect of the CO2 pressure on the thermal decomposition kinetics of calcium carbonate[END_REF] is ( )

a E A = = × (Dolli-
2/3 ( ) 1 f α α = - .

Position of the problem

Considering that the length scales in the direction parallel to the fault over which the thermo-poromechanical fields vary are much larger than in the direction normal to it, we analyse here a 1D problem. We consider a rapidly deforming and infinitely long shear band of thickness h consisting of fluidsaturated material. This shear band begins to undergo slip δ at a time 0 t = with an imposed overall slip-rate V in the x-direction, as shown in figure 1. Inside such a shear-band the pore pressure p P , the temperature T , and the velocity v are assumed to be functions only of time t and of the position z in the direction normal to the band (Fig. 2).

For the pore pressure and temperature range considered here for which the decarbonatation of CO 2 is active (T ≥ 700°C and P p ≥ 70MPa), water and carbon dioxide are in supercritical state. They are miscible and their compressibility, viscosity, density and thermal expansion coefficient are nearly equal. Thus, for simplicity, it is assumed that the porous space is saturated with a unique fluid, with the thermodynamic properties of water.

σ σ σ σ n τ τ τ τ P p0 , T 0 P p ,T h σ σ σ σ n τ τ τ τ q h q f v slip σ σ σ σ n τ τ τ τ σ σ σ σ n τ τ τ τ x z u(z,t) P p (z,t),T(z,t) q h (z,t) q f (z,t) σ σ σ σ n τ τ τ τ P p0 , T 0 P p ,T h σ σ σ σ n τ τ τ τ q h q f v slip σ σ σ σ n τ τ τ τ σ σ σ σ n τ τ τ τ x z u(z,t)
P p (z,t),T(z,t) q h (z,t) q f (z,t)

Governing equations

The governing equations are derived from the conservation and transport equations [START_REF] Vardoulakis | Bifurcation analysis in geomechanics[END_REF]. The reader can refer to Sulem and Famin (2009) and Sulem et al. (2009) for the details of the derivations.

Fluid mass conservation

The production/diffusion equation for the pore pressure is obtained from the mass balance equations and is expressed as the sum of four terms: the diffusion term, the thermal pressurization term, the term corresponding to the fluid released in the system by the chemical reaction and the term corresponding to the effect on pore pressure of inelastic porosity change:

1 1 1 p f d d f f P q m n T t z t t t ρ β ρ β β ∂ ∂ ∂ ∂ ∂ = - + Λ + - ∂ ∂ ∂ ∂ ∂ (3)
In the above equation, f ρ is the density of the saturating fluid;

β is the (elastic) volumetric pore fluid storage coefficient:

( )

f n n β β β = +
where n is the pore volume fraction (Lagrangian porosity), 1

f f f p T P ρ β ρ   ∂ =     ∂  
is the pore fluid compressibility and 1

n p T n n P β   ∂ =     ∂  
is the pore volume compressibility; f q is the fluid mass flux and is evaluated assuming Darcy's law for fluid flow with viscosity f η through a material with permeability f k :

f p f f f P q k z ρ η ∂ = - ∂ ; - f n n f λ λ β β - Λ = +
is the undrained thermo-elastic pressurization coefficient [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF], where 1

p f f f P T ρ λ ρ ∂   = -   ∂  
is the pore fluid thermal expansion coefficient and 1

p n P n n T λ ∂   =   ∂  
is the thermal expansion coefficient of the pore volume, which is equal to the thermal expansion coefficient of the solid fraction for thermo-poro-elastic materials;

-m d is the mass per unit volume of the fluid released in the system by mineral decomposition; Using the kinetic function for the considered chemical reaction (equation ( 2)) the rate of the released fluid is simply given by: ( )

0 d d m m m t t α ∂ ∂ = - ∂ ∂ (4)
where

3 0 CaCO 0 (1 ) m n χρ = -
is the total mass of calcite that can be decomposed (n 0 is the initial porosity, and χ is the ratio between the molar mass M CO2 of CO 2 (44g/mol) and the molar mass M CaCO3 of calcite (100g/mol)).

-

d n t ∂ ∂
is the rate of inelastic porosity change; this term includes the effect of co-seismic dilatancy, induced by fault surface irregularities and the effect of the decomposition of the solid phase. If we neglect the co-seismic dilatancy and assume that inelastic porosity change is only due to the change of solid mass, we can relate it to the mass of released fluid:

3 3 3 3 2 2 CaCO CaO CaCO CaCO CaCO CaO CaCO CO CaCO CO 1 1 1 1 with = d d m n m m t t t t M M M M ξ ρ ρ ξ ρ ρ ∂ ∂ ∂ ∂ = - - = ∂ ∂ ∂ ∂ - (5)
where ρ CaCO3 (resp. ρ CaO ) and M CaCO3 (resp. M CaO )

are the density and the molar mass of CaCO 3 (resp. CaO).

Energy conservation

Let E F be the rate of frictional heat produced during slip and E C the rate of heat used in the endothermic thermal decomposition. The equation of conservation of energy is expressed as

h F C q T C E E t z ρ ∂ ∂ = -- ∂ ∂ (6)
where C ρ is the specific heat per unit volume of the fault material in its reference state and q h is the heat flux. It is assumed that all the plastic work is con-verted into heat and thus

∂ = ∂ F v E z
τ . On the other hand, the rate of heat used in the chemical reaction can be expressed from the enthalpy change of the reaction and the rate of reacted fraction as

0 d C r T m E H t ∂ ∂ = ∆
. The heat flux is related to the temperature gradient by Fourier's law,

h T T q k z ∂ = -∂
where k T is the thermal conductivity of the saturated material. In equation ( 6) it is assumed that all heat flux is due to heat conduction neglecting heat convection by the moving hot fluid. This assumption is justified by the low pore volume fraction and the low permeability of fault gouges [START_REF] Lee | Frictional heating and pore pressure rise due to a fault slip[END_REF]. However, this assumption should be relaxed in case of a seismically-induced high permeability in the wall rocks of the fault.

Momentum balance

The 1D-momentum balance equation reads as

v z t τ ρ ∂ ∂ = ∂ ∂ (7)
As discussed by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF], the effect of even large accelerations like several times the acceleration of the gravity g is insignificant over the small length scales in the z-direction normal to the fault where the heat and fluid diffusion process are taking place during rapid slip and very high values of pressure and temperature gradients. For example, assuming an acceleration of 10g and a specific mass of the material of 2500kg/m 3 would result in a change of 2.5×10 -2 MPa/m for τ. The relevant length scale in zdirection is only few cm and thus the variation for τ can be neglected and mechanical equilibrium can be assumed:

0 z τ ∂ = ∂ (8)
Consequently, as the shear stress is constant in space, the Coulomb friction law cannot be assumed to be met in all deforming regions unless the pore pressure is also constant in space as it is the case in the undrained adiabatic limit. It is thus assumed that the frictional resistance is proportional to the mean effective stress inside the band:

/2 /2 1 ( ) ( , ) h n p h t f P t d h τ σ ξ ξ -   = -     ∫ ( 9
)
where f is the friction coefficient of the gouge and h is the width of the shear-band. From the equation of momentum balance we obtain a linear distribution for the velocity field corresponding to homogeneous straining:

/ 2 z v V h = (10) 4 NUMERICAL EXAMPLES

Thermal pressurization in faults at shallow depth

The first example presented here refers to Aigion fault (Cornet et al, 2004a,b). The experimental characterization program performed on fault zone cores extracted from Aigion fault at a depth of 760m has shown that the Aigion fault gouge is a clay-rock mixture [START_REF] Sulem | Experimental characterization of the thermo-poromechanical properties of the Aegion fault gouge[END_REF]. The mineral composition of the clayey gouge exhibits a large content of quartz (74%), some hydrous minerals (chlorite, 2-3% and illite 8-12%) and also a small content of albite (8-13%) and hematite (2-3%). Although the clay fraction is relatively small, it has a significant influence on the thermo-mechanical properties of the material. The clayey core of the fault has a very low fluid permeability of about 10 -19 m 2 and exhibits contractant volumetric behavior when heated. This behavior is typical of normally consolidated clays [START_REF] Hueckel | Thermoplastic modelling of undrained failure of saturated clay due to heating[END_REF]Pellegrini, 1991, Sultan et al. 2002.) The sensitive parameter for the description of the thermo-poro-mechanical coupling is the thermal pressurization coefficient of the material. It was found that at a nominal effective stress corresponding to 760m deep conditions, this coefficient is about 0.1MPa/°C [START_REF] Sulem | Thermo-Poro-Mechanical Properties of Clayey Gouge and Application to Rapid Fault Shearing[END_REF].

Assuming that the density of the overburden is 2500 kg/m 3 , we estimate a total vertical stress of about 19MPa for a depth of 760m. The hydrostatic initial pore pressure is 7.6MPa and the temperature recorded from thermal logs at 760m is 28°C.

In order to demonstrate the importance of the shear band thickness, two computations assuming a thickness of 1mm and 0.5mm respectively have been performed. A constant slip velocity of 1m/s is imposed. The friction coefficient of the gouge material is taken equal to f=0.5 [START_REF] Sulem | Experimental characterization of the thermo-poromechanical properties of the Aegion fault gouge[END_REF]). In the first one, the liquefaction limit ( p n P σ = , i.e. full fluidization of the gouge), is reached after 0.13s of shearing (Fig. 3a). The corresponding temperature reached in the centre of the sheared is 142°C. In the second case for which the shear band thickness which is twice smaller (h=0.5mm), the computed results show that the temperature increase is much higher and that vaporization of the pore fluid occurs (Fig. 3b) before liquefaction. This state is reached at t=1.4s for which in the centre of the band P p =17MPa and T=356°C. Note that fluid vaporization is possi-ble only if it is assumed that the soil is not fully saturated.

We observe that at shallow depth, the pore fluid pressurization limits the temperature increase so that no mineral decomposition occurs.

Fault in carbonate rocks at mid-seismogenic depth

We consider here a fault at 7km depth in a carbonate rock with an initial hydrostatic pore pressure at 70MPa, a geostatic stress at 180MPa, an initial temperature at 210°C. The thermal expansion coefficient is

5 1 2.4 10 ( C) n λ - - = × °.
For the pore volume compressibility, we account for the effect of effective stress change due to pore pressure rise. The expression of the pore volume compressibility is obtained from poro-elasticity theory and is given by [START_REF] Rice | Heating and weakening of faults during earthquake slip[END_REF] as:

( ) ( ) 1 1 n d s n n β β β = -+ (11)
where β s is the compressibility of the solid phase

( 5 -1

1.25 10 MPa -× for calcite) and β d is the drained compressibility of the porous rock. We consider here the empirical expression for the effective stress dependent drained compressibility as proposed by [START_REF] Wibberley | Earthquake slip weakening and asperities explained by thermal pressurization[END_REF] for the Median Tectonic line fault zone in Japan:

( ) ( ) ( ) 4
The initial porosity of the rock is taken equal to 0.03. The density of the porous rock is taken equal to 3 2.6g/cm ρ =

, the density of the solid phase (calcite)

is 3 2.71g/cm s ρ =
, and the specific heat is 2.7MPa/(°C) C ρ = . We assume a slip velocity of 1m/s. In the computations we account for the changes of the physical properties of the saturating fluid (here water) (Wagner and Pruss, 2002).

The decomposition of carbonate can induce substantial change in the porosity of the rock which affects the permeability. It is known that there is no unique relationship between porosity and permeability applicable to all porous media and that the geological evolution process of the pore space influence the permeability-porosity relationship. The empirical power law f k n α ∝ is commonly used for geomaterials. The exponent α characterizes the porosity sensitivity of permeability and can take values ranging from 1 to 25 [START_REF] David | Laboratory measurement of compaction-induced permeability change in porous rocks: Implications for the generation and maintenance of pore pressure excess in the crust[END_REF]. The high values of α correspond in general to rocks with a high porosity whereas for low porosity rocks law values of α are obtained. The commonly used cubic Carman Kozeny permeability law is assumed here to take into account the effect of porosity change due to mineral decomposition on the permeability of the rock:

3 2 0 0 0 1 1 f f n n k k n n   -   =     -     (13)
where k f0 is the reference permeability corresponding to the reference porosity n 0 . The assumed permeability law has a strong effect on the numerical results and there is a need for experimental data on permeability changes resulting from the particular process of carbonate decomposition. According to [START_REF] Milsch | Reactioninduced fluid flow in synthetic quartz-bearing marbles[END_REF], the fluid transport occurs as a local process on the grain scale and therefore it is expected that inside the slip zone, where intense grain comminution and production of very fine particles occur, the reaction induces strong changes of permeability. The recent high velocity shear experiments on Carrara marble of [START_REF] Han | Ultralow friction of carbonate faults caused by thermal decomposition[END_REF] have shown that the friction coefficient decreases rapidly to values as low as 0.06 due to the thermal decomposition of calcite induced by frictional heating. Here, we take 0.1 f = . We assume that the initial permeability of the intact medium is The computed results are presented in figure 4. The evolution in time of the temperature and the pore pressure in the centre of the band is plotted in figure 4a. Considering the constant slip velocity of 1 m/s, this graph can also be seen as the evolution of temperature and pore pressure with accumulated slip. The corresponding shear stress is plotted in figure 4b. These results show the coupling effect of heat which induces first a pore pressure increase. When the decomposition of the carbonate rock begins at about 700 °C, the temperature increase is drastically slowed due to the energy consumed in the endothermic chemical reaction.

Two competing effects act on the evolution of the pore pressure: on one hand the production of CO 2 induces an additional fluid mass and thus a pressurization of the pore fluid; on the other hand the in-crease of porosity due to the solid decomposition induces an increase of the permeability of the medium which limits the pressurization. The pore pressure in the centre of the band exhibits a maximum of about 145MPa which does not exceed the total normal stress acting (180MPa) on the band which means that full liquefaction is not reached. As mentioned above, through the friction law the shear stress is related to the mean effective stress inside the band (figure 4b). Consequently, the shear stress decreases rapidly during initial pressurization and then increases again. Thus, the mineral decomposition of the rock can be seen as a mechanism of fault weakening in a first stage then fault re-strengthening in a second stage. The accumulated mass per unit area of the fault plane of emitted CO 2 is plotted in figure 5a. After one second, the temperature reaches a quasi constant state and the production rate is almost constant. These results show that the mass of emitted CO 2 after 2s is about 0.018g/cm2. The corresponding rate of emitted CO 2 is about 180t/s per km 2 of fault zone. The corresponding solid mass change is shown through the evolution of the porosity in the centre of the shear band (Fig. 5b). The temperature field around the centre of the band is shown in figure 6 for 2s t = , with the corresponding value of the porosity. These plots show that the temperature increase is localized in the central zone of the band. The porosity is affected only in zones where the temperature exceeds 700°C which for the considered computation corresponds to a width of about 4mm.

The effect of mineral decomposition is also shown in figure 7 where the results of the above computation are compared to the ones obtained without CO 2 emission. In this latter case only thermal pressurization occurs. If no mineral decomposition occurs, the pore pressure increase is slower and for the considered parameters, the maximum pore pressure is about 120MPa and consequently the temperature increase is much stronger and melting of the rock can possibly occur (the melting temperature is about 1600°C for calcite).

CONCLUSION

In the last thirty years, the study of earthquake mechanisms has emphasized the major role of rock friction. Earthquakes appear to be the result of a frictional instability and occur by sudden slippage along a pre-existing fault. Therefore, the understanding of the nature of the friction on faults during the interseismic periods and during earthquakes is essential to assess the hazards of seismic rupture.

The presence of fluids greatly affects the frictional resistance of a fault. Pore fluid pressurization under shear heating is a mechanism of thermal softening. When thermal softening prevails against strain or strain rate hardening, dynamic plastic localization in the form of adiabatic shear banding is known to occur at high strain rates under quasi adia-batic conditions. The phenomenon of adiabatic shear banding is well known in metal plasticity for metal cutting and impact problems (e.g. [START_REF] Wright | The physics and mathematics of adiabatic shear banding[END_REF]. The extreme thinness of the sheared zone revealed by field observations of mature faults suggest that it is the result of a plastic instability in the form of shear localization as shown by [START_REF] Rice | Large scale localization: Thermal weakening of faults during seismic slip[END_REF] because uniform undrained adiabatic shearing of a layer is actually unstable. Frictional weakening mechanisms allow to explain the acceleration of sliding during nucleation of an earthquake and also the lack of pronounced heat outflow, and lack of at least shallow frictional melting, along major tectonic faults.

Shear heating can lead to fluid vaporization at shallow depth and mineral decomposition at midseismogenic depth. Indeed, thermal decomposition of minerals appears to be an important physical process in the phenomenon of thermal heating and pore fluid pressurization during seismic slip. The combined effects of frictional heating, temperature rise, endothermic mineral decomposition, pore pres-sure rise, porosity and permeability increase result in highly coupled and competing processes.

In this paper, the effect of carbonate decomposition is emphasized because carbonates are present in every fault zones from the ductile-brittle transition (~15 km) to the subsurface, and because, as revealed by field observations, there is evidence of CO 2 release in several active crustal faults. A first essential result of the analysis is that the endothermic reaction of calcite decomposition limits the co-seismic temperature increase to less than ~800 °C (corresponding to the initiation of the chemical reaction) within a carbonate shear band under rapid slip. Decarbonation is only one of the possible thermal reactions of mineral decomposition. Phyllosilicates are common secondary minerals in faults and their reaction of thermal dehydration is also endothermic. Therefore, it is likely that the reaction-induced temperature limitation obtained in our model may apply qualitatively to most fault zones. In this case, mature faults with a long history of slip and a large cumulated displacement are likely more prone to reaction-induced temperature limitation than recent faults, because of their larger content in volatile-rich secondary minerals. In the recent study of [START_REF] Brantut | Thermo-chemical pressurization of faults during coseismic slip[END_REF], the authors show that the temperature rise in the slipping zone is stopped while the dehydration reaction progresses. This would provide another explanation to the notorious absence of positive heat flow anomaly on active crustal faults such as San Andreas (Lachenbruch and Sass 1980): a large part of the heat produced by friction would be consumed by endothermic reactions. Another consequence is that friction melting is hampered by endothermic calcite decomposition in carbonate fault rocks, and probably in other faults containing a sufficiently high proportion of volatile-rich secondary minerals. This is consistent with the relative scarcity of pseudotachylytes in mature faults.

The second essential implication of the above analysis is that decarbonation is a source of CO2 that significantly increases the slip-weakening effect of thermal pressurization. The pore fluid pressure exhibits an initial phase of increase due to thermal pressurization, then a sudden acceleration of generated pore pressure when the solid decomposition is activated. However, the increase of permeability limits the pore pressure that reaches a maximum and then decreases. The numerical results reproduce this pore pressure pulse and the initial fault weakening followed by a re-strengthening of the shear stress.

Our model of thermal pressurization, taking into account these coupled processes, provides a robust framework for estimating the dynamic friction of faults and the energy balance of earthquakes.

Figure 1 .

 1 Figure 1. (a) The clayey core retrieved from Aigion fault (Greece); (b and c) Striation of the slip plane.

Figure 2 .

 2 Figure 2. Model of a deforming shear-band with heat and fluid fluxes

  Figure 3. Aigion fault; Temperature and pore pressure vs time at the center of the band; slip velocity=1m/s (a) Shear band thickness=1mm, liquefaction occurs at t=0.13s ; (b) Shear band thickness=0.5mm, vaporization of the pore fluid is reached at 1.4s.

Figure 4 .

 4 Figure 4. Fault at 7km depth: (a) Evolution of temperature and pore pressure in the centre of the shear band, (b) evolution of the shear stress.

Figure 5 .

 5 Figure 5. Fault at 7km depth: (a) Evolution of the mass of emitted CO2, (b) evolution of the porosity in the centre of the band.

Figure 6 .

 6 Figure 6. Fault at 7km depth: temperature and porosity fields after 2s slip.

Figure 7 .

 7 Figure 7. Fault at 7km depth: Effect of CO2 emission on the evolution of temperature (a) and pore pressure (b) in the centre of the shear band.

ACKNOWLEDGEMENTS

This work is dedicated to the memory of Ioannis Vardoulakis as recognition for a long standing cooperation and a deep friendship. I also wish to thank Nicolas Brantut, Vincent Famin, Jim Rice and Alexandre Schubnel for fruitful discussions.