Camille Parlangeau
email: camille.parlangeau@gmail.com

Olivier Lacombe

Sylvie Schueller

Jean-Marc Daniel

Jean- Marc Daniel

Inversion of calcite twin data for paleostress orientations and magnitudes: a new technique tested and calibrated on numerically-generated and natural data

come

Inversion of calcite twin data for paleostress orientations and magnitudes: A new technique tested and calibrated

on numerically-generated and natural data

Introduction

Defining deformation mechanisms in the upper crust as well as their driving stresses are key scientific and technical issues. Quantifying stresses allows to better understand the mechanical behavior of geological materials and to decipher tectonic mechanisms, from those related to plate motions at a large scale to those causing jointing and faulting or even microstructures at a smaller scale.

However, the way stress orientations and magnitudes actually evolve in naturally deforming rocks over time scales of several tens of millions years is still a pending question, and even sophisticated geomechanical models require stress constraints in order to be accurately calibrated. The planes of twinning are the large circles; they contain the twinning direction [ei : rj] ; for each twin plane, the arrow is parallel to the twinning direction ; its head indicates that the upper part of the crystal moves upward, toward the C axis, as a reverse microfault (After Turner et al., 1954).

In order to decipher the tectonic evolution and to provide constraints on the past states of stress, methods of paleostress reconstructions based on the mechanical interpretation of various structural or petrographic elements in natural rocks have been set out. Among these methods, those based on the microstructural study of twinning in minerals have proven to be efficient in reconstructing stress (orientations and/or magnitudes) sustained by rocks during their tectonic history. These methods make use of the property that twinning depends on the magnitude of the shear stress which has been applied to them [START_REF] Tullis | The use of mechanical twinning in minerals as a measure of shear stress magnitudes[END_REF]. E-twinning is a common mechanism of plastic deformation in calcite aggregates deformed at low pressure and low temperature. E-twinning occurs in the host crystal by an approximation to simple shear in a particular sense and direction along specific crystallographic e planes (Fig. 1). Since calcite is widely encountered in sedimentary basins and mountain belts and is among the most sensitive mineral for twinning, it can be seen as an important paleostress indicator for the upper crust (e.g., [START_REF] Lacombe | Calcite Twins, a Tool for Tectonic Studies in Thrust Belts and Stable Orogenic Forelands[END_REF]. Inversion of calcite twins for paleostress is thus a well-established approach to reconstruct paleostress orientations and magnitudes of differential stresses [START_REF] Amrouch | Stress and strain patterns, kinematics and deformation mechanisms in a basement-cored anticline: Sheep Mountain Anticline[END_REF][START_REF] Arboit | Determination of the tectonic evolution from fractures faults, and calcite twins on the southwestern margin of the Indochina Block[END_REF][START_REF] Kulikowski | Combining geophysical data and calcite twin stress inversion to refine the tectonic history of subsurface and offshore provinces: A case study on the Cooper-Eromanga basin[END_REF][START_REF] Lacombe | Paleostress magnitudes associated with development of mountain belts: Insights from tectonic analyses of calcite twins in the Taiwan Foothills[END_REF][START_REF] Lacombe | From paleostresses to paleoburial in fold-thrust belts: Preliminary results from calcite twin analysis in the Outer Albanides[END_REF][START_REF] Lacombe | Joint analyses of calcite twins and fault slips as a key for deciphering polyphase tectonics : Burgundy as a case study[END_REF][START_REF] Lacombe | From paleostresses to paleoburial in fold-thrust belts: Preliminary results from calcite twin analysis in the Outer Albanides[END_REF][START_REF] Rocher | Cenozoic folding and faulting in the south Aquitaine Basin ( France ): Insights from combined structural and paleostress analyses[END_REF][START_REF] Rocher | Chanical twin sets in calcite as markers of recent collisional events in a fold-and-thrust belt : Evidence from the reefal limestones of southwestern Taïwan[END_REF] in nature and is a promising tool to estimate principal stress magnitudes when combined with fracture analysis and rock mechanics [START_REF] Amrouch | Paleostress magnitudes in folded sedimentary rocks[END_REF]Arboit et al., in press;[START_REF] Lacombe | Comparison of paleostress magnitudes from calcite twins with contemporary stress magnitudes and frictional sliding criteria in the continental crust: Mechanical implications[END_REF][START_REF] Lacombe | Paleostress magnitudes associated with development of mountain belts: Insights from tectonic analyses of calcite twins in the Taiwan Foothills[END_REF] or stylolite roughness paleopiezometry [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF]. Among the available techniques, the Calcite Stress Inversion Technique, CSIT (Table 1) [START_REF] Etchecopar | Etudes des états de contrainte en tectonique cassante et simulations de déformations plastiques (approche mathématiques[END_REF] allows the determination of the five parameters of the deviatoric stress tensor, i.e., principal stress orientations and differential stress magnitudes, the latter being strongly dependent on the (still debated) existence of a Critical Resolved Shear Stress (CRSS) for twinning [START_REF] Burkhard | Calcite twins, their geometry, appearance and significance as stress-strain markers and indicators of tectonic regime: a review[END_REF][START_REF] De Bresser | Strength characteristics of the r, f, and c slip systems in calcite[END_REF][START_REF] Newman | The influence of grain size and grain size distribution on methods for estimating paleostresses from twinning in carbonates[END_REF][START_REF] Rowe | Palaeostress estimation using calcite twinning: experimental calibration and application to nature[END_REF][START_REF] Tullis | The use of mechanical twinning in minerals as a measure of shear stress magnitudes[END_REF].

Despite a wealth of successful regional paleostress reconstructions in polyphase tectonic settings (see Lacombe, 2010, and references therein), the application by [START_REF] Gągała | Performance maps -a tool for examination of reliability of procedures for automatic separation of heterogeneous fault / slip and calcite twin data[END_REF] of the CSIT to numerically generated calcite aggregates has cast some doubt on the ability of the technique to efficiently separate superimposed tensors with close principal stress orientations. Moreover, methodological uncertainties in stress parameters as determined by CSIT have never been properly stated and remain to be better defined.

Table 1: List of symbols used.

The goal of this study is to present and test a new scheme of inversion, called CSIT-2, partly inspired by the CSIT. The newly developed technique and its ability to efficiently separate superimposed stress tensors are tested on numerically-generated and naturally deformed monophase (only one tensor applied) and polyphase (two tensors applied) twin datasets. The results demonstrate the reliability of such an approach to derive the stress parameters of interest, hence to provide an efficient toolbox for tectonic studies.

A brief review of existing methods to derive stress from calcite twins

Since the pioneering work of [START_REF] Turner | Nature and dynamic interpretation of deformation lamellae in calcite of three marbles[END_REF], several methods of stress reconstruction have been developed on the basis of the analysis of calcite twin data [START_REF] Etchecopar | An inverse problem in microtectonics for the determination of stress tensors from fault striation analysis[END_REF][START_REF] Jamison | Use of calcite twin lamellae to infer differential stress[END_REF][START_REF] Laurent | Determining deviatoric stress tensors from calcite twins: Applications to monophased synthetic and natural polycrystals[END_REF][START_REF] Laurent | Stress tensor determination from the study of e twins in calcite: A linear programming method[END_REF][START_REF] Nemcok | A stress inversion procedure for polyphase calcite twin and fault/slip data sets[END_REF][START_REF] Pfiffner | Determination of paleo-stress axes orientations from fault, twin and earthquake data[END_REF]Yamaji, 2015aYamaji, , 2015b)).

The basis of the widely used paleopiezometric method of [START_REF] Jamison | Use of calcite twin lamellae to infer differential stress[END_REF] is that, in a sample without any preferred crystallographic orientation, the relative percentages of grains twinned on 0, 1, 2 or 3 twin plane(s) depend on the applied differential stress value (σ₁ -σ₃). Since this relationship has been experimentally calibrated, knowing these relative percentages in a sample and under the hypothesis of a constant CRSS for twinning, the order of magnitude of (σ₁-σ₃) can be estimated. The severe limitations of this method are that (1) it does not take into account the grain size dependence of twinning, (2) it applies only to uniaxial stresses, (3) it does not check before calculation whether twin data are related to one or more superimposed stress tensors and (4) it does not make it possible to relate differential stress estimates to a given stress regime since principal stress orientations are not determined [START_REF] Lacombe | Calcite Twins, a Tool for Tectonic Studies in Thrust Belts and Stable Orogenic Forelands[END_REF]. [START_REF] Rowe | Palaeostress estimation using calcite twinning: experimental calibration and application to nature[END_REF] determined empiric laws between twinning incidence, twin volume fraction, twin density and experimentally applied differential stress to derive a calcite twinning paleopiezometer. The two first parameters are dependent on the grain size. The twin density is more easily used to calculate the differential stress because it does not depend on grain size. However, the law has been calibrated for temperature above 400°C and clearly overestimates the differential stress values when applied to samples deformed at low temperature (Ferrill,1998). In addition, this method shares the same limitations as the [START_REF] Jamison | Use of calcite twin lamellae to infer differential stress[END_REF] technique, which consists in not checking whether twin data are related to one or more superimposed stress tensors and not making possible to relate the differential stress estimates to a given state of stress. Laurent et al.'s technique (1981) was the first to take into account both twinned and untwinned planes to find 5 parameters of the stress tensor. This technique uses the deviatoric stress tensor (based on the nearly null dependence of twinning on isotropic stress) and relies upon a Boolean calculation to determine the solution tensor. This method was seldom used because it is time consuming ( Laurent et al., 1990). [START_REF] Laurent | Determining deviatoric stress tensors from calcite twins: Applications to monophased synthetic and natural polycrystals[END_REF] also proposed a new method which determines the 5 parameters of the deviatoric stress tensor. This technique uses a non-linear equation to determine the deviatoric stress tensor and is faster than [START_REF] Laurent | Stress tensor determination from the study of e twins in calcite: A linear programming method[END_REF] technique. These methods have not given rise to numerous applications in the literature.

The CSIT [START_REF] Etchecopar | Etudes des états de contrainte en tectonique cassante et simulations de déformations plastiques (approche mathématiques[END_REF] method is to date the most used technique to retrieve the past stress tensors (e.g., [START_REF] Amrouch | Stress and strain patterns, kinematics and deformation mechanisms in a basement-cored anticline: Sheep Mountain Anticline[END_REF]Beaudoin et al. 2012[START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF][START_REF] Lacombe | From paleostresses to paleoburial in fold-thrust belts: Preliminary results from calcite twin analysis in the Outer Albanides[END_REF][START_REF] Lacombe | Paleostress magnitudes associated with development of mountain belts: Insights from tectonic analyses of calcite twins in the Taiwan Foothills[END_REF][START_REF] Lacombe | Comparison of paleostress magnitudes from calcite twins with contemporary stress magnitudes and frictional sliding criteria in the continental crust: Mechanical implications[END_REF][START_REF] Lacombe | From paleostresses to paleoburial in fold-thrust belts: Preliminary results from calcite twin analysis in the Outer Albanides[END_REF][START_REF] Rocher | Chanical twin sets in calcite as markers of recent collisional events in a fold-and-thrust belt : Evidence from the reefal limestones of southwestern Taïwan[END_REF][START_REF] Rocher | Cenozoic folding and faulting in the south Aquitaine Basin ( France ): Insights from combined structural and paleostress analyses[END_REF]. The inversion process is very similar to the one used for fault slip data [START_REF] Etchecopar | Etudes des états de contrainte en tectonique cassante et simulations de déformations plastiques (approche mathématiques[END_REF], since twin gliding along the twinning direction within the twin plane is geometrically comparable to slip along a slickenside lineation within a fault plane. The basic assumptions are: (1) the existence of a Critical Resolved Shear Stress (CRSS) for twinning and (2) a potential twin plane is twinned (respectively, untwinned) if the resolved shear stress applied on it is greater (respectively, lower) than the CRSS. The inversion process takes into account both the twinned and untwinned planes and provides the 5 parameters of the deviatoric stress tensor, i.e., principal stress orientations and differential stress magnitudes. It should be noted that the value of the CRSS depends on grain size as well as on internal grain deformation since calcite hardens once twinned [START_REF] Newman | The influence of grain size and grain size distribution on methods for estimating paleostresses from twinning in carbonates[END_REF][START_REF] Tullis | The use of mechanical twinning in minerals as a measure of shear stress magnitudes[END_REF][START_REF] Turner | Experimental deformation of calcicte crystals[END_REF]. However the evolution of CRSS with grain size has not been well constrained yet. The CSIT technique has been successfully tested on experimentally deformed natural samples [START_REF] Lacombe | Determination of deviatoric stress tensors based on inversion of calcite twin data from experimentally deformed monophase samples: preliminary results[END_REF][START_REF] Laurent | Determination of deviatoric stress tensors based on inversion of calcite twin data from experimentally deformed monophase samples. Part II. Axial and triaxial stress experiments[END_REF] and the results have shown maximum deviations of computed principal stress orientations compared to experimentally applied ones of 5-7° for monophase and 7-11° for polyphase cases. In 2009, Gągała questioned the ability of this technique to reliably reconstruct the orientations of the principal stress axes and the stress ratio of close superimposed stress tensors. The penalization function used to refine the tensor solution is considered to be too restrictive by [START_REF] Rez | Peek inside the black box of calcite twinning paleostress analysis[END_REF] and Yamaji (2015b), especially when dealing with natural samples where optical measurement bias may occur. For [START_REF] Rez | Peek inside the black box of calcite twinning paleostress analysis[END_REF], the penalization function of CSIT (see section 3.1) is strongly dependent on the compatible twinned and the incompatible untwinned planes with the stress tensor solution. In addition, the space of solutions with a penalization function of 0 is too large. So, they proposed a new penalization function with sharper maxima, depending on the number of compatible twinned planes, the number of compatible untwinned planes and the number of incompatible untwinned planes. Note however that the CSIT with this refined penalization function has never been applied to polyphase twin datasets. [START_REF] Nemcok | A stress inversion procedure for polyphase calcite twin and fault/slip data sets[END_REF] technique works for both calcite twin and fault-slip data. The approach consists in searching all 3D tensors which can activate or re-activate a twin plane. They analyze clusters of twin data to group them in sets based on their response to one or multiple stress tensors.

Each subset is tested to establish monophase solution twin datasets. Polyphase sets are divided into monophase subsets. Each subset is analyzed to obtain the reduced stress tensor (orientation of principal stress axes and stress ratio). The study of [START_REF] Gągała | Performance maps -a tool for examination of reliability of procedures for automatic separation of heterogeneous fault / slip and calcite twin data[END_REF] demonstrates that this method is however not suitable for treating polyphase data.

The most recent stress inversion technique to date has been proposed by Yamaji (2015b). This technique, based on the generalized Hough transform, provides like CSIT, the 5 parameters of the deviatoric stress tensor. The preliminary exploration of the extent to which calcite twinning may constrain stress (Yamaji, 2015a) demonstrates that twinned planes better constrain the stress tensor than untwinned planes and that differential stress estimates are poorly resolved for differential stresses greater than 50-100 MPa. This technique seems to separate superimposed stress tensors only if the intersection between their spherical caps is nonexistent or small. Spherical caps correspond to the data points on a unit sphere based on the five-dimensional stress space; the deviatoric stress tensors responsible for twinning are thus denoted by the size and position of the spherical cap (see Yamaji, 2015a for details). Comparing his new technique to existing ones, Yamaji (2015a) states that the CSIT appears to be unstable, mainly due to the selection method of the solution tensors, which is only dependent on the last twinned plane taken into account in each solution tensor tested (see the theory of CSIT below). The limitation of Yamaji's technique is that it has not to date been tested on experimentally or naturally deformed samples.

An outcome of this short review is that the analysis of calcite twin data makes it possible to reliably obtain the 5 parameters of the deviatoric stress tensor under specific assumptions in monophase twin datasets. However, the challenge of all techniques is to reliably separate and reconstruct superimposed stress tensors as from polyphase datasets. There is especially a consensus toward the difficulty to detect and separate tensors which are very close in terms of principal stress orientations and differential stress magnitudes. These questions are addressed hereinafter through the setup and calibration of a new inversion scheme, CSIT-2. Its application to synthetic (numerically generated) calcite twin data allows to define the methodological uncertainties and the applicability domains of the new technique as well as the influence of heterogeneities commonly found in natural samples. The technique is finally applied to a naturally deformed polyphase sample to check its ability to reconstruct the regional paleostresses.

A new calcite twin inversion scheme (CSIT-2): methodology

3.1 Theory and basic equations of CSIT [START_REF] Etchecopar | Etudes des états de contrainte en tectonique cassante et simulations de déformations plastiques (approche mathématiques[END_REF] As said above, the principle of CSIT is to invert a calcite twin dataset for stress. The basic underlying hypothesis is that a potential e-twin plane is twinned if and only if the resolved shear stress applied on the direction of the twinning exceeds the critical resolved shear stress value 𝜏 𝑎 (CRSS):

 If 𝜏 𝑠 ≥ 𝜏 𝑎 , then the plane is twinned. (1)
 If 𝜏 𝑠 < 𝜏 𝑎 , then the plane remains untwinned. (2) with 𝜏 𝑠 , the resolved shear stress applied along the gliding direction of the e-plane. The optimally oriented twin plane will be activated if the applied differential stress is equal to, or greater than, 2𝜏 𝑎 .

The basic assumptions are that strain and stress are coaxial (low strain conditions), the stress field is homogeneous at the grain scale and twinning is a non-reversible process.

The principle of the inversion of calcite twin dataset is to find a stress tensor (or several stress tensors) which verifies these two inequalities for the largest number of twinned planes and the whole set of untwinned planes. The solution has the form of a reduced stress tensor with 4 parameters: the orientations of principal stress axes (σ₁, σ₂, σ₃) and the stress ratio (Φ):

1 ≥ 𝛷 = 𝜎 2 -𝜎 3 𝜎 1 -𝜎 3 ≥ 0 (3)
The differential stress (𝜎 1 -𝜎 3 ) is normalized to 1, so the normalized resolved shear stress applied on each twin plane varies within [-0.5 ; 0.5].

The inverse problem consists in finding the stress tensor that best explains the spatial distribution of measured twinned and untwinned planes. The first step consists in an arbitrary choice of a percentage of twinned planes to be explained. The resolved shear stresses are calculated on the twinned and untwinned planes, which are ranked as a function of the decreasing resolved shear stress.

In theory, the solution tensor should meet the requirement that all the twinned and untwinned planes should be consistent with it. Thus, all twinned planes should sustain a resolved shear stress (τs) larger than that exerted on all the untwinned planes. The sorting allows to determine rapidly whether some untwinned planes are incompatible with the tensor (i.e., the resulting resolved shear stress is greater than that for some compatible twinned planes).

The second step of the process therefore consists in calculating a penalization function, f, ideally equal to 0, which is defined as:

𝑓 = ∑(𝜏 𝑠 𝑗 -𝜏 𝑠 𝑚𝑖𝑛 ) 𝑗=𝑛 𝑗=1 (4)
where 𝜏 𝑠 𝑚𝑖𝑛 is the smallest resolved shear stress applied on the twinned planes compatible with the tensor, and 𝜏 𝑠 𝑗 is the resolved shear stresses applied on the j untwinned planes such that 𝜏 𝑠 𝑗 > 𝜏 𝑠 𝑚𝑖𝑛 .

The penalization function increases if incompatible untwinned planes are incorporated in the solution.

The optimal tensor is obtained when the maximum number of twinned planes and the minimum number of incompatible untwinned planes are incorporated in the solution.

This process therefore yields the orientation of the 3 principal stress axes, the stress ratio and a nondimensional differential stress. This non-dimensional differential stress,

(𝜎 1 -𝜎 3 ) 𝜏 𝑎
, is such that,

(𝜎 1 -𝜎 3 ) 𝜏 𝑎 = (𝜎 1 -𝜎 3 ) 𝑛𝑜𝑟𝑚𝑎𝑙𝑖𝑧𝑒𝑑 𝜏 𝑠 𝑚𝑖𝑛 = 1 𝜏 𝑠 𝑚𝑖𝑛
Under the assumption of a known and constant CRSS, 𝜏 𝑎 , the actual differential stress is given by equation 5:

(𝜎 1 -𝜎 3 ) = 𝜏 𝑎 𝜏 𝑠 𝑚𝑖𝑛 (5)

New inversion scheme (CSIT-2): motivation and principle

If the basic principle of CSIT-2 remains more or less the same as the CSIT for the calculation of the penalization function, the major difference lies in the automatic detection of one or several tensors. The choice has been made to systematically cover the space with a regular interval of 10° for the 3 Euler's angles (defining the orientation of the stress tensor axes) and a fixed stress ratio value of 0.5. The total amount of tested tensors is of 5832 tensors in CSIT-2. The tensor that is retained by the CSIT method for further optimization is the tensor presenting the minimum value of the penalization function. This means that if several tensors with the same minimum penalization function exist, only one will be arbitrarily chosen for the inversion. However, in case of several tensors with the same low penalization function (ideally 0, corresponding to a solution with incorporation of no incompatible untwinned planes), it is interesting to keep these different tensors, which might correspond to superimposed stress tensors. Contrary to CSIT, the CSIT-2 allows to automatically detect the presence of one or several tensors. For that purpose, the following steps are performed (Fig. 2):

(1) The user defines a starting percentage of twinned planes to be explained (20% recommended).

(2) The program will test systematically different tensors with the 3 Euler's angles being sampled every 10° and a stress ratio of 0.5.

(3) The penalization function is calculated for each stress tensor (equation 4).

(4) The stress tensors with a penalization function larger than a limit fixed by the user are discarded.

(5) Each tensor is weighted by the number of its nearest neighbors, based on the calculation of the angular distance [START_REF] Yamaji | Distances for the solutions of stress tensor inversion in relation to misfit angles that accompany the solutions[END_REF]. This step is important to automatically detect the different clusters of tensors with low penalization function values.

(6) The tensors are sorted in descending order based on the number of nearest neighbors and then in ascending order based on penalization function.

(7) The "best" tensors are selected: the first tensor of the list is the reference one. The reference tensor has the largest number of neighbors and the lowest penalization function value. The other tensors are kept if the percentage of shared twinned planes does not exceed a percentage threshold chosen by the user.

(8) The selected tensors are then optimized for different percentages of twinned planes to be explained using the [START_REF] Rosenbrock | An Automatic Method for Finding the Greatest or Least Value of a Function[END_REF] optimization method.

Thus, the choice of the retained tensors is based on different criteria: percentage of twinned planes to be explained (step 1), penalization function (steps 4 and 6), angular distance (step 5) and similarity criterion (step 7). The starting percentage of twinned planes to be explained is generally 20%. This value has been chosen because for a lower percentage, optical measurement uncertainties and grain-scale heterogeneities as expected in naturally deformed samples could generate too much noise and thus too many stress tensors could satisfactorily match the solution. If a percentage higher than 20% is chosen there is a risk to miss potentially superimposed tensors by only detecting the best expressed tensor. This percentage of twinned planes to explain has been chosen after several numerical tests to configure the different parameters depending on the user (the starting percentage of twinned planes to explain, the angular distance and the percentage of shared twinned planes between two tensors to determine if they are similar or not). Finally, in case of superimposed stress tensors with a high degree of similarity, it is possible that the inversion process may select an average tensor (i.e., resulting from averaging the two applied tensors) in addition to the two applied tensors.

However, this average tensor is easy to identify and to discard because the technique also yields the two applied tensors.

In order to automatically determine the different "best" tensors, it is necessary to highlight the clusters of tensors with low penalization function values. For that purpose, several steps are needed.

The angular distance is calculated between all the selected tensors using the angular stress distance defined by [START_REF] Yamaji | Distances for the solutions of stress tensor inversion in relation to misfit angles that accompany the solutions[END_REF]. Each tensor is weighed by its penalization function as well as by the number of tensors within an angular distance of less than 30° [START_REF] Yamaji | Distances for the solutions of stress tensor inversion in relation to misfit angles that accompany the solutions[END_REF]. This 30° value has been retained because it involves a low variation of orientation of the principal stress axes and a possible variation of the stress ratio of less than 0.5. This step is of first importance and makes it possible to determine clusters with the highest density of tensors and a low penalization function value.

The first tensor on the sorted list is used as the reference tensor. It explains at least the chosen 20% of twinned planes with possible incorporation of incompatible untwinned planes. The set of twinned planes explained by other tensors are compared to the set explained by this reference tensor.

The parameter of similarity is defined as the number of common twinned planes between the first and the second tensor divided by the total number of twinned planes explained by the second tensor. A new tensor is retained if it involves less than 70% of similarities (for detailed explanations, see subsection 3.6). This makes it possible to detect very close tensors that consequently share a high percentage of twinned planes. The twinned planes explained by this second tensor are recorded and each other tensor is compared in the same way to determine a possible third tensor, and so on.

Following this similarity criterion, few tensors are selected. Then they are further tested against the twin dataset by increasing the percentage of twinned planes to be explained. For each percentage, an optimization is carried out [START_REF] Rosenbrock | An Automatic Method for Finding the Greatest or Least Value of a Function[END_REF]. The ultimate goal is to explain the largest number of twinned planes and untwinned planes with an as low as possible penalization function value. The solution tensors are retained on the basis of the trend of evolution of their parameters: orientations of principal stress axes, penalization function, 𝜏 𝑠 𝑚𝑖𝑛 value, stress ratio and number of incompatible untwinned planes against the increasing number of twinned planes to be explained. The first criterion considered is the stability of the principal stress axis orientations. If they are stable despite the increasing number of twinned planes to be explained, the second criterion is the stability of the penalization function together with the evolution of the percentage of incorporated incompatible untwinned planes. This last parameter is defined as being the number of incompatible untwinned planes divided by the number of twin planes taken into account by the tensor (compatible twinned planes + incompatible untwinned planes). The solution then corresponds to the percentage of twinned planes to be explained that is marked by the break in slope of the penalization function value. If this percentage corresponds to the break in slope of the incorporation of incompatible untwinned planes, then the stability of the stress ratio and of the resolved shear stress value should be taken into account to ultimately refine the solution.

Figure 3: Graphics showing the variation of the penalization function value with the stress ratio for a same orientation stress axis (σ1: N180-0 and σ3: N270-0) and an applied stress ratio of 0.3 that has been registered after the added loop (see details in the first paragraph part 3.2) for A) 20%, B) 70% and C) 80% of explained twinned planes.

An example in Figure 3 illustrates how the best solution (maximum number of twinned planes and minimum number of incompatible untwinned planes) is defined. Increasing the percentage of twinned planes to be explained does not affect here the principal stress orientations which remain very stable (Fig. 3A &B), but it affects the 𝜏 𝑠 𝑚𝑖𝑛 value (Fig. 3E) and thus the differential stress value (equation 5). Figure 3C represents the evolution of the penalization function with the increasing percentage of twinned planes to be explained. This curve is very similar to the curve describing incorporation of incompatible untwinned planes (Fig. 3F). In this example the best choice is the solution at 88% of explained twinned planes because the slope of the penalization function curve greatly increases above this percentage. It means that the incorporation of incompatible untwinned planes becomes too high beyond this percentage. Figures 3D and3E display the evolution of the stress ratio value and of the 𝜏 𝑠 𝑚𝑖𝑛 value used to calculate the differential stress (equation 5). For the stress ratio curve the solutions between 20% and 36% of explained twinned planes are not stable.

Above 36% of explained twinned planes, the stress ratio is stabilized at about 0.5. Figure 3G represents the pseudo-Mohr circle showing τs against σn and figure 3H displays the evolution of τs for untwinned and twinned planes as a function of the percentage of twinned planes to be explained (calculated based on the total amount of untwinned planes). Both figures provide visualization of 𝜏 𝑠 𝑚𝑖𝑛 . As shown, the value of 𝜏 𝑠 𝑚𝑖𝑛 corresponds to a slight incorporation of untwinned planes in the solution (< 10%, Fig. 3H). In addition, these incompatible untwinned planes are plotted in figure 3G close to line representing the retained 𝜏 𝑠 𝑚𝑖𝑛 ; this means that these incompatible untwinned planes sustain a low resolved shear stress compared to most of twinned planes consistent with the tensor.

The solution is thus defined with confidence.

Figure 4: Determination of the best solution tensors of CSIT-based inversion on the basis of the evolution of parameters used during inversion process. The evolution of the parameters is shown as a function of the % of twinned planes to be explained by the tensor solution. The applied tensor is a N-S strike-slip regime (σ1= N180-0; σ3= N270-0) with a stress ratio of 0.5 and a differential stress of 35 MPa (test M9 in Table 1).

During the initial detection step, the choice to fix the stress ratio at 0.5 (following in that Yamaji, 2015b) was made (1) to reduce the computation time, and (2) because the inversion tends in polyphase cases to be stuck at tensors with extreme stress ratio (0 or 1) that account for a larger number of twin data. However, to make sure this choice does not bias the determination of stress ratio, a loop has been added before the optimization (step 8) which increases the stress ratio value by step of 0.1 from 0 to 1. The first test is a monophase case with σ₁ = N180-0; σ₃ = N270-0, a stress ratio of 0.3, a differential stress of 50 MPa (Fig. 4). It can be observed that the inversion is not very sensitive to the stress ratio at low percentage of twinned planes to explain. This suggests that at 20% of explained twinned planes (step 1), considering a fixed stress ratio (0.5) will not lead to artificially focus tensor solutions toward tensors with stress ratio of 0.5 hence to discard (step 4) those with stress ratios different from this value .

With this initial detection step, the determination of stress tensors using CSIT-2 requires only a short computation time, about 15-20 minutes.

Creation of synthetic calcite twin dataset

Table 2: Configuration of data tested for monophase cases.

The present method is first tested with numerically generated calcite twin datasets (Table 2 and 3). The use of synthetic data with various (but controlled) sources of complexity allows to better evaluate the performance, applicability domain and limitation of the new technique. Indeed, natural data have expectedly issues (spatial distribution of grains with different sizes, heterogeneity of stress within the aggregates due for instance to stress concentration at grain boundaries, inheritance -growth twins or earlier deformation twins -due to regional tectonic history), so that it is only after the synthetic experiments have been completed that the applicability of the CSIT-2 technique to natural data will be evaluated.

For the simplest case (monophase synthetic sample with perfectly homogeneous grain size, table 2), one tensor is applied on randomly generated twin planes within grains with random orientations of optical C axes (Fig. 5A). Each grain is assigned a CRSS value roughly simulating its virtual size since twinning is grain size dependent (i.e, twinning is easier in large grains than in small grains). A 5-parameter (deviatoric) stress tensor is applied on the generated twin data, resulting in some e-twin planes being activated (twinned) depending on their orientation with respect to the applied stress (Fig. 5B).

Table 3: Configuration of data tested for polyphase cases.

For the case of a polyphase synthetic sample (table 3) with homogeneous grain size, and in order to account for strain hardening, a second tensor is applied with each previously twinned grain being now assigned a CRSS increased by a value of 2 MPa as proposed by [START_REF] Gągała | Performance maps -a tool for examination of reliability of procedures for automatic separation of heterogeneous fault / slip and calcite twin data[END_REF].

Figure 5: Wulff lower hemisphere stereographic projection of the random distribution A) C axes and B) poles of twin planes in numerically-generated datasets. The applied stress tensor corresponds to: N-S σ1 axis, E-W σ3 axis, a stress ratio equal to 1 and differential stress equal to 75 MPa.

In order to simulate a heterogeneous grain size in both synthetic monophase and polyphase calcite aggregates, the choice was made to consider two virtual classes of grains with different sizes, simulated as two classes of grains with different assigned CRSS for twinning.

Finally, in order to simulate optical bias that depends on several factors (e.g., angle between the thin section and the twin lamella, width and spacing of twin lamellae, extension of lamellae across grains), some twinned planes were randomly misclassified into untwinned planes. For a horizontal thin section with 100 randomly oriented grains (i.e 300 e-twin planes), and considering that a twin plane lying at an angle of 30° or less to the thin section cannot be measured using a classical U-Stage hence can be misclassified as untwinned, we estimated that the percentage of potentially optically biased twin planes is always lower than 11% with an average of 6% (Fig. 6).

Figure 6: Evaluation of optical measurement bias related to inaccurate observations of twin planes. A) For a synthetic sample of 100 grains and 300 e-twin planes, and a horizontal thin section, the orange circle represents the cone of potentially biased poles to e-twin planes lying at less than 30° of the thin section plane. B) Variation of the percentage of potentially biased e-twin planes (lying at less than 30° of a virtual thin section plane cutting across a random distribution of twin planes) within various tests (100 synthetic datasets have been drawn and tested).

Following Yamaji (2015a, b), we further considered a maximum percentage of misclassified untwinned planes of 20-25% in order to also account for the other potential sources of heterogeneity at the scale of the virtual aggregate. Note however that the weight of this bias on the quality of the result of stress inversion may be highly variable. For instance, if a randomly misclassified untwinned plane is close to the optimal position for a potential twin plane to twin regarding the stress criterion, it will receive a high resolved shear stress. Then it will appear as an incompatible untwinned plane incorporated in the solution even for very low percentages of twinned planes to explain. This will greatly influence the penalization function and hence the quality of the solution. This abnormal plane is seen as an error by the software and the penalization function is strongly impacted, which could lead to discard the stress tensor from the acceptable solutions. It is necessary to efficiently compare the 5 parameter stress tensors resulting from the inversion (called "solution tensor") to the initially applied tensor (called "applied tensor"), in order to determine the accuracy of the inversion method. [START_REF] Peeters | Analysis of Distance/Similarity Measures for diffusion tensor imaging, Visualization and Processing of Tensor Fields[END_REF] propose different kinds of distance calculations between 5 parameter tensors. The L²-norm was chosen for our purpose, because it is the most robust for all tensor configurations:

𝐿² -𝑛𝑜𝑟𝑚 (𝐴, 𝐵) = √∑ ∑(𝐴 𝑖𝑗 -𝐵 𝑖𝑗 ) 2 3 𝑗=1 3 𝑖=1 2 (6)
with A and B the tensors with 5 parameters and L²-norm the value of the distance in mega Pascal (MPa). Note that this equation strongly resembles the one used by Yamaji (2015b, §3.5) and differs only by a factor of √2/3. The equation can be written as a measure between tensors, σ A and σ B :

𝐿² -𝑛𝑜𝑟𝑚(𝐴, 𝐵) = √(𝜎 𝐴 -𝜎 𝐵 ): (𝜎 𝐴 -𝜎 𝐵 ) = √ 2 3 𝑑 (7)
where the colon denotes the double-dot product of tensors [START_REF] Yamaji | Distances for the solutions of stress tensor inversion in relation to misfit angles that accompany the solutions[END_REF]Yamaji, 2007), and d is Yamaji's (2015b) dissimilarity measure.

Figure 7 helps visualize how the L²-norm value evolves with different configurations of tensors, hence its sensitivity to changes in σ₁ orientation, stress ratio and differential stress. The L²norm value is clearly more sensitive to variations of differential stress than to variations of the stress ratio or principal stress orientations.

Tests for calibration of the new inversion scheme

For the first tests, one tensor is applied (so-called monophase dataset), with one grain size.

The grain size is simulated using the corresponding CRSS value. But in nature it is very rare to have a sample with just one grain size. To fit with reality, other tests have been carried out with two distinct grain sizes, but also with a dispersion around the two mean grain sizes (pseudo-Gaussian distribution of grain sizes associated to different CRSS values also following a pseudo-Gaussian distribution).

Taking into account the 10 MPa CRSS value commonly adopted for grains of size of about 300 m and deformed at 3% (see [START_REF] Lacombe | Calcite Twins, a Tool for Tectonic Studies in Thrust Belts and Stable Orogenic Forelands[END_REF], it has been decided to test a distribution of grain sizes as corresponding to a distribution of CRSS values centered on 10 MPa, with 60% of grains with a CRSS value of 10 MPa, 15% between 10 -13 MPa and 7 -10 MPa, and 5% between 5 -7 MPa and 13 -15

MPa (denoted homogeneous grain size hereinafter). This enables to test the sensitivity and robustness of the technique to slight variations in CRSS value or grain size within a defined grain size class.

Biases are also incorporated in order to simulate the natural heterogeneities of a rock sample as well as optical measurement mistakes due to the use of a U-stage. To incorporate heterogeneities, the program will randomly change a percentage of twinned planes (fixed by the user) into untwinned planes.

In a second part, CSIT-2 is tested on datasets with two applied tensors and the same scheme is carried out in order to determine the uncertainties in the results.

Figure 8: Definition of the terms "deviation", "error", "accuracy" and "precision".

In order to quantify and discuss the results obtained using this new technique, several tests on each configuration have been carried out in order to calculate the average value for each parameter of the deviatoric tensor (orientations, stress ratio and differential stress, Table 2 and3). About 20-30 tests were necessary to get the stability of the solution for monophase datasets and 30-40 tests for polyphase datasets. We define the precision of the technique as the degree to which successive tests carried out under unchanged conditions yield the same results. The precision is quantified through evaluation of the variability of the determined stress parameter values around the mean values over the number of tests required to reach solution stability; it is called deviation hereinafter and is related to reproducibility of the results (Fig. 8). The lower the deviation, the higher the precision. We further define the accuracy of the technique by the degree of closeness of the mean reconstructed stress parameter to the true applied stress parameter. The accuracy is quantified by the error between the mean stress parameter value and to the true applied value (Fig. 8). The lower the error, the higher the accuracy. Both error and deviation will be given either in absolute value or in percentage of the applied value. As a result, we infer that the maximum methodological uncertainty associated with the determination of each stress parameter (i.e., when the inversion technique is applied blindly by a potential user to a naturally deformed sample) expectedly corresponds to ± (error + deviation) as derived from the results on numerically generated twin datasets. The error on the stress ratio will be calculated on the maximum value of the stress ratio, i.e., 1; so that we get the same error percentage for a given error whatever the applied stress ratio. For example, if the error value on the stress ratio is 0.1 for an applied stress ratio of 0.5 or 1, the error percentage will be of 10% in both cases.

Evaluation of the similarity of applied tensors

The inversion process is expected to hardly detect and separate close tensors, because the percentage of common twinned planes can be very high in this case. For instance, the examination of the activated twinned planes, which are shared by both tensors with similar stress ratio but with 1

axes at 30° to each other, shows that, for a homogeneous grain size dataset, the degree (percentage) of similarity increases with the differential stress applied and can almost reach 70% (Table . 4). In order to test the ability of the technique to detect and separate close superimposed stress tensors, two stress tensors with a deviation of 30° of the σ₁ axis were applied. In fact in tectonic studies, two tensors with less than 30° difference in principal stress orientations may be considered as similar because of cumulated inaccuracies on sample orientations and thin section making, measurement bias as well as possible local natural reorientation of stress tensors in rocks. We also tested superimposed stress tensors with perpendicular horizontal σ₁ axes (N-S and E-W) as well as stress tensors differing by permutations of σ2 and σ3 stresses (which is a frequent situation encountered in regional paleostress reconstructions).

Table 4: Evaluation of the degree of similarity between applied tensors. The table reports the maximum number of twin planes which are twinned by the two applied stress tensors for various stress configurations. These numbers are derived from 100 tests for each configurations. Φ is the stress ratio and Δ1-3 is the differential stress. Stress axis orientations are given in azimut/dip (°).

The similarity percentage is calculated by dividing the number of twinned planes shared by both tensors by the number of planes twinned by the second tensor. This parameter is not symmetrical:

comparing A with B or B with A. It has been decided to compare the second tensor applied with the first one because the first tensor can activate all the possible twin planes (depending on the equations 1 and 2) which is not the case for the second applied tensor (for which a strain hardening is applied).

The second tensor can therefore activate twinning on a lower number of twin planes than the first applied tensor. It is this reason why the similarity criterion is calculated based on the second applied tensor. The configuration with a permutation of σ2 and σ3 stress axes and the configuration with two tensors the orientation of the maximum principal stresses of which differ by 30° are supposed to be the most complex configurations in terms of tensor detection and separation. As mentioned before, the percentage of similarities can reach 69.6% at high differential stress (configuration 7 in table1). This is the reason why the critical percentage of similarities accepted by the process of tensor selection in the new CSIT-2 inversion scheme is 70% by default. We consider that this value makes it possible to differentiate most of the tensors belonging to different tectonic phases even in case of large differential stresses, hence large percentages of similarities.

Results

Monophase twin dataset

In the following part, the tests start with the simplest case, a monophase dataset with homogeneous grain size and no bias. Then, in order to simulate an increasing demand for the technique, this dataset is then made heterogeneous in terms of grain size, and will further incorporate virtual measurement bias. The basic stress configuration involves a N-S trending 1 axis and a E-W trending 3 axis.

Monophase twin dataset with homogeneous grain size

The first tests were carried out using datasets with one applied tensor and one grain size (CRSS = 10 MPa for all grains). A single orientation of principal stress axes was tested, but for this configuration, different stress ratios (0, 0.5 and 1) and differential stresses (25 to 75 MPa) were applied (configurations M1, M2, M3, M4, M5 and M6 in table 2).

Results from the inversion process show that the orientation of principal stress axes, stress ratio and differential stress are well reconstructed (Figure 9). A slight discrepancy is observed for differential stress magnitudes; the error however remains lower than 2.2 MPa (3%, Figure 9).

Figure 9: Results of inversion of numerically-generated monophase datasets with homogeneous grain size (configurations M1 to M6 in table 2). Φ is the stress ratio and Δ1-3 is the differential stress in MPa.

Total error is L²-norm. The error on the orientation is calculated using the scalar product between σ1 axes of applied and solution tensors.

Additional tests were carried out following the same scheme but with a slight dispersion around the prevailing grain size (Gaussian distribution around a CRSS of 10 MPa, configurations M9, M10 and M12 in table 2). Results show that for the error (L²-norm) between the applied tensors and the solution tensors increases with the applied differential stress but remains very low and does not exceed 6 MPa (Fig. 10). The error on the orientation of principal stress axes is on average lower than 4° with a deviation of ± 2° (Fig. 10). For the stress ratio and the differential stress values the error is close to zero with a very low deviation (Fig. 10). Datasets with a slight (pseudo-Gaussian type) dispersion within the two considered grain size classes are further considered (configurations M12, M13 and M14 in table 2). The error (L²-norm) is lower than 5 MPa but increases with the differential stress, reaching the value of 8 MPa. There is therefore a loss of information about the whole 5 parameters of the tensor around 75 MPa of differential stress (Fig. 12A). In details, regarding the principal stress orientations and stress ratio, they are well found, and the related errors do not increase with increasing applied differential stress (Fig. 12B andC). The error on the differential stress is of 0-1 MPa (1-3% of the applied value) below 50 MPa and of 3 MPa (4%) at 75 MPa (Fig. 11). The deviation is of ±1-2 MPa (4-5%) for differential stress below 50 MPa but increases for differential stress at 75 MPa (±7 MPa / 9%). Compared with the previous results dealing with a single grain size class, both deviations and errors are higher, especially on the differential stress value. This is mainly due to the grain size dependence of the CRSS that has not been considered here since a 10 MPa value for the CRSS has been blindly used to calculate differential stress irrespective of the heterogeneous grain size.

Figure 12: Inversion results of numerically-generated monophase datasets M12, M13 and M14 (Table 2). A, B, C and D same key as in figure 10.

Separating the datasets into virtual homogeneous grain size classes strongly improves the results as shown by the tests on datasets with homogeneous grain size. Thus, as discussed later, it is of first importance to carry out inversion only after the dataset has been divided into nearly homogeneous grain size classes, each corresponding to a given value of the CRSS. This is the reason why all the tests shown hereinafter are only performed on datasets corresponding to a homogeneous grain size (pseudo-Gaussian dispersion around a single CRSS value of 10 MPa).

Monophase twin dataset with homogeneous grain size and optical bias

In order to test the effect of an increasing bias, 10% or 25% of twinned planes have been randomly changed into untwinned planes. Results are presented in figure 13. Each point corresponds to the result of 20-30 tests. A stress ratio of 0.5 and a differential stress of 50 MPa are imposed (configurations M15 and M16 in table 2). The error (L²-norm) for 10% or 25% of incorporated bias are very similar (Fig. 13A) but with a larger deviation for the results with 25% of bias. If the parameters of the deviatoric stress tensor are considered one by one, the principal stress orientations, stress ratio and differential stress are very similar to those of the applied tensor. Thereby, the error on σ1 orientation is of 3° with an deviation of ± 2° (Fig. 13B), there is a negligible error (<0.01) for the stress ratio but an deviation of less than ± 0.1 (10%) (Fig 13C ), and a very slight error (< 2 MPa) on differential stress values with an deviation of about 10% or less (Fig. 13D).

Figure 13: Inversion results of numerically-generated monophase datasets M15 and M16 (Table 2). A, B, C and D same key as in figure 10.

Figure 13E and F show that with 10% and 25% bias, a solution tensor with a low error (L²norm) is found to the cost of an increase of the penalization function and also of the incorporation of incompatible untwinned planes with the increase of the incorporation of bias. For the tests with 25% of bias, the solution have been determined on the basis of the breaks in slope of curves as previously discussed, but the values for the penalization function and the percentage of incompatible untwinned planes in the solution are greater than the one for 10% or even 0% of incorporated bias. This means that compared to the 0% bias configuration M10 a larger tolerance threshold toward increasing penalization function value for 10% or 25% of bias explains the similar values for the L²-norm.

Polyphase calcite twin dataset

Polyphase twin dataset with homogeneous grain size, no bias and different stress ratio

Results are shown for polyphase cases with two applied tensors (configurations B18 then B21, B19 then B22 and B20 then B23 in table 3) and with a Gaussian distribution of grain sizes (centered on 10 MPa). The first applied tensor is a N-S compression with an E-W extensional part (σ₁: N180-0; σ₃: N270-0) with a stress ratio of 0.3. The second applied tensor is turned of 30° between σ1 axes (σ₁: N30-0; σ₃: N300-0) with a stress ratio of 0.7. These tensors have been applied on twin data for an increasing differential stress (35, 50 and 75 MPa) (Fig. 14A). The results of the CSIT-2 inversion on theses 3 configurations show that there is a loss of information with increasing differential stress (Fig. 14B). Figure 14C shows that the orientation of the two applied tensors are well recovered with an error of 4-7° (and a deviation of 2-3°) for both tensors. The error on the stress ratio for the first tensor is negligible (< 0.1) with a deviation of ± 0.03-0.05 (± 3-5%) (Fig. 14D). For the second tensor, the error on the stress ratio is of 0.06-0.13 (6-13%) with a deviation of ± 0.04-0.06 (± 4-6%). The differential stress is well retrieved with an error of 0.05-2 MPa (0.15-6.5%) with a deviation of ± 1-2 MPa (± 2.9-5.7%) for low applied differential stress (<50 MPa) and of ± 5-6 MPa (± 6-8%) for an applied differential stress of 75 MPa (Fig. 14E). These results show that CSIT-2 is able to identify and separate superimposed stress tensors with very different stress ratios. Figure 14D shows that with increasing differential stress there is no loss of accuracy but a loss in the precision of the results.

Figure 14: Inversion results of numerically-generated polyphase datasets B18 then B21, B19 then B22 and B20 then B23 configurations (Table 3). A) Wulff stereo-diagrams (lower hemisphere) of twin plane poles generated for the three configurations described previously. The graphs show the results of the inversion with B) the results of the L²-norm error and deviation between the applied tensor and the solution tensor, C) the error and deviation on the σ1 orientations, D) the deviation of stress ratio values with the mean value and E) the deviation differential stress values with the mean value. The red line represents the reference (applied) stress ratio and differential stress values.

Polyphase twin dataset with homogeneous grain size and optical bias

Results for polyphase cases with a difference of 30° between σ₁ axes of both tensors are summarized in figure 15 (configurations B1 then B3 and B2 then B4 in table 3). The distribution of grain size (i.e. CRSS value) is a pseudo-Gaussian, centered on a CRSS value of 10 MPa. The stress ratio and the differential stress are the same for both tensors, that is, Φ = 0.5 and Δ₁₋₃ = 50 MPa.

Figure 15B shows that, for the first tensor applied the error is of 9-11 MPa for the L²-norm with 10 or 25% of optical bias. Adding optical bias (10-25%) has almost no impact on the differential stress parameter. For the second tensor, the error (L²-norm) is about 10 MPa. The first tensor has an error of 6° on the orientation of principal stress axes with a deviation of ± 2-3° (Fig. 15C). The second tensors has an error of 7-8° on the principal stress orientations with a deviation of ± 2-3° (Fig. 15C).

The accuracy for the second tensor is slightly lower than for the first one. The error on the stress ratio parameter is of 0-0.03 (3%) with a deviation of ± 0.03-0.09 (3-9%) for both tensors (Fig. 15D). The error on the differential stress parameter is of 1-7 MPa (2-14%) with a deviation of ± 5 MPa (10%)

for both tensors (Fig. 15E). The different parameters of the stress tensors are fairly well estimated.

Figure 15: Inversion results of numerically-generated polyphase datasets applying B1 and B3 configurations then B2 and B4 configurations (Table 3). A) the results of the L²-norm error and deviation between the applied tensor and the solution tensor, B) the error and deviation on the σ1 orientations, C) the deviation of stress ratio values with the mean value and D) the deviation differential stress values with the mean value. The red line represents the reference (applied) stress ratio and differential stress values.

As for the monophase tests with the increasing percentage of bias incorporation, CSIT-2 shows similar results in terms of L²-norm error for 10% or 25% of incorporated bias. A more important incorporation of incompatible untwinned planes and an increase of the penalization function value are accepted. Thus a value of 10% of bias will only be tested hereinafter.

Figure 16 : Inversion results of numerically-generated polyphase datasets applying B5 then B7, B1 then B3 and B6 then B8 configurations (Table 3). A, B, C, D and E same key as in figure 14.

Having in mind the expected loss of information related to high applied differential stress, the influence of the differential stress on the result of the inversion process was then tested. The applied parameters are two stress tensors with 30° of difference in orientation of the σ₁ axes, a stress ratio of 0.5 and 10% of optical bias (configurations B5 then B7, B1 then B3 and B6 then B8 in table 3). Three differential stresses are tested: 35, 50 and 75 MPa. Results show that the error increases with the applied differential stress (Fig. 16B). The error on the orientation of σ1 (Fig. 16C) is of less than 6°

(with a deviation of ± 2°) for the first tensor and of 10° (with a deviation of ± 5°) for the second tensor.

The error on the stress ratio value is close to zero (< 10 -1 ) for both tensors and the deviation is of ± 0.1 (10%) (Fig. 16D). The error on the differential stress value is of 1-3 MPa (3-8%) for the first tensor and the deviation is of 2-10 MPa (6-13%) (Fig. 16E). For the second tensor with an applied differential stress lower or equal to 50 MPa the error is of 2-3 MPa (7%) with a deviation of ± 3-5 MPa (9%) (Fig. 16E). At 75 MPa, the error increases for the second tensor. This error is of -15 MPa (-20%) for a deviation of ± 13 MPa (21%) (Fig. 16E). Looking at the different stress parameters, the differential stress is the less constrained parameter. For the second tensor applied, there is a strong error for the last tests with a clear underestimate of the differential stress. This last observation concurs with the results of Yamaji's (2015b) study.

The influence of the differential stress on the results of the inversion process is also tested on the following cases with two tensors with σ₁ and σ₃ axes rotated by 90° (configurations B5 then B12, B1 then B13 and B6 then B14 in table 3). Figure 15 shows the L²-norm variation with various applied differential stress: 35, 50 and 75 MPa and a stress ratio of 0.5 for both tensors. For both tensors, a similar trend is observed with an increase of the average L²-norm value (the error) up to 75 MPa of applied differential stress and of the deviation (Fig. 17B). For principal stress orientations (Fig. 17C), the error remains the same for both tensors and even with increasing differential stresses (of about 2-4°) with a slight deviation of ± 1°. The error on the stress ratio is the lowest (Fig. 17D), as well as for the associated deviation. Figure 17E also indicates that inversion provides a very good estimate of the differential stress until 50 MPa. There is an increase of the error between 50 and 75 MPa for both tensors which can be estimated approximately at -6 MPa (-10 MPa for the second tensor) with a very slight deviation of ± 5 MPa (7%).

Figure 17: Inversion results of numerically-generated polyphase datasets applying B5 then B9, B1 then B10 and B6 then B11 configurations (Table 3). A, B, C, D and E same key as in figure 14.

Other tests have been carried out for polyphase datasets with superimposed stress tensors (Fig. 18A) with similar σ₁ axes, stress ratio (0.5) and differential stress (35, 50 or 75 MPa) but with a permutation of σ₂ and σ₃ axes (configurations B5 then B12, B1 then B13 and B6 then B14 in table 3).

The error (L²-norm) is the same for both tensors (less than 10 MPa) and is stable with increasing differential stress (Fig. 18B). The error on principal stress orientation is low (Fig. 18C), less than 5°

for both tensor and the deviation is about ± 2-3°. The error on the stress ratio (Fig. 18D) is decreasing with increasing differential stress and the same for the deviation. The differential stress is perfectly found with an error of 0 MPa and a deviation of less than ± 5 MPa (7%) (Fig. 18E).

Figure 18: Inversion results of numerically-generated polyphase datasets applying B5 then B12, B1 then B13 and B6 then B14 configurations (Tables 3). A, B, C, D and E same key as in figure 14.

Additional tests were carried out to evaluate the effect of a difference of stress ratio between both applied tensors (configurations B5 then B15, B1 then B16 and B6 then B17 in table 3) (Fig. 19A). The first tensor has a N-S compression with an E-W extensional part and a stress ratio of 0.

The second tensor has been rotated with respect to the first tensor by 30° and corresponds to a stress ratio of 0.5. Results show that the total error (L²-norm) is increasing with the applied differential stress and the same for the deviation (Fig. 19B). For principal stress orientation, the error is of 4-7°

for a deviation of ± 3-4° for both tensors (Fig. 19C). The stress ratio is very well reconstructed with a negligible error (< 10 -2 ) and a negligible deviation for the first tensor and an error of -0.1 (-10%) and a deviation of less than ± 0.1 (10%) for the second tensor (Fig. 19D). The differential stress is well retrieved by the inversion process but there is a slight increase of the deviation width for both tensors and a slight error for the second tensor with an applied differential stress of 75 MPa (Fig. 19E). The error is of 1-7 MPa (2-14%) with a deviation of ± 2-3 MPa (5-7%) for differential stress below or equal to 50 MPa and about 9-10 MPa (12-14%) of error with a deviation of ± 6-7 MPa (8-9%) for applied differential stress values of 75 MPa.

Figure 19: Inversion results of numerically-generated polyphase datasets of B5 then B15, B1 then B16 and B6 then B17 configurations (Table 3). A, B, C, D and E same key as in figure 14.

In the last tests the ability of the method to differentiate two tensors with a same value of stress ratio but varying differential stresses is tested (Fig. 20). Thus, the first tensor with a N-S trending 1

has a differential stress of 50 MPa and the second tensor with N030° trending 1 has a differential stress of 35 or 75 MPa (configurations B1 then B7 and B1 then B8 in table 3). Tensors are very close and these configurations are the most difficult in terms of separation because of the overlap of twin data explained by the two tensors. As previously, the distribution of grain sizes (i.e. CRSS value) is a pseudo-Gaussian centered on a CRSS value of 10 MPa (Fig. 20B &B'). Figure 20B &B' show the distribution of shared twinned planes of both tensors. It can be observed that there is a higher number of shared twinned planes in figure 20B' than in figure 20B where the differential stress of the second tensor equals 75 MPa. This is explained by the fact that both tensors are very close and an increase of the differential stress value increases the number of common twinned planes (progressively overlapping spherical caps in the sense of Yamaji, 2015a). Figure 20A & A' show that the principal stress orientations and stress ratio are well estimated but, as observed previously, there is an increase of error on the differential stress value for the second tensor at high differential stress (Fig. 20A').

For the configuration shown in figure 20A, the error is of 4-7° with a deviation of ± 2-3° for principal stress orientations, less than 0.1 (10%) for the stress ratio, 1-3 MPa (2-8%) with a deviation of ± 2-3 MPa (4-8%) for the differential stress. For the configuration shown in figure 20A', the error is of 6° with a deviation of ± 2° for principal stress orientations, 0.02 (2%) for the stress ratio and about 2-5.5

MPa (4-7 %) with a deviation of ± 3-5 MPa (6-7%) for the differential stress.

Figure 20 : Inversion results of numerically-generated polyphase datasets B1 then B7 and B1 then B8 configurations (Table 3). A and A') Results of the inversion for each configuration. B and B') Histograms of the repartition of grain sizes and an example of stereographic projection (Wulff lower hemisphere) of twinned/untwinned plane poles.

Interpretation and discussion of the results of inversion of numerically generated datasets using CSIT-2

5.1 Overall efficiency of CSIT-2 to determine paleostress from monophase and polyphase datasets

Determination of the principal stress orientations

For monophase datasets with a perfectly homogeneous grain size principal stress orientations are fully retrieved (Fig. 9). When a slight dispersion around the mean grain size is added, the average error is about 4° with a deviation of ± 2° (Fig. 10). The error increases with two distinct grain sizes of about 4° (by considering a constant CRSS value of 10 MPa) (Fig. 11). With two distinct grain size classes with a pseudo-gaussian distribution, the error on the principal stress orientations is the same than previously (Fig. 12). The same observations can be made for the results of dataset inversion with a bias incorporation of 10 or 25% (by randomly changing twinned planes into untwinned planes) (Fig. 13). A slightly larger deviation in the principal stress orientations is observed with an increase of bias incorporation in datasets (25%). Thus, CSIT-2 can accurately and precisely reconstruct principal stress orientations from monophase datasets in all configurations even with a high incorporation of bias.

As a result, for a naturally deformed monophase sample with homogeneous grain size, the methodological uncertainty associated with the determination of principal stress orientations does not exceed ± 8°.

For polyphase datasets with one grain size class, the determination of the principal stress axes is also accurate and precise, the maximum error reaching 10° on average with a maximum deviation of ± 5° (Fig. 14151617181920). This higher error compared to monophase datasets is due to the fact that CSIT-2 (as for CSIT) is able to estimate parameters of the stress tensor (orientations of principal stress axes and stress ratio) while the fifth parameter (non-dimensional differential stress) is found later by considering a constant CRSS value of 10 MPa. In the cases of datasets with different grain sizes, applied tensors may consequently not cause twinning of some favorably oriented twin planes. These planes therefore act as misclassified twinned planes (e.g., bias) during the inversion process. In fact, small grains have a higher CRSS value and the applied resolved shear stress value may thus be too low to cause twinning along them. This can slightly affect the determination of principal stress orientation.

The incorporation of bias in monophase and polyphase datasets makes it slightly more difficult to accurately and precisely determine principal stress orientations (Fig. 13 and15). This is due to the strong influence of the penalization function in the choice of the tested tensors. If the penalization function value exceeds 0.5, the tensor is not selected (see section 3.2). Depending on the orientation of untwinned planes with respect to the tensor, this value can exceed 0.5, leading to an error on principal stress orientations. The error on principal stress orientations is of 6° and 10° for the first and the second tensors when they differ of 30° in the orientation of principal stress axes with a deviation of ± 5° (for 10% bias).

Gagala ( 2009) stated that the CSIT was not able to separate very close tensors. In contrast, CSIT-2 makes it possible to unambiguously and reliably reconstruct superimposed tensors with only 30° of difference between σ₁ axes, hence with a high percentage of shared twinned planes.

The orientations of the principal stress axes are thus the best constrained parameters, because the 3D spatial orientation of twinned planes provides strong constraints on possible orientations of stress axes. This observation was the basis of the [START_REF] Turner | Nature and dynamic interpretation of deformation lamellae in calcite of three marbles[END_REF] analysis of calcite twins.

In nature it is very rare to have a sample having recorded just one phase (one applied tensor).

Thus, the user will be more interested by the accuracy for polyphase cases (two or more applied tensors). If the tensors are very close in principal stress orientations (i.e. 30° rotation between principal stress axes between both tensors), the error is 10° with a deviation of ± 5° (Fig. 141516). Thus, the maximum uncertainty in calculated principal stress orientations is ± 15°. For cases with very different applied tensors in term of orientations (i.e. around 90° rotation between principal stress axes between both tensors), the error is 4° with and deviation of ± 2°, hence the maximum uncertainty is ± 6° (Fig. 17). In the case of σ2 -σ3 stress permutation between both applied tensors the error is 4° with a deviation of ± 4°, so the maximum uncertainty is therefore ± 8° (Fig. 18).

As a result, the maximum methodological uncertainty associated with determination of principal stress orientations from a naturally deformed sample with an homogeneous grain size using CSIT-2 can be considered ± 15°; the usual uncertainty is expected to be about ± 10°.

Determination of the stress ratio

In general, the stress ratio is accurately and precisely retrieved in monophase and polyphase datasets, with bias of 10 or 25% and high differential stress (Fig. 13 and15).

For monophase cases, the error is negligible with a deviation of ± 0.04, so the maximum uncertainty for monophase cases is ± 0.04 (4%) on the stress ratio value. For polyphase cases with 30° of difference between principal stress axes of the two applied tensors, the error is of 0.03 (3%)

with and deviation of ± 0.08 (8%), so the maximum uncertainty is of ± 0.1 (10%). For tensors with 90° of difference between principal stress axes the error is of 0.02 (2%) with a deviation of ± 0.05 (5%), so the maximum uncertainty is of ± 0.07 (7%). For all these configurations CSIT-2 is very accurate and precise for the stress ratio determination. For tensors with a permutation of σ2 and σ3 axes, the error is expectedly higher, 0.3 (30%) for low differential stress and decrease until 0.07 (7%)

for high differential stress with a deviation of ± 0.1 (10%). This last configuration is the most difficult one because at low applied differential stress the number of shared twinned planes is very high and the program will tend to find one tensor with an extreme value of stress ratio. Extreme stress ratio values (in practice Φ = 0) make it possible to take into account more twinned planes for two tensors with similar σ₁ axes. With increasing differential stress, the percentage of shared twinned planes between tensors decreases and allows to more accurately and precisely find the stress ratio values.

As a result, the maximum methodological uncertainty associated with the determination of the stress ratio from a naturally deformed sample with an homogeneous grain size using CSIT-2 can be considered ± 0.4 (40%); the usual uncertainty is expected to be about ± 0.1 (10%).

Determination of differential stress magnitudes

In monophase datasets with one grain size class and no bias (Fig. 9-10), differential stresses are retrieved with a deviation that does not exceed 3 MPa. The loss of information on the determination of differential stress, when the applied stress is greater than 50 MPa as highlighted by Yamaji (2015a), is not observed (maximum error of 2 MPa for an applied differential of 75 MPa).

However, this effect is clearly observed with monophase datasets with two different grain size classes (see §5.3) or with polyphase datasets .

It is noticeable that, among the 5 parameters of the deviatoric stress tensor, the differential stress is clearly the less constrained parameter. This is mainly due to the sensitivity of the CRSS to grain size and to the difficulty to take this effect into account. The error may reach 20% in polyphase cases at high applied differential stresses. Results from monophase datasets show a good estimate of the differential stress compared to polyphase datasets. For differential stress lower than 50 MPa the error does not exceed 2 MPa for monophase datasets without bias. The results of the analysis of monophase datasets with two grain size classes point toward the need for dividing the dataset into subsets of homogeneous grain sizes before inversion in order to better estimate the differential stresses.

In polyphase cases with one grain size class and incorporation of 10% of bias, the maximum misfits on the differential stress value is observed at high applied differential stress (75 MPa). This occurs in all cases except for the configuration with a rotation of 90° (Fig. 18) of the principal stress axes. This configuration gives the most accurate and precise values for all parameters of the deviatoric stress tensors (regarding polyphase cases). For cases with two applied tensors with a difference of principal stress orientation of 30°, the error is of 1-6 MPa (12%) for applied differential stress equal or below 50 MPa and of 3-15 MPa (20%) above 50 MPa (Fig. 16). The deviation can reach ± 2-9

MPa (18%) in the worst cases for differential stress applied equal or below 50 MPa and ± 10-13 MPa (17%) above 50 MPa applied.

As a result, the maximum methodological uncertainty associated with determination of differential stresses from a naturally deformed sample with a homogeneous grain size using CSIT-2 can be considered ± 37%; this rather high value is expected for high differential stresses. The usual uncertainty is expected to be about ± 25-30%.

Influence of heterogeneous grain size

Several configurations of grain sizes have been tested in order to determine the robustness of CSIT-2 to varying grain sizes. As observed when dealing with monophase datasets, a high variability of grain sizes does not affect the accuracy of the determination of principal stress orientations or stress ratio. However, and especially in the case of two grain size classes , there is a loss of information for differential stresses. This loss is higher at high differential stress. This is due to the application of a constant CRSS value of 10 MPa for all grains. As mentioned previously, several authors concurred to say that the CRSS value is strongly gain size dependent [START_REF] Newman | The influence of grain size and grain size distribution on methods for estimating paleostresses from twinning in carbonates[END_REF][START_REF] Rowe | Palaeostress estimation using calcite twinning: experimental calibration and application to nature[END_REF]. Thus, the CRSS value has to be adapted for each grain size. The main problem is our poor knowledge regarding the evolution of the CRSS value with grain size.

Principal stress orientations are generally accurately and precisely determined by CSIT-2, although there are few cases at low applied differential stress showing bad results (mainly on the dip value of principal stress axis). This is due to the fact that the software calculates the 4 parameters using (σ₁ -σ₃) = 1 without taking into account the different grain sizes. Thus, at low differential stress, twinning in small grains (with high CRSS values) cannot be activated due to the differential stress of the tensor even if they are favorably oriented regarding the stress tensor orientation. In these cases, having two different grain sizes can bias the results.

Influence of optical measurement bias

The influence of the incorporation of bias in the datasets appears to be complex. It depends on the orientation of misclassified planes regarding to the applied tensor. If a randomly misclassified plane has a low applied resolved shear stress value, then the effect of the misclassification on the results is insignificant. If the randomly misclassified plane has a high applied resolved shear stress value regarding the applied tensor, then there is a strong influence on the penalization function, which consequently increases. Because the selection of tested tensors relies partly on the penalization function value, if this value is too high (above the limit fixed by the user), the tensor will be discarded, even if this tolerance toward a value of 0.5 allows to be less restrictive than the previous method (CSIT). The CSIT-2 selection process may miss "good" tensors depending on the orientation of the misclassified untwinned planes. It is to date impossible to predict how it affects the orientations of principal stress axes and other parameters because it is depending on the resolved shear stress applied on this incompatible untwinned planes and how many there are. Few tests tend to show that the principal stress orientations can be strongly affected, especially the dip of the stress axes, the azimuths being well retrieved. The stress ratio and the differential stress are not affected by the inaccuracy of the principal stress orientations. They remain well retrieved.

Our results show that the CSIT-2 method is overall robust for twin datasets with up to 25% of incorporated bias. Considering the breaks in slope on the curves of the different parameters as a function of the increasing number of twinned planes explained (Fig. 3 and 13F) leads to accurate results even at such high bias (see sections 3.2 and 4 for more details), but at the cost of an increase of the penalization function (up to 1.5) and also of the incorporation of incompatible untwinned planes (number of incompatible untwinned planes over number of compatible twinned planes + number of incompatible untwinned planes, Fig. 13E, up to 1.5). However, for natural cases the amount of optical bias is not a priori known and an untwinned plane seen as incompatible for the inversion is not necessary an optical measurement bias. It is therefore recommended, when analyzing (blindly) a naturally deformed sample, to consider with care any solution tensor that would incorporate (more than) 30% of incompatible untwinned planes and a penalization function value that would exceed 1.5 (Fig. 13E). The geological setting helps determine whether the solution-tensor is significant or not.

Practical way to optimize data inversion

The optical bias, likely due to errors in measurements using a U-stage can have a significant effect (on the orientation of principal stress axes) despite the improvement of the technique. To circumvent this difficulty, an index of uncertainty in the untwinned status of a twin plane has been established for measurements of natural samples. This way, the user can double-check the status and, in case of doubt, even change the status of the plane of poor confidence in the dataset. Note however that this will be applied only on untwinned planes which appeared at low percentage of twinned planes to be explained and therefore on which the tensor exerts a high resolved shear stress value.

In order to minimize the optical bias and measurement errors and to get information from completely untwinned grains, twin data are currently being measured using EBSD rather than using a U-stage [START_REF] Parlangeau | Inversion of calcite twin data for paleostress (2): EBSD as a tool for data measurements[END_REF]. Optical bias consequently falls down to 0%. However, natural heterogeneities at the grain scale (stress shadows, heterogeneous grain-to-grain stress transmission) can of course not be circumvented and their weight on the quality of results can only be estimated by comparing natural (non-perfect) data and numerical data. The decrease of incorporated bias has a positive effect on the accuracy of the stress results (see comparison of monophase results without/with bias in section 5.3).

Overall, the least constrained parameter of the reconstructed deviatoric stress tensors is the differential stress magnitude, especially when applied differential stresses are greater than 75 MPa.

Combining inversion of calcite twin data with fracture analysis and rock experiments (e.g., [START_REF] Lacombe | Comparison of paleostress magnitudes from calcite twins with contemporary stress magnitudes and frictional sliding criteria in the continental crust: Mechanical implications[END_REF] may however help test the validity and the meaning of determined stress magnitudes against the mechanics of rock masses. Our study finally highlights that sorting the whole dataset into subclasses of homogeneous grain size before inversion is the best way to reliably define the 5 stress tensor parameters, especially differential stresses. That is the reason why it is strongly recommended to estimate the grain size distribution in a natural twin datasets and to apply separately the inversion technique to homogeneous grain size classes, provided the number of available twin is sufficient to secure a reliable stress tensor computation. In that case, provided each class yields similar accurate principal stress orientations and stress ratio, assigning a different value of the critical resolved shear stress to each class as a function of the mean grain size should provide more reliable (and expectedly consistent) estimates of differential stresses (see section 7). Taking into account the grain-size dependence of the CRSS in the inversion process is therefore a strong requirement; further work should focus on providing better constraints on the variation of the CRSS as a function of grain size (work in progress).

Comparison with CSIT

Few studies have focused on estimating uncertainties on the results of the former CSIT (Lacombe andLaurent, 1996, 1992;[START_REF] Laurent | Determination of deviatoric stress tensors based on inversion of calcite twin data from experimentally deformed monophase samples. Part II. Axial and triaxial stress experiments[END_REF]. In order to illustrate the improvement on stress determination brought by using CSIT-2 instead of CSIT, we provide hereafter a comparison of the stress results obtained from both techniques when applied to polyphase datasets since [START_REF] Ga̧gała | Reliability of selected procedures of stress inversion and data separation for inhomogeneous populations of calcite twins and striated faults: Insights from numerical experiments[END_REF] has casts some doubt on the ability of CSIT to well retrieve close superimposed stress tensors. The first configuration tested is the one applying a N-S compression with a stress ratio of 0.5 and a differential stress of 50 MPa, the second applied tensor is turner of 30° between σ1 axes with a stress ratio of 0.5 and a differential stress of 50 MPa (configurations B25 then B27 in table 3). The second configuration is exactly the same than the first one but with a differential stress of 75 MPa for both tensors (configurations B26 then B28 in table 3). The last tested configuration is applying a first tensor with a N-S strike-slip regime (σ₁: N180-0 and σ₃: N270-0) with a stress ratio of 0.5 and a differential stress of 35 MPa. The second tensor applied is also a N-S strike-slip regime but with a permutation of the σ₂ and σ₃ axes (σ₁: N180-0 and σ₃: N270-90) with a stress ratio of 0.5 and a differential stress of 35 MPa (configurations B24 then B29 in table 3).

Table 5: Results of the inversion using CSIT and CSIT-2 of the configuration B25 then B27. The L²norm is the distance between the applied tensor and the solution-tensor found by CSIT or CSIT-2.

For the first configurations B25 then B27 (Table 3), the results are shown in table 5. L²-norm shows that, again, CSIT-2 is more accurate than CSIT. CSIT finds a tensor which allows the permutation between σ₃ and σ₂ axes with a stress ratio very low of 0.03.

Table 7: Results of the inversion using CSIT and CSIT-2 of the configuration B24 then B29. The L²norm is the distance between the applied tensor and the solution-tensor found by CSIT or CSIT-2.

As a result, although some stress configurations are similarly retrieved by CSIT and CSIT-2.

CSIT-2 being always more accurate than CSIT. Some other more complex stress configurations are clearly badly retrieved by CSIT while being consistently retrieved by CSIT-2. This illustrates that using CSIT-2 will expectedly improve the reliably of stress determination in polyphase settings in forthcoming tectonic studies.

Application to naturally deformed polyphase samples

In order to test the CSIT-2 and its ability to derive geologically meaningful paleostresses, the technique is applied on samples collected from limestones of Tithonian-Aptian age (Maiolica formation) in the Monte Nero anticline located within the Umbria-Marches tectonic domain of the northern Apennines (Italy). The Monte Nero-Serra Santa is an arcuate fold with a backlimb dipping 30° towards W-SW and an intensively folded forelimb, where the strata are locally overturned (Fig. 21). Microstructural investigations revealed occurrence of several sets of joints and veins (mode I opening) that developed mainly before the onset of folding [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF]. 854 In the western limb of the fold, we selected a first sample containing two intersecting veins 855 belonging to successive sets recognized in the field. After unfolding, Vein 1 (Fig. 22A) belongs to a bed-perpendicular N050 vein set (N65-69S in present attitude), while Vein 2 belongs to a second set of bed-perpendicular veins striking N090-110° (N291-84N in present attitude). Vein 2 clearly cuts across Vein 1 and therefore unambiguously postdates it.

Figure 22: Inversion results of naturally deformed polyphase datasets. A) Picture of a thin section (crossed polars) showing crosscutting veins, together with stereographic projection of measured veins (Wulff's lower hemisphere equal area projection). B) Plots of C axes, and poles of twinned and untwinned planes from veins 1 and 2 together with distribution of related grain sizes. C) Results of inversion of calcite twins from vein 1 and vein 2 (for more details, see text). For the vein 1, the first tensor is oriented with σ₁ : N63-2 and σ₃ : 297-87 and the second tensor is oriented with σ₁ : N110-21 and σ₃ : 15-13. For the vein 2 small grains, the first tensor is oriented with σ₁ : N105-1 and σ₃ : 195-3 and the second tensor is oriented with σ₁ : N108-16 and σ₃ : 215-47. For the vein 2 big grains, the first tensor is oriented with σ₁ : N95-1 and 1227 σ₃ : 185-20 and the second tensor is oriented with σ₁ : N300-8 and σ₃ : 200-53. This chronology between the two vein sets is documented in many places at the scale of the fold and Inversion results shown in figure 22C illustrate that two tensors are found by the inversion process within each vein. Vein 1 records a N050 ± 6° compression with a stress ratio of 0.15, a differential stress of 63 MPa ± 3-6% and a N110 ± 6° compression associated with perpendicular extension, a stress ratio of 0.62 ± 0.1, a differential stress of 35 MPa ± 8%. Vein 2 yields a stress tensor similar to the second tensor from vein 1. The uncertainties associated with the previous results are estimated from the interpretation of the analyses on synthetic datasets. Results from the inversion process are consistent across grain size ranges.

A relative chronology between successive twinning stresses can be established through consideration of orientation of computed stress axes (as well as of stress ratio) with respect to vein orientation. The simple underlying idea is that a stress tensor determined from the calcite grains filling a vein with a 3 axis perpendicular to the vein strike is likely related to the vein opening while other tensors with stress axes inconsistent with the vein geometry likely reflect later, post-opening stress regimes [START_REF] Lacombe | Calcite Twins, a Tool for Tectonic Studies in Thrust Belts and Stable Orogenic Forelands[END_REF].

The compressional configuration of tensor 1 from Vein 1 with subhorizontal 1 axis lying parallel to the vein and subvertical 3 axis, is not at first glance consistent with opening of Vein 1.

But the low value of the stress ratio supports that 2 and 3 are very close in magnitudes hence prone to switch, so the computed 3 axis is consistent with vein opening in a strike-slip / compressional stress regime. This tensor should be considered as predating tensor 2, which is therefore likely postopening. The strike-slip configuration of tensor 2 from Vein 1 with subhorizontal 1 and 3 axes lying within and perpendicular to Vein 2, respectively, is in perfect agreement with opening of this vein. This means that Vein 1 recorded the stress regime responsible for opening of Vein 2 which postdates Vein 1. Similarly, the strike-slip configuration of tensor 1 from Vein 2 with subhorizontal 1 and 3 axes lying within and perpendicular to the vein, respectively, is in perfect agreement with opening of Vein 2, and should be considered as predating tensor 2 from Vein 2. Tensor 2 from Vein 2 is therefore a likely post-opening stress tensor, which reflects a later change from strike-slip to compressional stress configuration, the 1 axis remaining unchanged.

Interestingly, tensor 2 from Vein 1 and tensor 1 from Vein 2 are also very close in terms of stress magnitudes with only 0.1 of difference in stress ratio and 9 MPa maximum difference in differential stress magnitude (1-3).

On a methodological point of view, these results show that despite the natural heterogeneity of the calcite filling the veins, CSIT-2 captures the stress regimes consistent with opening of successive regional vein sets which developed during successive substages of the tectonic evolution of the Monte Nero anticline [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF]. The pre-folding N050 compression is related to the regional Apenninic compression which formed the Umbria-Marches domain and started during Tortonian times [START_REF] Barchi | The Umbria-Marche Apennines as a double orogen: Observations and hypotheses[END_REF][START_REF] Billi | Fracture-controlled fluid circulation and dissolutional weathering in sinkhole-prone carbonate rocks from central Italy[END_REF][START_REF] Marshak | Mesoscopic fault array of the northern Umbrian Appennine fold belt, Italy: geometry of conjugate shear by pressure-solution slip[END_REF][START_REF] Storti | The evolution of a model trap in the central apennines , italy : fracture patterns , fault reactivation and development of cataclastic rocks in carbonates at the narni anticline[END_REF][START_REF] Tavani | Stratigraphic versus structural control on the deformation pattern associated with the evolution of the Mt[END_REF]. The subsequent pre-folding N090 to N110 compression is interpreted as a perturbation of the Apenninic compression caused by a N-S striking fault at depth beneath the fold, likely inherited from Tethyan rifting, and which reactivation as a high-angle thrust controlled the structure of the Monte Nero anticline [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF].

Finally, differential stresses obtained from vein 1 (63 MPa ± 3-6% for the first tensor and 35 MPa ± 8-9%) and 2 (38-44 MPa ± 9-10%) are in good agreement with independent differential stress estimates based on tectonic stylolite roughness paleopiezometry at the same locations [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF]. Taking into account an average burial of the Maiolica formation of 2700 m at the time of deformation as derived from the stratigraphic succession, these values are also in good agreement with predicted ranges of differential stress values for compressional strike-slip stress regime (opening of Vein 1) and strike-slip stress regime (opening of Vein 2) at this depth based on differential stress/depth relationship in the crust as reported by [START_REF] Lacombe | Comparison of paleostress magnitudes from calcite twins with contemporary stress magnitudes and frictional sliding criteria in the continental crust: Mechanical implications[END_REF].

As shown in figure 23A, C axes and poles to twinned and untwinned planes for both veins show a non homogeneous spatial distribution. Even though the high level of consistency of the results argues in favor of the robustness of the technique to such potential bias, we also analyzed another sample containing a vein from the same set (N050 trending bed-perpendicular set) than V1 from the former sample that shows a more homogeneous distribution of C axes and poles to twinned and untwinned planes. This sample is located (Fig. 23A) in the western part of the Monte Nero anticline, also in the Massicio formation. The results of the inversion show a strike-slip configuration with a N45 tending σ₁ axis and a perpendicular horizontal σ₃ axis (Fig. 23B) consistent with the opening of the vein (N50-75S in present attitude). The maximum differential stress is 35.6 MPa (± 7%) with a stress ratio of 0.59.

Another tensor has been revealed by the inversion that corresponds to a N118 compression with a differential stress value of 73.9 MPa (± 7%) and a stress ratio of 0.44The 118° compression likely postdates the strike-slip regime with a N45 σ₁ axis which is consistent with vein opening.

In terms of orientations, both stress tensors are consistent with those determined independently from the nearby former sample, which supports the ability of the technique to provide consistent, regionally significant stress tensors. The consistency of the results from both samples further indicates that a slight preferred orientations of the C axes and poles to twinned and untwinned planes exerts a negligible bias on stress tensor determination using CSIT-2.

The 35.6 MPa and 73.9 MPa differential stress values associated with the N45 strike-slip and the N118 compressional stress regimes are different from the values obtained from the former sample for the corresponding tensors (63 MPa and 52/59 MPa, respectively). However, these values are associated with a much greater uncertainty than those derived from the tests on synthetic datasets because the low number of grains precluded any separation of the natural dataset into classes of homogeneous grain size, so that the 9 MPa value of the CRSS adopted for calculation is badly constrained. In spite of this additional effect, the values for the second sample remain of the same order than those derived from the first sample and fall within the range of differential stresses reported by [START_REF] Beaudoin | Fingerprinting stress : stylolite and calcite twinning paleopiezometry reveal the complexity of progressive stress patterns during folding -the case of the Monte Nero anticline in the Apennines, Italy[END_REF] in the Monte Nero anticline.

Conclusions

Inversion of calcite twin data is to date an efficient approach to determine both paleostress orientations and magnitudes in the upper crust. This paper proposes a new technique inspired from the CSIT [START_REF] Etchecopar | Etudes des états de contrainte en tectonique cassante et simulations de déformations plastiques (approche mathématiques[END_REF] but which circumvents its main limitations. The ability of the new technique to detect, to separate and to determine stress tensors from monophase and polyphase twin datasets including measurements errors or various grain sizes is demonstrated by numerous tests on synthetic (numerically generated) twin datasets. In contrast to [START_REF] Gągała | Performance maps -a tool for examination of reliability of procedures for automatic separation of heterogeneous fault / slip and calcite twin data[END_REF], [START_REF] Rez | Peek inside the black box of calcite twinning paleostress analysis[END_REF] and Yamaji (2015), who do not report any application of their techniques to natural data, this paper shows that apart from synthetic samples for which the physical problems related to stress transmission, grain rheological properties (e.g., elasticity) and grain-scale stress heterogeneities are neglected, the technique also yields reliable paleostresses from naturally deformed polyphase samples. It is however strongly advised to handle the results provided by any stress inversion techniques with care, and to

Figure 1 :

 1 Figure 1: Scheme of a twin lamella (pink) in a calcite crystal (yellow). B) Stereographic projection (lower hemisphere, equal area) of the e-twinning system of the calcite crystal. The optical axis C is vertical at the center of the diagram. The poles of the three twin plane sets are [e1], [e2] and [e3].The planes of twinning are the large circles; they contain the twinning direction [ei : rj] ; for each twin plane, the arrow is parallel to the twinning direction ; its head indicates that the upper part of the crystal moves upward, toward the C axis, as a reverse microfault(After Turner et al., 1954).

Figure

  Figure 2: Workflow of CSIT-2 for inversion of calcite twin datasets.

3. 4 Figure 7 :

 47 Figure 7: Evolution of the L²-norm value against variations of stress tensor parameters. A) Illustration of different tensor configurations and resulting L²-norm value. B) Sensitivity of L²-norm value to stress tensor parameters (principal stress orientation, stress ratio and differential stress).

Figure 10 :

 10 Figure 10: Inversion results of numerically-generated monophase datasets M9, M10 and M12. The graphs show A) the results of the L²-norm error and deviation between the applied tensor and the solution tensor, B) the error and deviation on the σ1 orientations, C) the deviation of stress ratio values with the mean value and D) the deviation differential stress values with the mean value. The red line represents the reference (applied) stress ratio and differential stress values.

Figure 21 :

 21 Figure 21: Simplified geological map and stratigraphic log of the Monte Nero anticline (northern Appenines) with location of the sample studied (modified after Beaudoin et al., 2016).

Figure 23 :

 23 Figure 23: Inversion results of naturally deformed polyphase datasets. A) Plots of C axes, and poles of twinned and untwinned planes with distribution of related grain sizes. B) Results of inversion of calcite twins (for more details, see text). The first tensor is oriented with σ₁ : N298-8 and σ₃ : 68-77. The second tensor is oriented with σ₁ : N45-12 and σ₃ : 310-24.

  However, the 10 MPa value of the CRSS was still blindly considered for the inversion process, in order to test the robustness of the technique with respect to significantly varying grain sizes in an aggregate. Figures11A and Bshow the application of inversion process for applied tensors with variable differential stress of 50 MPa and 75 MPa and two perfectly distinct grain sizes. The solution tensor is nearly similar to the applied tensor for 50 MPa, whereas it differs more importantly from the applied tensor with 75 MPa, especially in terms of differential stress. According toYamaji (2015), there is a loss of information on the differential stress value given by the inversion if this value is greater than 50-100 MPa. For these applied tensors, we further tested a constant virtual CRSS not equal to 10 MPa but which instead corresponds to the CRSS applied for each grain size weighted by the number of grains with each size class. The results are worse in terms of differential stresses (43.2 MPa instead of 50 MPa or 41.7 MPa instead of 75 MPa). This confirms that a good estimate of the CRSS value does not depend solely on the average of CRSS values on measured grains but also of the spatial distribution of twinned planes consistent with the applied tensor.

	Figure 11: Inversion results
	of numerically-generated-
	monophase datasets A) M7
	and B) M8 (Table 2). C)
	Histogram showing the
	number of grains of each
	size class used during simu-
	lation.

4.1.2 Monophase twin dataset with heterogeneous grain size

Datasets with virtually heterogeneous grain sizes are simulated using two separate CRSS values. These values are 5 MPa (2/3 of grains) and 15 MPa (1/3 of grains) (configurations M7 and M8 in table

2

).

Acknowledgements:

B. Célérier, K. Amrouch, A. Yamaji and the editor J.-P. Avouac are thanked for their constructive comments that greatly improved the manuscript. We also thank a lot E. Kholer for his help in drafting this manuscript. Authors are most obliged to N. Beaudoin who provides the natural samples from the Monte Nero anticline.

value for solution-tensors found by CSIT-2 are way lower than the one calculated using CSIT.

Especially for the configuration B27 (Table 3) where CSIT finds a stress ratio of 0.06 instead of 0.5 which significantly increases the L²-norm value. The stress axis orientations are retrieved by both CSIT and CSIT-2. The differential stress is also well retrieved by both techniques with less error on CSIT-2 results than using CSIT.

For the second configurations B26 then B28 (Table 3), the results are shown in the table 6. This configuration is different from the first one only about the differential stress value which is higher, 75 MPa for both tensors. As said previously, there is a loss of the information about differential stress between 50-75 MPa of applied differential stresses. Again, the distance between the applied tensors and the solution tensors is lower for the results using CSIT-2 than for CSIT. Especially for the first tensor (B26 in table 3) for which CSIT strongly overestimates the differential stress. The tensor found is clearly intermediate between the first and the second applied tensors.

Table 6: Results of the inversion using CSIT and CSIT-2 of the configuration B26 then B28. The L²norm is the distance between the applied tensor and the solution-tensor found by CSIT or CSIT-2.

For the last configurations B24 then B29 (Table 3), the results are shown in table 7. This configuration is the most difficult to retrieve for CSIT-2 as shown in figure 18, especially for the stress ratio value. In table 7, the L²-norm distance calculated between applied and solution tensors is therefore statistically valid at the regional scale. These vein sets are associated with sets of bedperpendicular stylolites of tectonic origin and are parallel to the teeth they exhibit (peaks oriented N050 and N090-110). This observation supports the fact that both fracture sets developed under subhorizontal compression (1 axis) before folding.

We carried out U-stage measurements of calcite twins within these two veins from 3 mutually perpendicular thin-sections. Calcite grains display thin (0.5 m) and rectilinear e-twins, which cross the crystals over all their length (Fig. 22A). These characteristics indicate that calcite grains probably underwent weak strain, 3-4% at most, under conditions of temperature lower than 150°C-200°C [START_REF] Ferrill | Calcite twin morphology: a low-temperature deformation geothermometer[END_REF]. For each vein, we measure and report the orientations of C axes and of twinned and untwinned planes and we estimate grain sizes using the 3 thin-sections (Fig. 22B) to avoid as much as possible bias related to 2D observations. We could define a single grain size class for Vein 1 and two grain sizes for Vein 2. One should keep in mind that because these classes must contain a sufficient number of grains for the inversion to be statistically valid (at least 30 grains), they must not be too narrow. Inversion is carried out for each grain size class independently. For Vein 1 the CRSS considered is of 9 MPa. For Vein 2, the CRSS value is of 10 MPa for the lower grain size (200-400 µm) and 7.5 MPa for the higher grain size (400-600 µm). Note that although the grain size class of Vein 1 is similar to the smaller grain size class of Vein 2, the adopted CRSS value are slightly different (9 and 10 MPa, respectively). This is because the CRSS is also dependent on internal twinning strain, which is roughly estimated using twin density to be higher in Vein 2 than in Vein 1. This way to take into account strain hardening is similar to [START_REF] Rocher | Intraplate paleostresses reconstructed with calcite twinning and faulting: improved method and application to the eastern Paris Basin (Lorraine, France)[END_REF] and [START_REF] Amrouch | Stress and strain patterns, kinematics and deformation mechanisms in a basement-cored anticline: Sheep Mountain Anticline[END_REF].

Taking into account the crosscutting relationships between the veins (Fig. 22A), twinned calcite from Vein 1 expectedly recorded (at least) the state of stress responsible for its opening, but also the later state of stress responsible for the formation of Vein 2. In contrast, calcite from Vein 2 should not have recorded the state of stress responsible for the earlier opening of Vein 1 but should have expectedly recorded the state of stress responsible for its own formation.

have a critical look at the results, making the best use of available geological information (e.g., relative chronology) to guide the inversion procedure and separation of superimposed stress tensors.

The main technical limitations of paleostress reconstructions from calcite twins to date, whatever the inversion technique is used, are related to (1) potential errors in optical measurements;

(2) the influence of grain size and grain size distribution in the deformed aggregates, and (3) our insufficient knowledge of the variations of the CRSS value for twinning with grain size. While the first is likely to be overcome by twin data acquisition using EBSD [START_REF] Parlangeau | Inversion of calcite twin data for paleostress (2): EBSD as a tool for data measurements[END_REF] and the second is possibly circumvented by defining several grain size classes and treating them separately, understanding the influence of the grain size distribution and characterizing the way the CRSS varies with changing grain size are challenges for forthcoming studies.