New aspects of tectonic stress inversion with reference to the TENSOR program

Abstract Analysis of tectonic stress from the inversion of fault kinematic and earthquake focal mechanism data is routinely done using a wide variety of direct inversion, iterative and grid search methods. This paper discusses important aspects and new developments of the stress inversion methodology as the critical evaluation and interpretation of the results. The problems of data selection and separation into subsets, choice of optimization function, and the use of non-fault structural elements in stress inversion (tension, shear and compression fractures) are examined. The classical Right Dihedron method is developed in order to estimate the stress ratio R, widen its applicability to compression and tension fractures, and provide a compatibility test for data selection and separation. A new Rotational Optimization procedure for interactive kinematic data separation of fault-slip and focal mechanism data and progressive stress tensor optimization is presented. The quality assessment procedure defined for the World Stress Map project is extended in order to take into account the diversity of orientations of structural data used in the inversion. The range of stress regimes is expressed by a stress regime index R’, useful for regional comparisons and mapping. All these aspects have been implemented in a computer program TENSOR, which is introduced briefly. The procedures for determination of stress tensor using these new aspects are described using natural sets of fault-slip and focal mechanism data from the Baikal Rift Zone.

Analysis of fault kinematic and earthquake focal mechanism data for the reconstruction of past and present tectonic stresses are now routinely done in neotectonic and seismotectonic investigations. Geological stress data for the Quaternary period are increasingly incorporated in the World Stress Map (WSM) (Mtfller & Sperner 2000;MUller et al. 2000;Sperner et al. 2003).
Standard procedures for brittle fault-slip data analysis and stress tensor determination are now well established (Angelier 1994;Dunne & Hancock 1994). They commonly use fault-slip data to infer the orientations and relative magnitude of the principal stresses.
A wide variety of methods and computer programs exist for stress tensor reconstruction. They are either direct inversion methods using least square minimization (Carey-Gailhardis & Mercier 1987;Angelier 1991;Sperner et al. 1993) or iterative algorithms that test a wide range of possible tensors (Etchecopar et al. 1981) or grid search methods (Gephart 1990b;Hardcastle & Hills 1991;Unruh et al. 1996). The direct inversion methods are faster but necessitate more complex mathemat-ical developments and do not allow the use of complex minimization functions. The iterative methods are more robust, use simple algorithms and are also more computer time intensive, but the increasing computer power reduces this inconvenience.
This paper presents a discussion on the methodology of stress inversion with, in particular, the use of different types of brittle fractures in addition to the commonly used fault-slip data, the problem of data selection and the optimization functions. Two methodologies for stress inversion are presented: new developments of the classical Right Dihedron method and the new iterative Rotational Optimization method. Both methods use of the full range of brittle data available and have been adapted for the inversion of earthquake focal mechanisms. The interpretation of the results is also discussed for two important aspects: the quality assessment in view of the World Stress Map standards and the expression of the stress regime numerically as a Stress Regime Index for regional comparisons and mapping. be obtained by contacting the first author. A program guideline is also provided with the program package.

Stress inversion methodologies
Stress analysis considers a certain volume of rocks, large enough to sample a sufficiently large data set of slips along a variety of different shear surfaces. The size of the volume sampled should be much larger than the dimensions of the individual brittle structures. For geological indicators, relatively small volumes or rock (100-1000 m 3) are necessary to sample enough fault-slip data, while for earthquake focal mechanisms, volumes in the order of 1000-10 000 km 3 are needed.
Stress inversion procedures rely on Bott' s (1959) assumption that slip on a plane occurs in the direction of the maximum resolved shear stress. Inversely, the stress state that produced the brittle microstructures can be partly reconstructed knowing the direction and sense of slip on variably oriented fault planes. The slip direction on the fault plane is inferred from frictional grooves or slickenlines. The data used for the inversion are the strike and dip of the fault plane, the orientation of the slip line and the shear sense on the fault plane. They are collectively referred to as fault-slip data. Focal mechanisms of earthquakes are also used in stress inversion. The inversion of fault-slip data gives the four parameters of the reduced stress tensor: the principal stress axes o-1 (maximum compression), o-2 (intermediate compression) and o-3 (minimum compression) and the Stress Ratio R = (o-2o-3)/(o-1 -o-3). The two additional parameters of the full stress tensor are the ratio of extreme principal stress magnitudes (o'3/o-1) and the lithostatic load, but these two cannot be determined from fault data only. We refer to Angelier (1989Angelier ( , 1991Angelier ( , 1994 for a detailed description of the principles and procedures of fault-slip analysis and palaeostress reconstruction. We are aware of the inherent limitations of any stress inversion procedures that apply also to the discussion proposed in this paper (Dupin et al. 1993;Pollard et al. 1993;Nieto-Samaniego & Alaniz-Alvarez 1996;Maerten 2000;Roberts & Ganas 2000).
The question was raised as to whether fault-slip inversion solutions constrain the principal stresses or the principal strain rates (Gephart 1990a). We will not discuss this question here, and leave readers to form their own opinions on how to interpret the inversion results. The brittle microstructures (faults and fractures) are used in palaeostress reconstructions as kinematic indicators. The stress inversion scepticals (e.g. Twiss & Unruh 1998) argue that kinematic indicators are strain markers and consequently they cannot give access to stress. Without entering in such debate, we consider here that the stress tensor obtained by the inversion of kinematic indicators is a function that models the distribution of slip on every fault plane. For this, there is one ideal stress tensor, but this one is only certainly active during fault initiation. After faults have been initiated, a large variety of stress tensors can induce fault-slip by reactivation.

Stress and strain relations
In fault-slip analysis and palaeostress inversion, we consider generally the activation of pre-existing weakness planes as faults. Weakness planes can be inherited from a sedimentary fabric such as bedding planes, or from a previous tectonic event. A weakness plane can be produced also during the same tectonic event, just before accumulating slip on it, as when a fault is neoformed in a previously intact rock mass. The activated weakness plane F can be described by a unit vector n normal to F (bold is used to indicate vectors). The stress vector o-acting on the weakness plane F has two components: the normal stress ~, in the direction of n and the shear stress -r, parallel to F. These two stress components are perpendicular to each other and related by the vectorial relation tr----v + ~'.
The stress vector o-represents the state of stress in the rock and has o-l, o-2 and o-3 as principal stress axes, defining a stress ellipsoid. The normal stress v induces a component of shortening or opening on the weakness plane in function of this sign. The slip direction d on a plane is generally assumed to be parallel to the shear stress component ~-of the stress vector o-acting on the plane. It is possible to demonstrate that the direction of slip d on F depends on the orientations of the three principal stress axes, the stress ratio R -= (o-2o-3)/(o-1 -o-3) and the orientation of the weakness plane n (Angelier 1989(Angelier , 1994. The ability of a plane to be (re)activated depends on the relation between the normal stress and shear stress components on the plane, expressed by the friction coefficient: = a tan lv ~ If the characteristic friction angle 4~ of the weakness plane F with the stress vector o-acting on it overcomes the line of initial friction, the weakness plane will be activated as a fault. Otherwise, no movement will occur on it. This line is defined by the cohesion factor and the initial friction angle 4)o.

Data types and their meaning in stress inversion
Determination of palaeostress can be done using two basic types of brittle structures: (1) faults with slip lines (slickensides) and (2) other fracture planes (Angelier 1994;Dunne & Hancock 1994).
In the following discussion as in program TEN-SOR, under the category 'fault', we understand only fault planes with measurable slip lines (slickensides), while we refer to all other types of brittle planes that did not show explicitly traces of slip on them as 'fractures'. Brittle structures other than the commonly used slickensides can also be used as stress indicators (Dunne & Hancock 1994). They are generally known as 'joints' or 'fractures'. The term 'joint' formally refers to planes with no detectable movements on them, both parallel and perpendicular to the plane (Hancock 1985). Others use the term joint to refer to large tensional fractures and in this case, tension fractures are mechanically similar to joints (Pollard & Aydin 1988). To avoid confusion, we use here the general term 'fracture' for all planar surfaces of mechanical origin, which do not bear slip lines on the fracture surface. Brittle fractures as understood here bear some information on the stress stage from which they are derived (Dunne & Hancock 1994). In order to assess the relationship between fracture planes and stress axes, it is important to determine their genetic classes. They may form as tension fractures (tension joints)  For faults, one can measure directly the orientations of the fault plane and the slip line, and determine the slip sense (normal, inverse, dextral or sinistral) using the morphology of the fault surface and secondary structures associated as in Petit (1987). In vectorial notation, this corresponds to n (unit vector normal to the fault plane) and d (slip direction, defined by the orientation of the slip line and the sense of slip on the plane). These form the so-called fault-slip data.
Data on slip surface and slip direction can also be obtained indirectly by combining several brittle structures, as suggested by Ragan (1973). Conjugate sets of shear fractures nl and n2 can be used for reconstructing the potential slip directions on the fracture planes, and to infer the orientations of principal stress axes. In conjugate fracture systems, 0-1 bisects the acute angle nlAn2, 0-2 is determined by the intersection between nl and n2, and 0-3 bisects the obtuse angle. The slip directions da and d2, respectively on shear planes defined by na and he, are perpendicular to o-2 in a way that the two acute wedges tend to converge. Similarly, the slip direction on a shear plane without observable slip line can be inferred if tension fractures are associated at an acute angle to the shear plane (Ragan 1973). For shear plane ns and associated tension fracture nt, 0-1 is parallel to the tension fracture, 0-2 is determined by the intersection between ns and nt, and 0-3 is perpendicular to the tension fracture (parallel to nt). The slip direction ds on the shear plane defined by ns is perpendicular to o-2 so the block defined by the acute angle nsAnt tend to move towards the block defined by the obtuse angle.
In these two cases the slip direction d on the shear planes n is reconstructed and these can be used in the inversion as additional fault-slip data.
For the stress inversion purpose, we distinguish between three types of brittle fractures.
• Tension fractures (plume joints without fringe zone, tension gashes, mineralized veins, magmatic dykes), which tend to develop perpendicular o-3 and parallel to 0-1. The unit normal vector nt represents an input of the direction 0-1. • Shear fractures (conjugate sets of shear fractures, slip planes displacing a marker), which form when the shear stress on the plane overcomes the fault friction. It corresponds to the input of a fault plane n~, but without the slip direction d. • Compression fractures (cleavage planes), which tend to develop perpendicular to 0-1 and parallel to 0-3. The vector nc represents an input of the direction 0-1.
For the stylolites, it is more accurate to use the orientation of the stylolite columns as a kinematic indicator for the direction of maximum compression (0-1) instead of the plane tangent to the stylolite seam (i.e. input of direction 0-1). Earthquake focal mechanisms are determined geometrically by the orientations of the p-and tkinematic axes bisecting the angles between the fault plane and the auxiliary plane. They can be determined also by the orientation of one of the two nodal plane (nl or n2) and the associated slip vector (da or de), or by the orientation of the two nodal planes (ni and n2) and the determination of the regions of compression and tension (i.e. input of either nl and n2 or p and t).
In addition to the orientation data, it is also important to record qualitative information for all fault-slip data: the accuracy of slip sense or fracture type determination (slip sense confidence level), a weighting factor (between 1 and 9, as a function of the surface of the exposed plane), whether the fault is neoformed or has been reactivated, the type and intensity of slip striae, the morphology and composition of the fault or fracture surface, estimations of the relative timing of faulting for individual faults, based on cross-cutting relationships or fault type.

Data selection and separation into subsets
After measuring the fault and fracture data in the field or compiling a catalogue of earthquake focal mechanisms, the data on brittle structures are introduced in a database. This raw assemblage of geometrical data forms a data set (or pattern if we follow the terminology of Angelier 1994). The raw data set is used as a starting point in stress inversion. For rock masses that have been affected by multiple tectonic events, the raw data set consists of several subsets (or systems) of brittle data. A subset is defined as a group of faults and fractures (in fault-slip analysis) or of focal mechanisms (in seismotectonic analysis) that moved during or have been generated by a distinct tectonic event. Moreover, the movement on all the weakness planes of the subset can be fully described in a mechanical point of view by the stress tensor characteristic of the tectonic event.
For an appropriate constraint on a stress tensor, data subsets should be composed of more than two families of data. Using the definition of Angelier (1994), a family is a group of brittle data of the same type and with common geometrical characteristics. In stress inversion, movement on a system of weakness planes is modelled by adjusting the four unknowns of the reduced stress tensor. Therefore, the stress tensor will be better constrained for data subsets with the largest amount of families of data of different type and orientation. This concept is used in the diversity criteria for the quality ranking procedure, described later.
In raw data sets, it is frequently observed that all the brittle data do not belong to a single subset as defined above. This is often due to the action of several tectonic events during the geological history of a rock mass. But the observed misfits can also be the consequence of other factors such as measurement errors, the presence of reactivated inherited faults, fault interaction, non-uniform stress field and non-coaxial deformation with internal block rotation (Dupin et al. 1993;Pollard et al. 1993;Angelier 1994;Nieto-Samaniego & Alaniz-Alvarez 1996;Twiss & Unruh 1998;Maerten 2000;Robert & Ganas 2000).
The frequent occurrence of multiple-event data sets and the numerous possible sources of misfits have important implications in fault-slip analysis and palaeostress reconstruction. It necessitates the separation of the raw data sets into subsets, each characterized by a different stress tensor. This check is necessary even in the case of a singleevent data set, as there are several possibilities to create outliers. Errors might occur during field work (e.g. uneven or bent fault planes, incorrect reading), during data input (incorrect transmission), or during data interpretation (e.g. not all events are detected due to the lack of data). A certain percentage of misfitting data (c. 10-15%) is normal, but they have to be eliminated from the data set for better accuracy of the calculated results.
In the iterative approach for stress tensor determination, data are excluded on the basis of a misfit parameter that is calculated for each fault or fracture as a function of the model parameters (o-1, o-2, o-3 and the stress ratio R) that best fit the entire set of data. A first stress model is determined on the raw data set, and then the data with the largest misfit are separated from the raw data set. After a first separation, this process is repeated and the original data set is progressively separated (split) into a subset containing data more or less compatible with the stress model calculated, and non-compatible data which remain in the raw data set. After the separation of a first subset from the raw data set, this process is repeated again on the remaining data of the raw data set, to eventually separate a second subset. We will discuss this procedure in more detail later when presenting the Right Dihedron and the Rotational Optimization methods.
In summary, data separation is performed during stress inversion as a function of misfits determined with reference to the stress model calculated on the data set. This is done in an interactive and iterative way and the two processes (data separation into subset and stress tensor optimization for that subset) are intimately related.
The first step of the selection procedure starts in the field. Faults or fractures of the same type, with the same morphology of fault surface, the same type of surface coating or fault gauge are likely to have been formed under the same geological and tectonic conditions. Already in the field, they can be tentatively classified into different families, and families associated into subsets.
Cross-cutting relations might help to differentiate between different families of faults, but it is not always easy to interpret these relations in terms of successive deformation stages. If relations between two types of structures are consistent and systematic at the outcrop scale (i.e. a family of normal faults systematically younger than a family of strike-slip faults), they can be used to differentiate fault families in the field and to establish their relative chronology. But this is not enough to differentiate brittle systems of different generation, as data subsets should be composed of several families of structures in order to provide good constraints during stress inversion. As much qualitative field information and as many observations as possible of the relations between pairs of fault or fractures are needed. But they are often insufficient alone to identify and differentiate homogeneous families of fault and fractures related to a single stress event.
The starting point for fault separation and stress inversion can be the first separation performed in the field. If this is not possible, a rapid analysis of the p and t axes associated with the faults and fractures can help differentiate between different families of faults and fractures, based on their kinematic style and orientation (e.g. normal, strike-slip and reverse faulting, tension and compressionat joints). This might be helpful if the measured faults and fractures belong to deformation events of markedly different kinematic styles. The separation done in the field has to be checked and refined during the inversion. In most cases, however, the selection of data relies mostly on an interactive separation during the inversion procedure.
However, we want to warn of pure automated data separation, because this might result in completely useless subsets. For better clarity we illustrate this risk with a 2D data set (depending on two parameters) instead of a 4D example (depending on four parameters), as it would be necessary for fault-slip data (three principal stress axes plus stress ratio). Figure 1 shows how the automated separation of a data set with two clusters results in three subsets, all of them representing only parts of the two clusters of the original data set. The same happens if multi-event fault-slip data sets are separated automatically. This problem can be handled by first making a rough 'manual' separation and then doing the first calculation. In the 2D example of Figure 1 this can easily be done by separating the data of the two clusters into different subsets. The manual separation of fault-slip data is not so straightforward because it is a 4D problem. The first requirement for a successful separation are field observations which indicate the existence of more than one event and can be used to discriminate the character of the different events. The best indicators for a multi-event deformation are differently oriented slip lines on one and the same fault plane which, in the best case, show different min- Separation of an artificial data set into subsets using an automated procedure. The original data set shows two clusters (a), but the calculated mean M~ lies in between them (b). Separation according to the deviation from M1 results in a data set compatible to M1 (c) and two other data sets with 'own' means M2 and M3 (d). Thus, the automated separation leads to three subsets (instead of two) with three different means Ml-M 3. The mean of the largest subset M1 (c) has nothing in common with the means of the two clusters in the original data set (a). eralization, so that the mineralization can then be used as a separation criterion. Additionally to field observations, obviously incompatible data can be separated through careful inspection of the data plots. As shown in Figure 2 some of the fault planes might have similar orientations, but different slip lines (e.g. down-dip versus strike-slip movements; see encircled data in Fig. 2). Separating these data into different subsets can serve as a starting point for the assignment of the other data, so that consistent subsets emerge (Fig. 2). These subsets are then checked for compatibility by feeding them into a computer method like the Right Dihedron and the Rotational Optimization implemented in the program TENSOR.
In this paper, we strongly advise making an initial separation in function of field criteria, a careful observation of the data plots and p-t analysis. The raw data set or the preliminary subsets should then be further separated into subsets while optimizing the stress tensor using successively the improved Right Dihedron method and the Rotational Optimization method in an interactive way.

Improved Right Dihedron method
The well-known Right Dihedron method was originally developed by Angelier & Mechler (1977) as a graphical method for the determination of the range of possible orientations of crl and o'3 stress axes in fault analysis. The original method was translated in a numeric form and implemented in different computer programs. We discuss here a series of improvements that we developed to widen the applicability of this method in palaeostress analysis. These new developments concern (1) the estimation of the stress ratio R, (2) the complementary use of tension and compression fractures, and (3) the application of a compatibility test for data selection and subset determination using a Counting Deviation.
Although the Improved Right Dihedron method will still remain downstream in the process of palaeostress reconstruction, it now provides a preliminary estimation of the stress ratio R and data selection. It is typically designed for building initial data subsets from the raw data set, and for making a first estimation of the four parameters of the reduced stress tensor. The Improved Right Dihedron method forms a separate module in the TENSOR program. The original method is described first briefly, before focusing on the improvements.

General principle The Right Dihedron method is
based on a reference grid of orientations (384 here) pre-determined in such a way that they appear as a rectangular grid on the stereonet in lower hemisphere Schmidt projection. For all fault-slip data, compressional and extensional quadrants are determined according to the orientation of the fault plane and the slip line ( Fig. 3a-c), and the sense of movement. These quadrants are plotted on the reference grid and all orientations of the grid falling in the extensional quadrants are given a counting value of 100% while those failing in the compressional quadrants are assigned 0%. This procedure is repeated for all fault-slip data ( . u e l f l l~d This method is particularly suitable for the stress analysis of earthquake focal mechanisms. For faultslip data, it gives only a preliminary result, as it does not verify the Coulomb criteria. Problems occur for data sets with only one orientation of fault planes (and with identical slip direction). In this case, 0-1 and o'3 deviate by 15 ° from the 'correct' position, if they are placed in the middle of the compressional/extensional quadrants. For conjugate faults the position in the middle of the 0%/100% area corresponds with o-1 and o-3 (Fig.  4). Moreover, it can be applied only when the sense of movement is given; otherwise the compressional and extensional quadrants are undetermined.
The graphical method gives a range of possible orientations of o-1 and o-3. These correspond to the grid orientations on the resulting counting net with respectively values of 0% and 100%. The mean orientations of these reference points give respectively the most probable orientations for 0-1 and o-3. Because 0-1 and o-3 are determined independently, they are not always perpendicular to each other. They can be set orthogonally by choosing either o'1 or o-3 fixed and rotating the other axis around a rotation axis defined as the normal to the plane containing 0-1 and 0-3. The orientation of the intermediate 0-2 axis can then easily be deduced.
One problem with this method is that it does not determine the stress ratio R (R = (0-2 -o-3)/(0-1 -0-3)) and that 0-1 and 0-3 are undefined when the extreme values on the counting net do not reach 0% and 100%.
Estimation of stress ratio R The stress ratio R, defined as equivalent to (o-2 -o-3)/(o-1 -o-3) is one of the four parameters determined in the stress inversion, with the three principal stress axes o-1, o-2 and o-3. Until now, the Right Dihedron method has given only an estimate of the orientations of o-1, o-2 and o-3, but not of the stress ratio R. A careful observation of the Right Dihedron counting nets, however, shows that their patterns differ as a function of the type of stress tensor (extensional, strike-slip or compressional). We therefore investigated a way to express this pattern as a function of a parameter that would be an estimation of the stress ratio R.
A way to do this is to compare the orientation of the previously determined o-2 axis with the distribution of counting values on the average counting net. Using the position of the o-2 axis on the counting grid, we found that a good estimation of the stress ratio R can be obtained with the empiric relation: where S2val is the counting value of the point on the reference grid nearest to the orientation of 0-2. This formula is only valid for large fault populations with a wide variety of fault plane orientations.
The accuracy of this method for the estimation of the R ratio has been validated using models with synthetic sets of faults obtained by applying different stress tensors on a set of pre-existing weakness planes of different orientation and computing the shear stress component ~-on the plane and the friction angle ~b. The differents sets are then submitted to stress analysis using the Improved Right Dihedron method. The values of the stress ratio R obtained are in general close to the ones used to produce the models, within a range of R _+ 0.1 (Fig. 5). Similarly, the orientation of the stress axes generally match within a few degrees the ones used to generate the synthetic sets. Experience gained using this method in conjunction with other methods of direct inversion (like the Rotational Optimization method described later) on a large number of sites shows that the stress ratio R esti- . Relation between stress ratios R used to produce models of synthetic fault sets by applying different stress tensors on a set of 152 pre-existing weakness planes of different orientation, and R values obtained by analysing these sets with the Improved Right Dihedron method. Models were produced with R = 0, 0.1, 0.3, 0.5, 0.7, 0.9 and 1.0. for extensional, strike-slip and compressional stress regimes. 0-1 located within the cone angle /3 and o-3 at an acute angle (-</3) to the tension plane. The orientations between the cone angle /3 and the surface generated by the revolution of a line inclined at an angle /3 from the fracture plane are considered as intermediate.
It is possible to define on the counting nets, areas in compression (value 0), in extension (value 100), and intermediate areas (value 50). The individual counting nets are summed up and averaged as in the case of faults, to obtain the average counting net. This allows the combined analysis of the three different types of data (slickensides with known sense of movement, tension and compression fractures, Fig. 3).
For the value of the angle/3 we use the common initial friction angle of 16.7 ° given by Byerlee (1978). The use of this value is justified by the shear stress/normal stress relations in the initial friction law (Jaeger 1969). If the angle between the tension fracture or the normal to the compression fracture and 0-1 is larger than 16.7 ° , the resolved shear stress is theoretically high enough to cause slip on that plane, which is in contradiction to the postulated nature of that plane. Other values for the angle /3 can be used, without modifying the general principle.
When working with a database composed only of fracture data, this method provides a way to estimate the stress ratio R, while this could not be determined by formal inversion. mated using the Improved Right Dihedron method is generally close to the one obtained by the Rotational Optimization method (see hereafter).
We conclude that for all types of stress tensors (extensional, strike-slip and compressional) with 0-1 > o-2 > 0-3 and 0.25 < R < 0.75, the Improved Right Dihedron method successfully estimates the four parameters of the stress tensor. In the extreme case of flattening (0-2 = 0-3, R = 0) or constriction (0-1 = 0-2, R = 1) only one of the extreme values will be well defined (0% for flattening, 100% for constriction); the other one will have medium values and a circular distribution. Therefore only 0-1 is well defined when R ~ 0 and o-3, when R 1.

Use of compression and tension fractures as palaeostress indicator
The Right Dihedron method also allows the use other types of brittle data such as compression and tension fractures as defined above for estimating the four parameters of the stress tensors (Fig. 3d, e). For tension fractures, 0-3 is considered as oriented within a cone angle of/3 degrees around the normal to the plane and o.1 is located at an acute angle (-</3) of the tension plane. The opposite is true for compression fractures, with Counting Deviation Another important improvement of the original Right Dihedron method is the development of a parameter for estimating the degree of compatibility of the individual counting nets with the average counting net of the subset. It relies on the calculation of a Counting Deviation CD (expressed in %) for each datum by comparing its counting net with the average counting net. Fault or fracture data with a low CD value contribute in a positive way to the average counting net (reinforce the extreme counting values) and the data with higher CD values contribute in a negative way (weaken the extrema). The principle of Counting Deviation and its use in compatibility testing are presented hereafter with reference to Figure 6. The first (and forward) step is to compute the average counting net by summing the values of each point on the reference grid for all the individual counting nets (Fig. 6a-c), taking into account the weighting factor associated with each datum (Fig. 6d). After this operation, the second step is performed in a reverse way. Each individual counting net is compared with the average counting net. As a result, a differential counting net is obtained for each datum by subtracting for each orientation of the reference grid the coun-  6. Use of the Counting Deviation CD to filter the data. As in Figure 3, data from three different types of brittle structure (a--e) are used to produce the average counting net (d), with consideration of the weighting factor. For all data, each counting value of the individual counting nets is multiplied by the weighting factor (4 for the normal fault and 1 for the other data). The stronger influence of the normal fault on the resulting average counting net is clearly visible (d). The differential counting nets (e) are obtained by subtracting respectively the counting values of the individual counting nets (a, b or c) from the ones of the average counting net (d), without application of the weighting factor. The Counting Deviation (CD) is the average counting value for the 384 points on the differential counting nets, expressed as a percentage (e). On the average counting net diagram (d), the CD values are displayed in a histogram, weighted according to the weighting factor associated with the individual data. The average counting deviation for the whole data set is 25 _+ 12.3% (1~). In this case, the reverse fault has a CD value above the average CD + 2cr and can be eliminated from the data set. If the process is repeated on the two remaining data, the dextral strike-slip fault will in turn be eliminated.

X4 X1 X1
D/Resulting average counting net and histogram of counting deviations • ~L~. N ting values of the individual nets from the those of the average counting net (Fig. 6e). The Counting Deviation (CD) is the average of all the 384 values on the differential counting net. As the counting values in the counting net are expressed as percentages, the unit for the Counting Deviation is also the per cent. The compatibility of individual data with the entire data subset is estimated by the dispersion of the CD values away from the arithmetic mean of all the CD values. The homogeneity of the data set is expressed by the standard deviation of the CD values: smaller standard deviations suggest that most data of the subset are compatible with the final counting net and hence possibly belong to the same kinematic event. In the TENSOR program, the fault planes can be rejected if their CD values are larger than the arithmetic mean +10-or +20, depending on the choice of the user. This is used to control the compatibility of each data with the final result, and to reject it from the subset if necessary.
Using this system, a first solution is computed on the initial data set (single-event deformation) or on the manually separated subsets (multi-event deformation), and the incompatible faults are rejected. The procedure is repeated several times until well-defined areas of 0 and 100 values are obtained in the counting net. (An example of progressive separation is given in Fig. 8a and Table 1 for fault-slip data set and in Fig. 9a and Table 3 for a focal mechanism data set.) The final result and the corresponding subset will serve as a starting point for the rotational optimization.
In summary, the Improved Right Dihedron method allows a first estimation of the orientations of the principal stress axes and of the stress ratio R, and a first filtering of compatible fault-slip data. The selected fault-slip population and the preliminary tensor can be used as a starting point in the iterative inversion procedures like the Rotational Optimization method described hereafter.

Rotational Optimization method
The Rotational Optimization method presented here is a new iterative inversion procedure. As for all iterative inversion methods, it is based on the testing of a great number of different stress tensors, with the aim of minimizing a misfit function. In principle the whole range of orientations for the three stress axes and the stress ratio R has to be tested to find the minimum value of the misfit function (grid search). Considering these four parameters, the problem is close to a four-dimensional one, with an additional constraint that the three principal stress axes have to be orthogonal. This leads to a large number of different configurations of the stress tensor to be tested with the data set.
To find the solution in a more direct way, we developed the Rotational Optimization method and propose to initiate the search procedure using the stress tensor estimated with the Right Dihedron method. It allows restriction of the search area during the inversion, so that the whole grid does not have to be searched.
In the following, we introduce a few classic theoretical notions, discuss in detail the misfit functions to minimize, and present the Rotational Optimization procedure with the help of a natural example from the Baikal Rift Zone in Siberia.

Shear stress construction and misfit parameters
The basic idea in the direct inversion is to test a number of stress tensor on all faults (and fractures) of the data set by computing the direction, sense and magnitude of the shear stress ('r) acting on the plane, the magnitude of the related normal stress (On), the characteristic dihedral angle 20 and the friction angle ~b.
In TENSOR, this is done using the method of Means (1989). In the inversion of fault-slip data, the isotropic part of the stress tensor is missing. Thus the reduced stress tensor is identical with the deviatoric part of the stress tensor. The magnitudes are expressed in a relative way because the absolute values cannot be determined using geological data only (see Angelier 1989). By convention, the magnitude of 0-1 is fixed at 100 and the magnitude of 0-3 at 0 (in abitrary units). The relative magnitudes are therefore in the range 0 = 0-3 -< 0-2 -< 0-1 = 100. The magnitude of the 0-2 axis is fixed by the stress ratio R, as a function of the magnitudes 0-1 and 0-3.
The occurrence of slip on a pre-existing rock discontinuity is governed by friction laws. On a Mohr diagram, the corresponding point is enclosed between the initial friction curve and the maximum friction line (failure curve). The relation between the minimum shear stress ('r) and the normal stress (v) is approximately linear (Angelier 1989). Hence, the sliding criteria can be simplified by assuming that the friction angle ~b must be greater than the initial friction angle and smaller than the maximum friction angle. The friction criteria of Byerlee (1978) are used here as default values (initial friction angle ~bo of 16.7 ° and maximum angle (~)max of 40.4°). If ~b < ~bo, no reactivation of pre-existing discontinuities will occur; if ~b > ~bmax failure will occur with the development of a new fracture.
Minimization functions A great advantage of the iterative approach for stress tensor inversion is that the complexity of the function to minimize is not a limiting factor and the function can be easily changed without changing the algorithm.
In general, the minimization function has the following form: where w(i) is the weight of the individual data and fj(i) is the function that has to be minimized.
In the following, by fault-slip data we understand fault planes (defined by the unit vector n normal to the slip plane) and associated slip vectors d on the slip planes, either observed or reconstructed as explained above. Faults with associated slip striae are also known as slickenside. In this category, we include also fault planes with slip direction only (direction of movement unknown). For those fanlt-slip data, the most classic misfit parameter is the angular deviation between the slip vector d and the shear stress t (slip deviation a), computed in function of the stress tensor and the orientation of the slip plane n:

fa(i) = a(i) (= function F1 in TENSOR program)
For slickensides with known sense of movement, the slip vector is uniquely defined and a ranges between 0 and 180 °. When the sense of movement is not known, only the orientation of the slip vector is defined, and two opposite directions are possible. In this case, the maximum slip deviation is at 90 ° from the slip direction. It is generally assumed that for a slip deviation of a --< 30 ° the observed fault movement is compatible with the theoretical shear vector.
It is also common to use a least-square function for the minimization. This implies a Gaussian-type distribution of individual misfits. A good example is function $4 Angelier (1991), already proposed by the same author in 1975: These two functions are still insufficient in the case of newly formed conjugate fault systems. In this case 0-1 and 0-3 might lie anywhere in the plane perpendicular to the two conjugate fault sets, because the slip deviation will be minimum for any configuration of stress axes. As long as o-2 is parallel to the intersection line of the two fault systems, the orientation of o-1 and o-3 will not influence the slip deviation. Thus additional constraints are necessary, taking into account the ability of the fault to slip (Angelier 1991). This can be either the friction angle ~b or the shear stress magnitude Irl on the fault plane, which both should be maximized. As the mechanical properties of the faulted rocks are generally unknown in standard palaeostress investigation, it is more appropriate to use I~1 to express the tendency of the fault to slip (see also Morris et al. 1996). Similarly, the normal stress component on the slip surface (v) also influences the ability of the fault to slip. It should be minimum in order to lower the fault friction and hence to favour slip. These three parameters can be taken into account simultaneously by combining in the same function the slip deviation a and the normal stress magnitude I vl that should be minimized, with the shear stress magnitude Irl that should be maximized. This is done in function f3(i): f3(i) = (f2(i) × 360) + ((Tinv(i) + Iv(i)l -29.7)/p) The first term is for the minimization of the deviation angle a and corresponds to the f2(i) defined above, multiplied by 360 to ensure that this term remains dominant with regards to the second term. With magnitudes of o1 and 0-3 fixed respectively at 0 and 100, the minimum possible value for the expression (Tiny(i) + Nnorm(i)) is 29.7 on a Mohr circle construction. To allow this term to converge to zero in the most favourable cases, 29.7 is substracted from the result of this expression. The proportionality factor p also allows the first term to be kept dominant with regard to the second term. In TENSOR, it is set by default to 5, but can be modified.
The TENSOR program can also use fracture data to constrain the stress tensors, and is not restricted to the analysis of slickenside data only. To implement this, other minimization functions have been developed.
For shear fractures (slip surface with observed movement but without slip lines), the second term of function f3(i) can be used as defined above: f4(i) = ((Tinv(i) + I~(i)l -29.7)/p) For tensional fractures (plume joints, mineralized veins, magmatic dykes), the normal stress Iv(i)[ applied to the fracture surface should be minimal to favour fracture opening. Simultaneously the shear stress magnitude It(i)[ should also be minimal to prevent slip on the plane: ~(i) = (Iv(i)l + I~(i)l)/p) For compressional fractures, the normal stress applied to the fracture should be maximal (0-1 -Iv(i)l), while the shear stress [r(i)l should be minimal: It is important to be aware that the stress ratio R can only be calculated with the functions that minimize the slip deviation a (fl-3(i)). This restriction has been implemented in the Rotational Optimization procedure in TENSOR. When working with compression or tension fractures only, a rough estimation of the stress ratio R can be obtained by the Right Dihedron method. This value is then maintained during the following procedures in the Rotational Optimization.
In program TENSOR, all functions described above should used for data sets containing only the corresponding data type. For mixed data sets (containing different types of data together) a composite function has been implemented. It allows the use of mixed data sets, by combining the optimization procedure for slickenside data and tension, shear or compression fractures. The contribution of each type of data is adapted in a way that they all have to be minimized to improve the quality of the tensor. In the general optimization function Fj, the individual functions f(i) are adapted to the type of data as follows: • for fault planes with slip lines: f3(i) • for shear fractures: f4(i) • for tension fractures: f5(i) • for compression fractures and stylolites: f6(i) This composite function (F5 in TENSOR) has been proved very efficient in palaeostress inversion of mixed data sets.
Interactive Rotational Optimization and kinematic separation of fault-slip data Stress tensor determination using the Rotational Optimization procedure consists in a controlled 4D grid search involving successive rotations of the tensor around the three principal stress axes (o-1, 0-2 and 0-3) and equivalent testing of the stress ratio R. For each stress axis, the rotation angle is determined for which the misfit function has its minimum value. The minimum value of the misfit function for each run is determined by taking the minimum of the polynomial regression curve adjusted to the results of each test, in a graphic representation with the rotation angle as abscissa and the value of the misfit function as ordinate (Fig. 7). The tensor is then rotated accordingly. The same procedure is repeated successively for the following stress axis. After rotation around the three stress axes, the value of R is determined in a similar way by testing a range of possible values of R. Each run involves first the adjustment of the stress axes, then the R ratio. The process is repeated several times until the tensor is stabilized, so that further rotations of the stress axes or modifications of the stress ratio do not improve the results.
The starting point for the first run is the results from the Improved Right Dihedron analysis (Fig.  8a). The tensor is rotated successively around each stress axis within a range of _+45 ° in steps of 22 to 5 ° and the full range of R values is checked in steps of 0.25 to 0.12 between 0 and 1. The configuration (o-1-3, R) which gives the lowest value for the optimization function is used as the starting point for the next run with smaller rotation angles and smaller steps for R. During the following runs, these values are progressively narrowed to -+-5 ° for the stress axes rotation and ---0.1 for the stress ratio check.
The palaeostress tensor obtained is defined for the population of fault-slip data for which it was computed. However, the raw population of faultslip data measured in the field is usually not homogeneous and not all fault-slip data can be attributed to a single stress tensor. This results from the fact that fault-slip data do not always fulfil Bott's (1959) basic assumption, which suffers from a series of limitations. The main limitations are: inhomogeneous stress field, pre-existing anisotropies, interaction between different faults or segments of a fault zone, asymmetrical stress tensor and thus rotation of the entire rock body or of internal blocks relative to the stress field (Dupin et al. 1993;Pollard et al. 1993;Twiss & Unruh 1998). Additionally, the fault pattern can be complicated by the existence of two or more subsequent deformation events. Therefore, palaeostress analysis involves the separation of fault-slip data into populations that can each be characterized by a unique stress tensor.
A maximum slip deviation a of 30 ° is set as upper limit for defining whether a fanlt-slip datum is compatible to a stress tensor. When using the composite function (function F5 in TENSOR), the value of the function is also used to check if fractures defined as shear, tensional or compressional fractures are compatible with the tensor. This function is defined in such a way that it reaches the value of 0 for best-fit situations and may reach a maximum of 20 for compressional and tensional fractures and 22 for shear fractures.
In the rotational optimization procedure, the separation is performed progressively during the inversion. This can be done after each optimization run rather than after the final determination of the tensor. To ensure an efficient fault separation, the Rotation angle Fig. 7. Principle of the 4D Rotational Optimization. The stress tensor is rotated successively around the o-l, 02 and o'3 axes by _ 12.5 and 45 °, a polynomial regression curve is computed by least square minimization using values of the optimization function (F5 here). The minimum of the function is found then additional rotations of ---5 ° from the angle corresponding to the minimum value of the function are again performed. A new regression curved is computed and the minimum is taken to define the rotation angle to apply to the initial tensor (+34 ° for the rotation around o-l). The equivalent rotation is performed and the procedure is repeated for the next two axes. A similar procedure is done for optimizing the stress ratio R, by checking a range of possible values between 0 and 1. The range of rotation angle sand values of R is progressively decreased to provide finer constraints during the next runs.
procedure of tensor optimization and fault separation is executed successively with a decreasing value of the maximum slip deviation (e.g. from 4 ° to 30°; Fig. 8 and Table 1). Because a stress tensor is adjusted on a particular fault population, a new stress tensor has to be recalculated after each modification of the fault population. Using this procedure, the first tensor obtained corresponds to the subset, which is represented in the original population by the greatest number of faults. The rejected faults are then submitted to the same procedure, to extract the next subset. Figure 8 and Table 1 present an example of progressive stress tensor optimization and fault separ-ation using successively the Right Dihedron and the Rotational Optimization methods on a natural fault-slip data set measured along a border fault of the Central Baikal basin in the area of Zama (Delvaux et al. 1997b(Delvaux et al. , 1999. This fault data set contains a minority of reverse or thrust faults related to an older brittle event, and a majority of normal or oblique-slip faults related to Late Cenozoic extension. We want to point out again that an uncontrolled automated separation presents a great risk because it might lead to useless results (see remarks above). Thus, we would recommend first doing a rough separation using field observations, careful inspec-  (Delvaux et al. 1997b(Delvaux et al. , 1999. The initial data base contains 54 fault-slip data, some related to an older compression stage but the majority of them related to Late Cenozoic extension. All stereograms have Schmidt lower hemisphere projections. The different parameters used for estimating the quality results are reported in Table 1. (a, 1 and 2) Initial data base; (3-6) progressive fault separation using the Right Dihedron, leading to a first subset and starting tensor for subsequent rotational optimization. (b, 7) After successive optimizations and elimination of incompatible data (8-12), the final solution is obtained for the first separated set, represented by the greatest number of fault-data (13 14). For the remaining 25 fault-slip data that were excluded from the initial data base (e, 15 and 16), a new series of separation was done using Right Dihedron (17), leading to a second subset which was again progressively optimized and separated (18 and 19) until the final solution is reached for the second set (20). Cross-cutting relations and fault plane morphology observed in the field show that the second set corresponds to an older stress state, compatible with the Early or Mid-Palaeozoic local stress field (Delvaux et al. 1995b) and that the first set is compatible with the Late Cenozoic stress field (Delvaux et al. 1997b).
tion of data plots, and then only starting the optimization procedure. At the end of the procedure, when two or more subsets are separated from the original fault population, it is necessary to test the stress tensors of each subset on the total fault population, to check if the fault separation has been done in an optimal way. Effectively, the faults-slip data are considered compatible with a stress tensor as soon as the deviation angle a is less than 30 °. In the process of stress tensor inversion and fault separation, a given fault-slip datum is attributed to the first population for which it is compatible with the characteristic stress tensor. But it can also be compatible with the stress tensors of the other subsets, sometimes with an even smaller deviation angle a. During the final cross-check (not shown in the example of Fig.   8), each fault datum is attributed to the population for which the composite function F5 is minimum. Because each subset was modified, a new optimization run has to be performed on each subset. As explained above, the deviation angle a alone is not enough to unambiguously discriminate between faults of different data subsets or between focal planes of the same focal mechanism, when they all fit the basic requirement that o~ < 30 ° and 0 > 0o. The criteria on the slip deviation c~ and the friction angle 0 are used to eliminate fault-slip data or focal planes from the data set, while the composite function is generally used as a minimization function during the inversion and for the finer discrimination between data subsets or focal planes.

Interpretation of results and quality assessment
The quality of the results obtained by stress inversion is dependent on a number of factors such as the number of data per subset, the slip deviation, the type of data, and even on the experience of the user. Thus, for a rating of the quality and a better comparability of the results it is necessary to introduce a quality ranking. A first quality ranking scheme was proposed by Delvaux et al. (1995b). It was later improved by Delvaux et al. (1997b).
In the first release of the World Stress Map (WSM; Zoback 1992), the quality ranking proposed for geological indicators was defined in a way that all stress tensors obtained from the inversion of Quaternary-age fault-slip data obtained the best quality. In co-operation with the WSM project (Sperner et al. 2003) and on the basis of these earlier approaches, the quality ranking scheme included in the TENSOR program was further This data set contains fault-slip data belonging to two different stress stages. The subset represented by the largest number of data will be extracted first. The CLw parameter progressively increases due to the elimination of data with less constrained slip sense determination. For the first subset, the values of Slen increases and limit the final quality of the tensor QRt to C quality, while the WSM quality QRw remains to B. This is justified by the fact that the retained fault population consists mainly in two conjugated sets of fault-slip data. As a result, the stress ratio R is not well constrained.
developed to meet the more strict requirements of the new release of the WSM project.
In accordance with the new ranking scheme for the WSM, the quality ranges from A (best) to E (worst), and is determined as a function of threshold values of a series of criteria. The threshold values have been chosen more or less arbitrarily and then tested with field data and checked to find whether they gave a good range of quality according to our 'feeling' of the results. A given quality is assigned if the threshold values corresponding to that particular rank are met for all the criteria.
For the WSM quality ranking scheme the following parameters are used, with threshold values as defined in Table 2: n, number of fault-slip data used in the inversion n/nt, ratio of fault-slip data used relative to the total number measured C~w, deviation between observed and theoretical slip directions Clw, slip sense confidence level for individual faults, applied for the TENSOR program: Certain As discussed in Sperner et al. (2003), the diversity of orientation of fault planes and of slip lineations is also an important parameter to consider. It is clear that the stress tensors are better constrained if the fault planes and the slip lines have a large variety of orientations rather than if they are all parallel to each other. The distribution of orientation data can be expressed by the normalized length of a vector obtained by addition of unit vectors representing the poles of the fault planes p or the slip direction s: where Px, Py, Pz and sx, Sy, s z are respectively the Cartesian co-ordinates of the unit vector p and s (expressed as direction cosines). The normalized lengths range between 1 for unimodal populations (all poles/slip lines parallel) to --+0.6 for homogeneously distributed populations.
Unfortunately the computation of such distribution has to be done on the original fault-slip data, which are generally not given in published results. For this reason, it was decided not to use these two criteria in the definition of the WSM quality rank for geological indicators (here named QRw, Table 2).
In the TENSOR program, the diversity criteria using Plen and Slen were also implemented and combined with the five first parameters, to determine the TENSOR quality rank QRt (Table 2). In this way, Plen and Slen are additional criteria which lower the quality of the result if the diversity of orientation of data is insufficient regarding to the WSM quality rank QRw.
A quality index (QI) expressed numerically has been previously proposed by Delvaux et al. (1995bDelvaux et al. ( , 1997b, based on the following formula: QI = n 2 X (n/nt)/OZw with threshold values of QI -> 1.5 for A quality, 1.5 > QI -> 0.5 for B quality, 0.5 > QI -> 0.3 for C quality and QI < 0.3 for D quality. However, this system proved to be unsatisfactory for a reliable quality assessment and we prefer to abandon it. Table 2. Threshold values as defined in Sperner et aL (2003)

Stress regime index
The stress regime can be expressed numerically by the stress regime index R' defined in Delvaux et al. (1997b). The main stress regime is a function of the orientation of the principal stress axes and the shape of the stress ellipsoid: extensional when 0-1 is vertical, strike-slip when o-2 is vertical and compressional when o'3 is vertical. For each of these three regimes, the value of the stress ratio R is fluctuating between 0 and 1. When the value is close to 0.5 (plane stress), the stress regimes are said 'pure' extensional/strike-slip/compressional. The transition between the three regimes is expressed by opposed values of R. An extensional regimes with R = 1 is equivalent to a strike-slip regime with R = 1. Similarly, a strike-slip regime with R = 0 is equivalent to a compressional regime with R --0.
To facilitate the representation of the range of stress regimes, Delvaux et al. (1997b) defined a stress regime index R' which expresses numerically the stress regime as follows: • R' = R when o-1 is vertical (extensional stress regime) • R' = 2 -R when 0-2 is vertical (strike-slip stress regime) • R' = 2 + R when 0-3 is vertical (compressional stress regime).
In the WSM, the naming of the stress regimes is correspondingly: normal faulting/strike-slip faulting / thrust faulting regime. The R' index forms a continuous progression from 0 (radial extension = flattening) to 3 (radial compression = constriction), while R is successively evolving from 0 to 1 in the extensional field, 1 to 0 in the strike-slip field and again 0 to 1 in the compressional field. R' has values of 0.5 for pure extension, 1.0 for extensional strike-slip, 1.5 for pure strike-slip, 2.0 for strikeslip compressional and 2.5 for pure compression. The intermediate stress regimes are sometimes also named 'transtension' for the transition between extension and strike-slip and 'transpression' for the transition between strike-slip and compression. The stress tensors are displayed in map view by symbols representing the orientation and relative magnitude of the horizontal stress axes, as suggested by Guiraud et al. (1989) and further developed by Delvaux et al. (1997b).
In regional studies it is sometimes necessary to obtain the mean regional stress directions of a series of closely related sites. In Delvaux et al. (1997b), we introduced the concept of a 'weighed mean tensor', without defining it clearly. This was based on the calculation of mean SHmax, Shmin and o-v stress axis by vectorial addition, taking into account the number of fault-slip data used for the stress inversion of each tensor. Similarly, the mean stress ratio was computed as the average stress ratio index R' defined above. In this procedure, the orientations of 0-1, 0-2 and 0-3 axes are assessed to SHmax, Shmin and Sv as a function of the stress regime (extensional, strike-slip or compressional), expressed as azimuth and plunge. Now, we prefer to use simply the average SHmax azimuth (as defined in the World Stress Map) and the average stress regime index R' as defined above. These two parameters describe fully the orientations of the stress axes (in terms of SHmax and Shmin, assuming Sv vertical), and the stress regime (R' = 0-1 for extensional regime, 1-2 for strike-slip regimes and 2-3 for compressional regimes).

Applications
Inversion of earthquake focal mechanism data Earthquake focal mechanisms are defined by two orthogonal nodal planes, one of them being the plane that accommodated the slip during seismic activation (fault plane) and the other being the auxiliary plane. In the absence of seismological or geological criteria, both nodal planes are potential slip planes and they cannot be discriminated. The Right Dihedron method is particularly well adapted for the stress inversion of focal mechanisms, as it also uses two orthogonal planes to define compressional and extensional quadrants. The Improved Right Dihedron method is useful for a first estimation of the stress tensor, and also for an initial separation of mechanisms from the database. The preliminary stress tensor and filtered focal mechanism data set are used as a starting point in the Rotational Optimization procedure ( Fig. 9 and Table 3).
In the Right Dihedron method, both nodal planes of an incompatible focal mechanism are eliminated simultaneously as they both have the same Counting Deviation values. In the Rotational Optimization procedure, all the nodal planes, which have slip deviations greater than the threshold value of, for example, 40 ° at the beginning of the analysis procedure or 30 ° at the end, are eliminated. The two nodal planes have generally different values for the slip deviation ~. If one of the planes has c~ greater than the threshold value, it is considered as an auxiliary plane and is eliminated. If both c~ values are higher than the threshold value, the entire focal mechanism is eliminated. If both o~ values are lower than the threshold value, the two nodal planes are kept for further processing. This results in the progressive selection of the probable fault plane for each focal mechanism.
The final choice of the presumed fault plane for  Fig. 9. Example of progressive kinematic separation and stress tensor optimization for earthquake focal mechanism data from the Central Baikal region (taken from Delvaux et al. 1999). All stereograms with Schmidt lower hemisphere projections. The different parmeters used for estimating the quality results are reported in Table 3. The initial data base is composed of 24 focal mechanisms, each expressed by two focal planes (a, 1 and 2). The Right Dihedron procedure first eliminates incompatible focal mechanisms (3-5) and produces a starting solution for the Rotational Optimization (5 and 6). The Rotational Optimization progressively improves the tensor and selects one focal plane for each mechanism on the basis of the value of the composite function, and eventually eliminates focal planes whose slip deviation a is >30 ° (b, 9-12).
each mechanism is done as a function of the value of the composite function during the rotational optimization. As explained above, this function combines the minimization of the misfit angle a and the maximization of the shear stress 171 on every plane. The composite function is therefore also convenient for the determination of the fault and auxiliary planes of focal mechanism if both have a slip deviation a lower than 30 °. The fault plane will be the one with the smallest value of the composite function. In this selection, any additional constraint (geological or seismological) on the fault plane has to be considered as dominant over the selection which minimizes the composite function. It should be kept in mind that the orientation of the activated fault plane can be influenced by pre-existing planes and might be more unfavourably oriented than a new plane. A special option has been implemented in the Rotational Optimization module of TENSOR that rejects automatically the non-compatible mechanisms and discriminates between the fault plane and the auxiliary plane of a focal mechanism, for any given stress tensor.
A similar (but less automated) procedure was used in Petit et al. (1996) for stress inversion of focal mechanisms from the Baikal Rift Zone. They compared their results with those obtained using the Carey-Gailardis & Mercier (1987) method on the same data set. It proved that, in most cases, the two methods give similar results. The differences result mainly from the different selections of fault and auxiliary planes, not from the inversion itself.

Tests with synthetic fault sets
To validate the results obtained with the Right Dihedron and the Rotational Optimization methods, we performed a series of tests using a synthetic set of fault and fracture data representing a strike-slip regime with north-south compression, vertical intermediate axis and R = 0.4. This was produced by applying the corresponding tensor on a set of 152 pre-defined weakness planes of various orientations distributed homogeneously on aster- Table 3. Value of the parameters used for estimating the quality of the results of progressive stress tensor optimisation and focal mechanism data separation using successively the Right Dihedron and Rotational Optimisation for the example in Figure 8 Right Dihedron Step For the Rotational Optimization, the procedure started on all focal planes (nt = 48) then one focal plane is selected for each mechanism at step 8 (nt = 24). In step 10, mechanisms previously excluded by the Right Dihedron procedure are reincorporated (n = 16-22). The highest possible quality for focal mechanisms is quality B because of the uncertainty in the differentiation between actual focal plane and auxiliary plane for each mechanism. The values of Plen and Slen increase with increasing separation, reflecting a decreasing in diversity of focal planes and associated slip directions. In this example the quality index (QI) regularly increases eonet ( Fig. 10-1). To simulate the different types of brittle structures, we used an initial friction angle of 10 °, a cohesion factor of 10 (arbitrary units) and a ratio o-3/o-1 of 0.2. On a Molar diagram ( Fig. 10-2), the individual data are plotted in different areas, corresponding to the types of structures generated: 88 faults, 9 tension fractures, 24 compression fractures and 31 shear fractures. When working separately with all the faults, tension fractures and compression fractures, the stress tensors obtained with both methods are always close to that used for generating the data set ( Fig.  10-3 to 10--6). The orientations of the stress axes always fit closely to those used for generating the data set. The R value obtained for the set of 88 faults is 0.51 with the Right Dihedron method and 0.4 for the Rotational Optimization method (exactly the original value). For the tension and compression fracture sets, the R value is estimated only with the Right Dihedron method and the same value is used in the Rotational Optimization. The R values obtained (0.75 for the set of 9 tension fractures and 0.35 for the set of 24 compression fractures) reflect the weaker constraint on R using fracture instead of fault data. The relatively small number of data in the tension fracture set might also influence the accuracy of the result.

Regional applications
The program TENSOR is now in use in more than 30 different laboratories worldwide and has been used in a large spectrum of geotectonic settings, for the late Palaeozoic to the Neotectonic period. These include the East African rift (Delvaux et al. 1992(Delvaux et al. , 1997aDelvaux 1993b,), the Baikal Rift Zone (Delvaux et al. 1995b(Delvaux et al. , 1997bPetit et al. 1996;San'kov et al. 1997), the Dead Sea Rift (Zain Eldeen et al. 2002), the Oslo graben (Heeremans et al. 1996), the Apennines (Cello et al. 1997;Ottria & Molli 2000), the Neogene evolution of Spain (Stapel et al. 1996;Huibergtse et al. 1998) and the Altai belt in Siberia (Delvaux et al. 1995a(Delvaux et al. , 1995cNovikov et al. 1998). Inversion of focal mechanisms has also been performed using this program (Petit et al. 1996).

Conclusions
New aspects of the tectonic stress inversion have been discussed. The separation of data into subsets during the inversion is an integral part of the stress analysis process. It has to be performed with much care and has to rely first on an initial 'manual' separation on the basis of careful field observations.  Not only faults with slip lines (slickensides) can be used in the stress inversion, but also fractures (tension, shear and compressional) and focal mechanisms. The improved Right Dihedron method allows (1) estimatation of the stress ratio R, (2) the use of tension and compression fractures in combination with slickenside data, and (3) an initial separation of the data to be performed on the basis of the Counting Deviation. The stress tensor and prelimi-nary data set obtained after the use of the Right Dihedron are used as a starting point for the Rotational Optimization. This method is based on a controlled grid search with rotational optimization of a range of misfit functions. Of the different misfit functions proposed, the composite function allows simultaneous minimization of the slip deviation a for slickensides, maximization of shear stress "r for fault planes and shear fractures, minimization of normal stress v for tension fractures and maximization of normal stress v for compression fractures and stylolites. A quality ranking scheme was developed in collaboration with the World Stress Map to estimate the quality of the results obtained. It is based on a series of parameters, for which threshold values are pre-defined. The type of stress regime is also expressed numerically by a stress regime index R'. The application of the Right Dihedron and Rotational Optimization is illustrated by a synthetic data set and natural data sets of fault-slip data and focal mechanism data from the Baikal Rift Zone.
All the aspects discussed here were implemented in the TENSOR program developed by the first author (available on request). This program is a tool for controlled interactive separation of faultslip or focal mechanism data and progressive stress tensor optimization using successively the Improved Right Dihedron method and the iterative Rotational Optimization method.