Fuzzy uncertainty representations of co-seismic displacement measurements issued from SAR imagery

Yajing Yan, Gilles Mauris, Member, IEEE, Emmanuel Trouvé, Senior Member, IEEE, Virginie Pinel

Abstract-An emerging way to reduce the geodetic parameter uncertainty is to combine the large numbers of data provided by satellite synthetic aperture radar (SAR) images. However, the measurements by radar imagery are subject to both random and systematic uncertainties. Thus, mathematical theories which are adequate for each type of uncertainty representation and handling have to be selected. Probability theory is known as the adequate theory for uncertainties corresponding to random variables, but questionable for systematic uncertainties, arising from information incompleteness. Fuzzy theory being a generalization of interval mathematics, is more adapted to such uncertainty. Moreover it provides a bridge with probability theory by its ability to represent a family of probability distributions. Therefore, we consider here the conventional probability and the fuzzy approaches for handling the random and systematic uncertainties of D-InSAR and SAR amplitude image correlation measurements. The applications are performed on the measurement of the displacement field due to the 2005 Kashmir earthquake. The fuzzy approach, being free from distribution and independence hypotheses, gives the most pessimistic uncertainty assessment, while the conventional probability approach gives the most optimistic uncertainty assessment. As confirmed by the Monte-Carlo simulation applied to an Earth deformation model, the actual uncertainty should be situated between the fuzzy and conventional uncertainties.

Index Terms-measurement uncertainty, fuzzy/possibility theory, remote sensing, SAR images, ground displacement, coseismic observations

I. INTRODUCTION

Geodetic data, such as satellite images (radar and optic), are important remote sensing sources of information for ground displacement measurement with great accuracy over large area. So far, with the increasing number of operational sensors, large volumes of SAR images acquired in different modes, ascending and descending passes at various incident angles and frequencies, are available. Moreover, the launching in the coming years of the future satellite generation Sentinel will provide a large number of free SAR data [START_REF] Attema | Mission Requirements Document for the European Radar observatory Sentinel-1[END_REF]. Consequently, using large number of geodetic measurements in order to accurately determine the displacement field is becoming more and more frequent in geophysics, especially to better constraint the geophysical modeling [START_REF] Pedersen | Fault slip distribution of two June 2000 Mw6.5 earthquake in South Iceland estimated from joint inversion of InSAR and GPS measurements[END_REF] [3] [4] [START_REF] Grandin | September 2005 Manda Hararo-Dabbahu rifting event, Afar (Ethiopia): Constraints provided by geodetic data[END_REF]. In this context, one Laboratoire d'Informatique, Systèmes, Traitement de l'Information et de la Connaissance, Université de Savoie -Polytech Annecy-Chambéry, BP80439, F74944, Annecy-le-Vieux Cedex, France. Email: yajing.yan,gilles.mauris,emmanuel.trouve@univ-savoie.fr ISTerre, CNRS, IRD R219, Université de Savoie, Campus Scientifique, 73376, Le Bourget du Lac Cedex, France. Email: virginie.pinel@univsavoie.fr important purpose of geodetic data processing is to reduce parameter uncertainty by an adequate combination of all the available measurements.

Once gross errors have been eliminated, measurement uncertainties come from random and systematic effects. Concerning the measurements by radar (SAR) imagery, the uncertainties arise from noise sources of radar instrument, on the path of radar wave propagation, at the reflecting surface, as well as uncertainty sources introduced by data processing [START_REF] Sudhaus | Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland[END_REF]. On one hand, random uncertainty exists due to decorrelation noise, since there are usually some backscattering property changes on the ground between two subsequent SAR acquisitions. On the other hand, systematic uncertainties can be induced by atmospheric disturbances depending on the state of atmosphere and the ground surface at the time of the two SAR acquisitions. Also, it can result from the imprecision of orbit auxiliary information, Digital Elevation Model (DEM) errors, as well as from the imperfect corrections during data processing, which deviate the data by a constant or a ramp from the true value.

To model such systematic uncertainties coming from limited knowledge (called also epistemic uncertainties in many cases), probability theory is questionable, and thus fuzzy/possibility theory has been proposed in ( [START_REF] Mauris | Fuzzy modeling of measurement data acquired from physical sensors[END_REF], [START_REF] Mauris | Expression of measurement uncertainty in a very limited knowledge context: a possibility theory-based approach[END_REF]) and further developed by a few authors in a general measurement context [START_REF] Urbanski | Fuzzy approach to the theory of measurement inexactness[END_REF] [9] [10] [START_REF] Mencattini | How to process the random part of RFVs: comparison of available methods and new proposals[END_REF], and also by a few authors in geosciences [START_REF] Verhoest | A possibilistic approach to soil moisture retrieval from ERS synthetic aperture radar backscattering under soil roughness uncertainty[END_REF] [13] [14] [START_REF] Yan | Fuzzy vs probability uncertainty analysis of seismic displacement measurements issued from d-insar and sar image correlation measurements: Application to the kashmir earthquake[END_REF]. Thus, it is worthwhile to study the consideration of the most suitable uncertainty theories according to the application specificities. This paper is a first contribution to such issue for the measurement of displacement field by a joint inversion of D-InSAR and SAR amplitude image correlation measurements by a least squares adjustment.

This paper is organized as follows. In section II, the conventional probability approach and the fuzzy approach for uncertainty representation and propagation in the linear case are detailed. The available data in the considered application and their associated uncertainties are described in section III. Then, the two approaches are applied to the measurement of the Earth's surface 3D displacement field due to the 2005 Kashmir earthquake (magnitude of 7.6). The behaviours of each uncertainty approach are highlighted through intercomparisons of results. In section IV, a model based reference of the 3D displacement field is established by using a mechanical Earth deformation model in order to evaluate the results obtained previously. The interest of the fuzzy approach is highlighted with respect to the results provided by the physical model. Finally, some conclusions and perspectives are drawn.

II. UNCERTAINTY REPRESENTATION APPROACHES

In geodetic practice, there are always multiples sources of uncertainty in the considered measurements, which leads to complex characteristics for the associated uncertainties. Probability theory is the adequate theory for uncertainties corresponding to random variables, the latter being described by one probability distribution (often a Gaussian one) or more simply by the first two moments, i.e. the mean and the variance. Systematic uncertainties, arising from information incompleteness, are often described by an interval and thus cannot be associated to one single probability distribution. Fuzzy theory can be considered as a generalization of interval mathematics and provides a bridge with probability theory by its ability to represent a family of probability distributions [START_REF] Dubois | Probability-possibility transformations, triangular fuzzy sets and probabilistic inequalities[END_REF]. Therefore, hereafter, we consider the conventional probability and fuzzy approaches for handling the random and systematic uncertainties encountered in SAR displacement measurements.

A. Conventional probability approach

The standard reference in uncertainty modeling is the Guide to the Expression of Uncertainty in Measurement (GUM) edited by an international consortium of legal and professional organizations [START_REF] Bipm | IUPAC, IUPAP, and OIML, Evaluation of measurement data -Guide to the expression of uncertainty in measurement[END_REF]. The GUM groups the occurring uncertain quantities into Type A and Type B. Uncertainties of Type A are determined with the classical statistical methods, while Type B is determined by other ways, e.g. experience and knowledge about an instrument. Both types of uncertainty can have random and systematic components. In fact, the GUM proposes to treat both uncertainties (random and systematic) in a stochastic framework and considers variances to describe them and processes them with the law of propagation of variance, generally assuming independence. Applying this approach to linear inversion by a least squares adjustment, the uncertainties are propagated as follows. Let us consider U a set of m quantities u k (x 1 , x 2 , . . . , x n ) (e.g. 3D (E, N, Up) Earth surface displacement) which are linear combinations of the vector of n measured variables R = (r 1 , r 2 , . . . , r n ) (e.g. projections of the 3D displacement) with coefficients vector A = (a 1,k , a 2,k , . . . , a n,k ), (k = 1, . . . , m). Thus:

u k = n i=1 a i,k r i : U = A T R (1) If the variance-covariance matrix of R is denoted by Σ R , Σ R =      σ 2 1 σ 12 σ 13 • • • σ 21 σ 2 2 σ 23 • • • σ 31 σ 32 σ 2 3 • • • . . . . . . . . . . . .      (2)
then, the variance-covariance matrix Σ U of U is given by

Σ U i,j = n k=1 n l=1 a i,k Σ Rk,l a l,j : Σ U = A T Σ R A (3)
When R and U are related by a n × m matrix P which represents the geometrical relationship between the measured variables and the displacement to be determined, the inversion by the generalized least squares method [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] [19] leads to:

A T = (P T Σ -1 R P ) -1 P T Σ -1 R (4)
thus the nominal value U is given by

U = (P T Σ -1 R P ) -1 P T Σ -1 R R (5)
and the uncertainty Σ U is given by

Σ U = (P T Σ -1 R P ) -1 (6)
with Σ R the variance-covariance matrix of R.

To determine confidence intervals for the parameters, the GUM suggests to use a Gaussian distribution (justified by the central limit theorem), and for other distributions to apply Monte Carlo simulations. This well known approach is fully justified in cases of a lot of data having independent random uncertainties but questionable for systematic uncertainties often dependent and far from a Gaussian representation. Therefore, it generally leads to an over-optimistic assessment of the uncertainties.

B. Fuzzy/possibility approach

The possibility theory, first introduced by L. Zadeh in 1978 [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF], is associated with the theory of fuzzy sets by the semantics of uncertainty that it gives the membership function. A possibility distribution π is a mapping from a set to the unit interval such that π(x) = 1 for some x belonging to the set of reals (Fig. 1). A possibility distribution π 1 is called more specific (i.e. thinner in a broad sense) than π 2 as soon as ∀x ∈ ℜ, π 1 (x) ≤ π 2 (x) (fuzzy set inclusion). The more specific π is, the more informative it is. If for one x, π(x) = 1 and for all y = x, π(y) = 0, then π is totally specific (fully precise and certain knowledge). If for all x, π(x) = 1, then π is totally non specific (complete ignorance) (Fig. 1). In fact, a numerical degree of possibility can be viewed as an upper bound to a probability degree [START_REF] Dubois | Probability-possibility transformations, triangular fuzzy sets and probabilistic inequalities[END_REF]. Namely, with every possibility distribution π one can associate a non-empty family of probability measures dominated by the possibility measure: P(π) = {P, ∀A, P (A) ≤ Π(A)}. This provides a bridge between probability and possibility. There is also a bridge with interval calculus. Indeed, a unimodal numerical possibility distribution may also be viewed as a nested set of coverage intervals, which are the α-cuts of π: [x α , x α ] = {x, π(x) ≥ α}. Obviously, the confidence intervals built around the same point x 0 are nested. It has been proven in [START_REF] Mauris | A fuzzy approach for the expression of uncertainty in measurement[END_REF] that stacking the coverage intervals of a probability distribution F on top of one another leads to a possibility distribution (denoted π x0 having x 0 as modal value). In fact, in this way, the αcuts of π x0 are identified with the confidence interval I ⋆ β of probability level β = 1α around the nominal value x 0 . In this way a probability distribution can be represented by an equivalent possibility distribution. Moreover, a possibility distribution can be used to represent a family of probability distributions by taking the largest 1α confidence intervals obtained from each probability distribution of the family. This is useful to represent uncertainty when only partial probability knowledge is available. For example if the measured variable r is known to be bounded and unimodal (with mode r 0 ) and with σ r as standard deviation, then the maximum specific possibility distribution is a triangular possibility distribution with the mode r 0 as vertex and with [r 0 -√ 3σ r , r 0 + √ 3σ r ] as support [START_REF] Mauris | Expression of measurement uncertainty in a very limited knowledge context: a possibility theory-based approach[END_REF]. To represent the uncertainty in the fuzzy approach by a single parameter (in a way similar to the variance in the probability approach), the full width at half maximum (denoted h 0.5 ) of the possibility distribution which corresponds to a α-cut level of 0.5, is commonly used. In summary, a possibility distribution can model both random and systematic uncertainties in a unified modeling. Afterwards, the possibility distributions are propagated in the least squares adjustment using fuzzy arithmetic based on Zadehs extension principle [START_REF] Zadeh | Fuzzy sets as a basis for a theory of possibility[END_REF]:

U = (P T Σ -1 R P ) -1 P T Σ -1 R ⊗ R ( 7
)
where ⊗ is the fuzzy multiplication matrix operation in which the conventional sum and scalar product are replaced by the corresponding fuzzy operations [START_REF] Klir | Fuzzy sets and fuzzy logic: Theory and application[END_REF]. In this principle, the variables are considered as non interactive variables; this corresponds somehow to consider a total dependence between variables. Consequently, uncertainty propagation by the fuzzy approach leads to an over-pessimistic assessment of the uncertainties.

C. Displacement uncertainty analysis by the two approaches

In this article, to combine the available measurements, we apply the Generalized Least Squares (GLS) method [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF] [19] which gives the nominal displacement value and its corresponding variance from the variances of the measurements provided by the SAR image processing. For the fuzzy approach, we use the same GLS matrix A T to obtain the forward model. But we build the possibility distributions of the displacement from the value r 0 and their associated uncertainty σ r0 , considering that σ r0 contains both random and systematic components. Moreover, the measurements are considered as bounded (this is the case in the considered context), thus we represent them by a symmetric triangular fuzzy distribution with support [r 0 -√ 3σ r0 , r 0 + √ 3σ r0 ]. Therefore, the full width at half maximum is √ 3σ r0 . Let us remark that with a Gaussian assumption of standard deviation σ r0 , the value corresponding to the 0.5 α-cut the width of the Gauss equivalent possibility distribution is equal to 1.35 σ r0 , i.e. the width of the 50% confidence interval of a Gaussian distribution with σ r0 as standard deviation.

III. DISPLACEMENT MEASUREMENTS DUE TO THE 2005 KASHMIR EARTHQUAKE

A. Description of the available data

The sub-pixel image correlation and the differential interferometry (D-InSAR) are two conventional techniques used to extract displacement measurements from SAR data. The subpixel image correlation computes the offsets in range (line of sight) and azimuth (along the trajectory of satellite) directions on amplitude images, with a sub-pixel accuracy. It is widely used to measure the displacement of great magnitude [START_REF] Pathier | Displacement field and slip distribution of the 2005 Kashmir earthquake from SAR imagery[END_REF] [24] [START_REF] Wright | Magma-maintained rift segmentation at continental rupture in the 2005 Afar dyking episode[END_REF]. The D-InSAR, using the phase information, provides the displacement in range direction with an uncertainty in the order of centimeter, even millimeter [26] [19]. This technique is usually applied to measure the displacement of small magnitude. In case of a strong earthquake induced by a rupture of a fault, in near field of the fault, the measurements from sub-pixel image correlation can provide reliable displacement information. While in far field of the fault, the measurements from D-InSAR are taken as accurate sources [START_REF] Pathier | Displacement field and slip distribution of the 2005 Kashmir earthquake from SAR imagery[END_REF].

In this article, a series of co-seismic ENVISAT images from October 2004 to June 2006 are used to map the deformation due to the 2005 Kashmir earthquake. 22 measurements from sub-pixel image correlation and 5 measurements from D-InSAR, are available respectively. In near field of the fault, because of coherence loss, phase information cannot be extracted by D-InSAR, thus there is no D-InSAR measurements available in this area. These measurements can be classified in four families according to their acquisition geometry: ascending range (Asc. Rg), ascending azimuth (Asc. Az), descending range (Des. Rg) and descending azimuth (Des. Az). In each family, in first approximation, all the measurements are considered as corresponding to the same displacement (in the same direction) because the incident angle is the same for all the measurements.

For measurements from sub-pixel image correlation, the uncertainty parameter is the so called "pseudo-variance" provided by the ROIPAC software [START_REF] Rosen | Updated Repeat Orbit Interferometry Package released[END_REF]. In fact, this pseudovariance includes both random and systematic components and thus it is not a conventional statistical variance. The random part comes from the decorrelation noise present in the data, while the systematic part comes from the default of correlation method and the DEM error estimation. However, the possible systematic uncertainty due to the imperfect data processing is not included in the pseudo-variance.

For measurements from D-InSAR, the uncertainty parameter corresponds to the variance of the phase value estimated from the coherence [START_REF] Trouvé | Improving Phase Unwrapping Techniques by the Use of Local Frequency Estimates[END_REF]. It characterizes only random variations in the phase value. However, systematic uncertainty due to phase unwrapping errors, atmospheric impact, etc., is probably present in the measurement, but has not yet been taken into account in our approach. 2 shows two examples of the measurements issued from both sub-pixel image correlation and D-InSAR. The two profiles are issued from displacement images. The first profile (Fig. 2 (a)) is located in an area far from the fault and the second profile (Fig. 2 (b)) passes across the fault. In theory, the same displacement values should be found by these two techniques where both measurements are available, as they measure exactly the same quantity of displacement. However, with the presence of uncertainty in both measurements, a discrepancy of displacement values is observed on both profiles. On one hand, a more or less significant fluctuation of displacement value is observed in the sub-pixel image correlation measurements, which complies with the presence of random uncertainty. Near the fault, the fluctuation is small. While in the area far from the fault where the displacement magnitude is small, the fluctuation becomes significant. On the other hand, there is a small shift between the displace-ment values estimated by both techniques. Consequently, it is probable that systematic uncertainty is also present in one or the other measurement. However, it seems that the random uncertainty is more important than the systematic uncertainty in the considered measurements.

B. 3D Displacement field retrieval by linear inversion

The different measurements from sub-pixel image correlation and D-InSAR are different projections of the 3D displacement at the Earth's surface (E, N, Up) in the SAR image range and azimuth directions. Consequently, the 3D displacement field can be constructed from at least 3 different projections by a linear inversion. In this case, R corresponds to different measurements from sub-pixel image correlation and D-InSAR. P corresponds to the projection vectors matrix. U denotes the 3D displacement with 3 components E, N, Up. To solve this linear inverse problem, the GLS method is used. In order to highlight the behaviours of conventional and fuzzy approaches, three levels of comparisons are considered: between displacement values, between uncertainty parameter values and between possibility distributions. Moreover, the effect of uncertainty reduction due to adding D-InSAR measurements is analyzed.

Regarding the nominal displacement value, the result obtained with D-InSAR measurements is shown in Fig. 3 (a). In fact, the results in cases with and without D-InSAR measurements are globally consistent, with an average difference in the order of millimeter. Regarding the uncertainty parameter values, σ and h 0.5 , 4 cases are considered: with and without D-InSAR measurements in both conventional and fuzzy approaches. The evolution of these values varies from one case to another (Fig. 3 (b) and(c)). Uncertainty parameter values in the fuzzy approach are always larger than those in the conventional approach. Adding measurements from D-InSAR reduces the uncertainties in both conventional and fuzzy approaches. In order to understand the spatial evolution of uncertainty in both approaches, the ratio ( 1.35σ h0.5 ) of the conventional uncertainty (Fig. 3 (b)) on the fuzzy uncertainty (Fig. 3 (c)) is shown in Fig. 4. A geographic effect is observed and it corresponds to the distribution of the number of available measurements. In the darker area, the difference between the two uncertainty parameters is large, because in this area, there are more measurements available. In the conventional approach, the output uncertainty is reduced. In the fuzzy approach, on the contrary, the output uncertainty remains constant or increases slightly. Consequently, the difference between the two uncertainty parameters increases in this area.

Fig. 4 provides a summarized vision of the relationship between the two uncertainty parameters. For each point, the comparison is also perfomed between distributions: possibility distribution corresponding to the fuzzy uncertainty and equivalent possibility Gaussian distribution corresponding to the conventional uncertainty. Fig. 5 shows un example of possibility distributions for the point B. As expected, the conventional possibility distribution is more specific than the fuzzy based one.

IV. COMPARISON WITH A MODEL BASED 3D DISPLACEMENT

Since there is no ground truth available to validate the results obtained in the previous section, we compare them with a model based reference of the 3D displacement constructed by means of a well known mechanical deformation model. Blue: fuzzy approach. Red: conventional approach. The black line represents the width of the 0.5 α-cut, the uncertainty parameter used in the fuzzy approach.

Indeed, this model is widely used since many years by geophysicists [START_REF] Clarke | Geodetic investigation of the 13 May 1993 Kozani -Grevena (Greece) earthquake[END_REF] [30] [31] [32] [23] [START_REF] Atzori | Finite fault inversion of DInSAR coseismic displacement of the 2009 L'Aquila earthquake (central Italy)[END_REF]. Thus it is worthwhile to confront the corresponding results with the ones obtained with our approaches.

A. Construction of a model based reference

For the model based reference, the 3D displacement is deduced by 2 steps. First, the fault geometry parameters are optimised by nonlinear inversion of a mechanical deformation model describing the surface displacement field induced by a fault rupture in depth, by looking for a global minimum misfit with the surface measurements from the same input data sets as used in the previous section. At this step, the input data sets are subsampled by quadtree algorithm, thus the number of points in each measurement is greatly reduced in order not to burden the inversion. As a result, each point used in the inversion represents the deformation in an area more or less large. Second, the 3D displacement is derived from the fault geometry parameters obtained at the first step through the forward deformation model [START_REF] Okada | Surface deformation due to shear and tensile faults in a halfspace[END_REF]. In this method, all the points used in the inversion contribute in a joint way and the fault geometry parameters correspond to a global deformation model, which is quite different from the method used in the previous section. Therefore, the displacement field obtained here is very smooth, which reveals the global deformation behaviour. By contrast, in the previous section, the neighboring points are independent, the retrieved 3D displacement reveals essentially the local deformation behaviour and is sensitive to noise.

Because of the complexity of the nonlinear inversion of the deformation model, we cannot propagate uncertainties associated with input measurements through the deformation model directly. In order to estimate the uncertainty associated with the derived 3D displacement, synthetic noise is simulated and added to the input measurements. Measurements from subpixel image correlation are mainly influenced by white gaussian noise, while measurements from D-InSAR are essentially influenced by spatially correlated noise due to atmospheric impact [START_REF] Sudhaus | Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland[END_REF]. The characteristics of these two types of noise are estimated from preseismic data sets using both sub-pixel image correlation and D-InSAR on pairs of images without deformation. Then, 1000 simulations of white gaussian noise or spatially correlated noise with the same characteristics are carried out for each measurement according to the data type. 1000 sets of fault geometry parameters are obtained, from which 1000 values of the 3D displacement are retrieved for each pixel. Then the mean and the standard deviation are calculated, and the corresponding triangular possibility distribution is built for each pixel according to the way described in Section II-B.

The fuzzy uncertainty is used in the model based reference in order to take into account the lack of complete probability knowledge in modeling. Although the model used in this section is not perfect, it is commonly used and validated by most of the geophysicists since more than 20 years. Therefore, the results obtained in the previous section are worthwhile to be compared to the results obtained by this method.

B. Comparison with the model based reference

In order to facilitate the comparison, a profile, crossing the fault, as shown in Fig. 3 (a), is performed and the comparison is realised along this profile. The results are plotted in Fig. 6. The uncertainty (1.35σ for the conventional approach, h 0.5 for the fuzzy approach and the model based reference) is here represented by uncertainty bar and the jumping of displacement value (at about 23.5 km from the reference point A) indicates the position of the fault rupture.

According to Fig. 6, there is a good agreement between the results obtaind in the previous section and the model based reference. For most of the pixels, there is a good superposition, especially in the field far from the fault. Near the fault, for the model based reference, the maximal displacement is under estimated, which results in the shift of nominal displacement value with respect to that retrieved in the previous section. However, taking into account the uncertainty associated with the nominal displacement value, the difference is small. Thus, the good agreement, especially with the fuzzy approach, allows validating the results obtained in the previous section. It is highlighted that the uncertainty in the conventional approach is optimistic, while the uncertainty in the fuzzy approach is pessimistic but closer to the model based reference in this comparison.

Moreover, in order to make use of the richness of the possibility distribution, possibility distributions obtained in the previous section and with the model based reference are compared. Here, besides point B, we consider another point C situated on the fault (as shown in Fig. 3 (a)). The possibility distributions are shown in Fig. 7.

In Fig. 7 (a), there is a good intersection between the possibility distributions obtained in the previous section and with the model based reference. In particular, the intersection with the possibility distribution of the fuzzy approach (blue) is more significative than with the equivalent Gaussian distribution (red). Therefore, the possibility distribution of the fuzzy approach is more relevant than the equivalent Gaussian distribution. Fig. 7 (b) confirms this conclusion. In this case, the difference of results with respect to the model based reference is larger, thus there is no significative intersection with the equivalent Gaussian distribution, but a small one with the possibility distribution of the fuzzy approach.

V. CONCLUSION

In this article, two approaches based respectively on the probability and fuzzy theories have been applied to represent uncertainty in measurement of the 3D displacement field due to the 2005 Kashmir earthquake. The 3D displacement field is estimated by a least squares based linear inversion using the measurements from sub-pixel image correlation and D-InSAR. The uncertainty is propagated in the conventional approach and in the proposed fuzzy approach respectively. The results are compared with a model based reference issued from a nonlinear inversion of a well known mechanical Earth deformation model.

According to the foundations of probability and fuzzy theories, the uncertainty in the considered probability approach is under estimated, while the uncertainty in the proposed fuzzy approach is over estimated. In the context of our measurements, the actual uncertainty should be situated between these two uncertainties. The more the assumptions of Gaussian distribution and independence on uncertainties associated to input measurements, the closer to the conventional uncertainty it is. Otherwise, it is closer to the fuzzy uncertainty. Although the proposed fuzzy approach is pessimistic in uncertainty assessment, it is more robust against the assumptions of Gaussian distribution and independence on uncertainties made in the conventional approach. In fact, the fuzzy and conventional approaches provide respectively the most pessimistic and the most optimistic uncertainty assessment from almost similar simple analytical computations. Monte-Carlo simulations applied to the Earth deformation model provide a more relevant assessment, but with much more computational complexity and they are extremely time consuming.

In order to improve the uncertainty analysis, the part of systematic uncertainties not yet taken into account in the SAR measurements has to be identified and represented in an appropriate way. Besides this subject, the fusion of the results of the direct SAR measurements with the ones provided by the mechanical deformation model, is also a part of ongoing works.

Fig. 1 :

 1 Fig. 1: Example of fuzzy possibility distributions.

Fig. 2 :

 2 Fig. 2: Examples of displacement measurements. Red: displacements estimated by sub-pixel image correlation. Green: displacements estimated by D-InSAR

Fig. 3 :Fig. 4 :

 34 Fig. 3: (a) 3D displacement values (b) conventional uncertainty parameter values (σ) (c) fuzzy uncertainty parameter values (h 0.5 ) estimated in the linear case. The profile illustrated in (a) Up component, is used to compare the displacement values and uncertainties obtained in both linear and nonlinear case (Fig.6). Point B and C are used for possibility distribution analysis (Fig.5, Fig.7).

Fig. 5 :

 5 Fig. 5: Comparison of possibility distributions on point B.Blue: fuzzy approach. Red: conventional approach. The black line represents the width of the 0.5 α-cut, the uncertainty parameter used in the fuzzy approach.

Fig. 6 :

 6 Fig. 6: Profile of displacements and uncertainties obtained by linear inversion and with the model based reference (a) East component (b) North component (c) Up component.

(Fig. 7 :

 7 Fig. 7: Comparison of possibility distributions. blue: fuzzy approach. red: conventional approach. black: reference.

ACKNOWLEDGMENT

This work was supported by the EFIDIR project (http://www.efidir.fr) (ANR-07-MDCO-004) and the PAKSIS program CATTELL, granted by the French National Agency (ANR). The authors wish to thank the ANR for their support. The ENVISAT data sets are available thanks to ESA project category 1, No. 3802.


David Amitrano

Variability in the power-law distributions of rupture events How and why does b-value change

Keywords: -key Failure materials, 81.40.Np --key dynamic critical behavior, 64.60.Ht

Rupture events, as the propagation of cracks or the sliding along faults, associated with the deformation of brittle materials are observed to obey power-law distributions. This is verified at scales ranging from laboratory samples to the Earth's crust, for various materials and under various loading modes. Besides the claim that this is a universal characteristic of the deformation of heterogeneous media, spatial and temporal variations are observed in the exponent and tail-shape. These have considerable implications for the ability and the reliability of forecasting large events from smaller ones. There is a growing interest in identifying the factors responsible for these variations. In this present work, we first present observations at various scales (laboratory tests, field experiments, landslides, mining induced seismicity, crustal Earthquakes) showing that substantial variations exist in both the slope and the tailshape of the rupture event size distribution. This review allows us to identify potential explanations for these variations (incorrect statistical methods, heterogeneity, stress, brittle/ductile transition, finite size effects, proximity to the failure). A possible link with the critical point theory is also drawn showing that it is able to explain a part of the observed variations considering the distance to the critical point. Using numerical simulations of progressive failure we investigate the role of mechanical properties on the powerlaw distributions. The results of simulations agree with the critical point theory for various macroscopic behaviors ranging from ductility to brittleness providing a unified framework for the understanding of power-law variability observed in rupture phenomena.

Introduction

The deformation of brittle materials is associated with damage (i.e. micro-cracks nucleation and propagation) and friction (i.e. sliding on preexisting cracks or geological faults). These processes occur by the accumulation of discrete events of size spreading over several orders of magnitude. This is observed at scales varying from that of laboratory samples to that of the Earth's crust, including underground excavations, quarries, cliffs, landslides and volcanoes. Scaling invariance is observed for temporal (e.g. the distribution of waiting time), spatial (e.g. fractal structure of the rupture) and size distributions, displaying the complexity of the rupture dynamics. Here we focus on the ubiquitous observation that the size distribution of the rupture events obeys a power-law characterized by its exponent, called b-value in the Earth sciences field. Despite the fact that the mean value of this exponent is robustly observed to be near 1 for Earthquakes, substantial variations exist in both the slope (e.g. [START_REF] Rosso | Depinning of elastic manifolds[END_REF]) and the shape of the distribution for largest events, i.e. the tail-shape (e.g.

Present address: ISTerre, University J. Fourier, Grenoble, BP 53, 38041, Grenoble Cedex 9, France [START_REF] Main | Apparent breaks in scaling in the earthquake cumulative Frequency-Magnitude distribution: Fact or artifact?[END_REF]) of the rupture event size distribution. Since the early works of Mogi [START_REF] Mogi | Magnitude frequency relations for elastic shocks accompanying fractures of various materials and some related problems in earthquakes[END_REF] and Scholz [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF] on brittle rocks, there is a constant effort to understand the origins of these fluctuations and whether material properties and mechanical loading conditions are able to explain them.

The dynamics of fracturation and damage during mechanical loading usually displays a power-law distribution of rupture event sizes.

N (> s) = c.s -β ( 1
)
where s is a measure of the rupture event size, N (> s) is the number of events with size larger than s, c and β are constants. N (> s) is the empirical cumulative distribution function cdf . The frequency size relationship can also be examined using the probability density function, pdf , which is obtained by derivating the cumulative distribution function. In the case of power-law, the exponent of the pdf equals -β -1.

p(s) = c .s -β-1 (2)
In numerical simulations, the rupture event size can be estimated by the number of broken bonds during an avalanche for fracture models, or by the area of a slip zone in friction models. In laboratory experiments, when using acoustic emission (AE) monitoring for observing the damage dynamics, s can be estimated by the signal maximum amplitude A or energy E. The absolute value of the β exponent depends on the metrics used for characterizing the event size s. In AE monitoring, theoretical considerations, well verified empirically allows us to consider that E scales with A 2 [START_REF] Lockner | The role of acoustic emission in the study of rock fracture[END_REF], so that both N (> A) and N (> E) obey power-law,

N (> A) = c A .A -βA
(3)

N (> E) = c E .E -βE (4)
but the exponents will be different in the two cases, with the scaling relationship β A = 2β E . In a log-log representation, such power-law distributions appear linear and exponent β is given by the slope of the line. logN (> s) = logcβ.logs [START_REF]Failure by damage accumulation in rocks[END_REF] This distribution exhibits remarkable similarity with the Gutenberg-Richter empirical relationship observed for earthquakes [START_REF] Gutenberg | Seismicity of the earth and associated phenomena[END_REF].

logN (> M) = a -bM (6)
Where N (> M ) is the number of earthquakes with a magnitude larger than M . The magnitude M is theoretically and empirically related to the energy E released by the seismic rupture [START_REF] Hanks | A moment magnitude scale[END_REF][START_REF] Scholz | The mechanics of earthquakes and faulting[END_REF] M = 2 3 logE -11.8 [START_REF] Amitrano | Seismic precursory patterns before a cliff collapse and criticalpoint phenomena[END_REF] where E is in ergs. So the values of b, β A and β E are linearly related b = 3 2 β E = 3 4 β A . Figure 1 shows an example of Magnitude and Energy distributions observed for California Earthquakes with the corresponding exponents. The former relation between the exponents is empirically well verified.

As power-laws indicate scale invariance and because of the similarities in the underlying physics, AE of rocks observed in the laboratory has been considered as a smallscale model of the seismicity of the Earth's crust [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF]. Observations of both earthquakes and AE show variations of the b-value in time and space domains. The purpose of this work is to identify possible origins of these variations and to provide a physical model able to explain them. The paper is organized as follows: i) A review of laboratory and field observations of b-value variations and outcomes of numerical models able to produce power-law distributions of rupture events, providing a short list of potential explanations of the b-value variations. ii) Detailed analysis of each potential explanation and of their effect on b-value. iii) Numerical modeling showing the effect of the mechanical properties on both the power-law exponent and the tail-shape. Analysis of the agreement with critical point theory.

2 Observed variations in the power-law distributions of rupture events

Experimental and Field observations

Mogi [START_REF] Mogi | Magnitude frequency relations for elastic shocks accompanying fractures of various materials and some related problems in earthquakes[END_REF], suggested that the b-value depends on the material heterogeneity, a low heterogeneity leading to a low exponent. Scholz [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF] observed that the b-value decreased before the stress peak of rock specimens compressed in the laboratory and argued for a negative correlation between b-value and stress. The relationship between the b-value and the fractal dimension D of AE source locations has also been investigated in the laboratory [START_REF] Lockner | Precursory AE patterns leading to rock fracture[END_REF] and showed a decrease of b-value in conjunction with the strain localization, i.e. to a decrease of D-value. These patterns have been observed essentially for brittle rocks (e.g. [START_REF] Lei | Detailed analysis of acoustic emission activity during catastrophic fracture of faults in rock[END_REF][START_REF] Hirata | Fractal structure of spacial distribution of microfracturing in rock[END_REF]) and for other heterogeneous materials such as chipboard, fiberglass [START_REF] Guarino | An experimental test of the critical behaviour of fracture precursors[END_REF] or paper [START_REF] Rostil | Statistics of acoustic emission in paper fracture: Precursors and criticality[END_REF] providing arguments for the existence of a relationship between the b-value decrease to the progressive spatial clustering of the damage.

Temporal fluctuations of the b-value have also been observed for earthquakes leading some authors to consider the b-value decrease as a precursor to major earthquakes being considered as macro-failure events [START_REF] Ch Lin | Temporal b-Value variations throughout a seismic faulting process: The 2008 Taoyuan earthquake in Taiwan[END_REF][START_REF] Smith | The b-value as an earthquake precursor[END_REF][START_REF] Tsukakoshi | Decreased b-value prior to the m 6.2 northern miyagi, japan, earthquake of 26 july 2003[END_REF][START_REF] Urbancic | Space-time correlations of b-values with stress release[END_REF]. Observations of the seismicity emitted by a rock cliff before a collapse revealed that such a decrease of the b-value was associated with a power-law acceleration of the event rate and of the released energy rate [START_REF] Amitrano | Seismic precursory patterns before a cliff collapse and criticalpoint phenomena[END_REF]. Considering that such a pattern fits with the critical point theory, the authors suggested that the apparent b-value decrease could be induced by a cut-off of increasing size as the failure approaches rather than by a real change of the power-law exponent.

For brittle rocks, the b-value has been reported to depend on the mechanical loading conditions, i.e. type of faulting for Earthquakes [START_REF] Schorlemmer | Variations in earthquake-size distribution across different stress regimes[END_REF], the confining pressure [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b-value[END_REF], or the roughness of sliding surfaces [START_REF] Sammonds | Evolution of microseismicity during frictional sliding[END_REF] in the laboratory. All these observations agree with a dependence of the b-value on friction and/or mean stress. The distinction between these two factors remains difficult because they are observed to be correlated [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b-value[END_REF], as friction decrease when the confining pressure increases [START_REF] Jaeger | Fundamentals of rock mechanics[END_REF].

Spatial variations of the b-value have been studied by calculating this parameter on small areas for tracking spatial variations of either the stress or the mechanical properties of the Earth's crust (e.g. [START_REF] Montuori | Crustal heterogeneity highlighted by spatial b-value map in the Wellington region of New Zealand[END_REF][START_REF]Mapping spatial variability of the frequencymagnitude distribution of earthquakes[END_REF][START_REF] Wyss | Mapping active magma chambers by b values beneath the Off-Ito volcano, japan[END_REF]). A decrease of the b-value with increasing depth has been observed for earthquakes in California [START_REF] Mori | Depth dependence of earthquake frequency-magnitude distributions in california: Implication for rupture initiation[END_REF][START_REF] Gerstenberger | A systematic test of the hypothesis that b varies with depth in california[END_REF]. Similar depth dependence has been observed in the western Alps [START_REF] Sue | Mechanical behavior of western alpine structures inferred from statistical analysis of seismicity[END_REF] and for induced seismicity related to Aswan Lake in Egypt [START_REF] Mekkawi | A Long-Lasting relaxation of seismicity at Aswan reservoir, Egypt, 1982-2001[END_REF] indicating a possible relationship between b-value and the mean stress and/or the brittle/ductile transition. Schorlemner et al [START_REF] Schorlemmer | Microseismicity data forecast rupture area[END_REF] observed that the M=6 Parkfield earthquake fits well an area of particularly low b-value, suggesting that this parameter could be used to identify areas prone to the failure.

A few studies display variations in the shape of the distribution tails, showing that the power-law trend is not always respected for the largest events of the catalogue (e.g. [START_REF] Main | Apparent breaks in scaling in the earthquake cumulative Frequency-Magnitude distribution: Fact or artifact?[END_REF]). The distribution tail can be characterized alternatively by a lack (cut-off e.g. [START_REF] Sue | Mechanical behavior of western alpine structures inferred from statistical analysis of seismicity[END_REF] ) or by an excess (oversizing e.g. [START_REF] Bossu | Complexity of intracontinental seismic faultings: The gazli, uzbekistan, sequence[END_REF]) of large events compared with the powerlaw trend. This discrepancy with the power-law trend is of major importance regarding the ability to estimate the probability of the largest events by extrapolation of the observed distribution for the smaller ones.

Numerical simulations

Many studies use numerical simulation in order to recover power-law distributions of failure events (see [START_REF] Alava | Statistical models of fracture[END_REF] for a review). A common outcome of these simulations is that power-law distributions emerge spontaneously from elastic interaction within heterogeneous media [START_REF]Emerging complexity in a simple model of the mechanical behavior of rocks[END_REF]. Early models based on cellular automata simulating the behavior of a network of sliding blocks connected by springs [START_REF] Burridge | Model and theoretical seismicity[END_REF][START_REF] Chen | Self-organized criticality in a crack-propagation model of earthquakes[END_REF][START_REF] Olami | Self-organised criticality in a continuous, nonconservative cellular automaton modeling earthquake[END_REF] aimed to reproduce the behavior of a fault. These models reach a metastable state characterized by a power-law distribution of avalanches with an exponent remaining constant. The value of the exponent may change depending on the rules of interaction and of the non-conservative part of the energy balance.

Fiber bundle and lattice solids models, with of without electrical-mechanical analogy, have been extensively used for characterising the emergence of power-law in rupture phenomena, particularly as the macro-failure is approached (e.g. [START_REF] Pradhan | Dynamic critical behavior of failure and plastic deformation in the random fiber bundle model[END_REF][START_REF] Zapperi | Fracture statistics in the threedimensional random fuse model[END_REF][START_REF] Zapperi | Plasticity and avalanche behaviour in microfracturing phenomena[END_REF]). Due to its numerical simplicity, this type of model allowed the authors to investigate the behavior of large size models and in particular to identify the role of the finite size of the system and of the distance to the critical point on the emergence of power-law distributions. They showed in particular that power-law distributions of rupture avalanches are fully developed only in the vicinity of the critical point, i.e. the macrofailure, as described in critical point theory [START_REF] Sornette | Critical phenomena in natural sciences, chaos, fractals, self organization and disorder: Concepts and tools[END_REF]. Variations in the power-law shape and tail are shown to be related to the presence of an exponential cut-off that depends on the finite size of the model and diverges as the critical point is approached. The finite rate of the loading has also been proposed to induce a cut-off of the power-law distribution for largest size events [START_REF] Rosso | Depinning of elastic manifolds[END_REF][START_REF] Picallo | Energy dissipation statistics in the random fuse model[END_REF]. The evoked mechanism is that the loading rate should be infinitely small, or small enough, for allowing the long range interaction that permits the emergence of avalanches at all scales. In the case of a too fast loading, viscous dissipation could contradict the full range elastic interaction, limiting the emergence of largest rupture events. Departure from pure power-law may also originate from avalanches that overlap.

Models including more realistic rules of stress redistribution, energy dissipation, strain softening or hardening (e.g. [START_REF] Ben-Zion | Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive evolutionary changes, and different dynamic regimes[END_REF][START_REF] Ka Dahmen | Micromechanical model for deformation in solids with universal predictions for Stress-Strain curves and slip avalanches[END_REF]) allow investigation of the conditions of appearance of cut-off or characteristic size failure events. More complex models including tensorial description of the stress state have been proposed [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b-value[END_REF][START_REF] Amitrano | From diffuse to localized damage through elastic interaction[END_REF][START_REF] Tang | Numerical simulation of progressive rock failure and associated seismicity[END_REF] that simulate both the progressive damage of brittle solids, the damage localization process and, for some of them, the power-law distributions of rupture events. This type of model has been shown to produce critical behavior associated with damage localization [START_REF] Girard | Failure as a critical phenomenon in a progressive damage model[END_REF]. The tensorial description of the elastic interaction allowed to draw the connection between the power-law distributions of rupture and the scale invariance of the strain field. This argues for the importance of long range correlation in the emergence of power-law distribution of rupture events.

Potential explanations for power-law variations

Incorrect statistical methods

The main statistical bias is obviously related to the number of events, i.e. the size of the catalogue, used for calculating the exponent of the power-law distribution. Many studies, essentially field works, estimated the b-value on samples as small as a few hundred events. This is generally due to the difficulty of reaching a large number of well-identified seismic events in certain field conditions. Systematic tests based on random sampling from powerlaw distribution [START_REF] Pickering | Sampling power-law distributions[END_REF] have shown that the confidence of the b-value estimated by least squares regression decreases as the sample size increases. It becomes lower than 0.1 for a sample larger than 200. We should note that studies performed in the framework of statistical physics are generally based on several thousand events thus avoiding sampling biases.

The second bias is related to the method used for fitting a power-law trend onto the event size distribution. In the framework of statistical physics, the powerlaw trend is generally estimated by a linear regression in a log-log plane. When restricted to the power-law trend, this method provides an objective estimation of the exponent giving the same weight to all the size classes, without taking into account the number of events in each class. More elaborate methods combining maximum-likelihood fitting methods with goodness-of-fit tests based on the Kolmogorov-Smirnov statistics have been proposed to estimated both the exponent and the range of the power-law trend [START_REF] Clauset | Power-Law distributions in empirical data[END_REF]. In the field of Earth Sciences the maximum likelihood method (MLM) proposed by Utsu and Aki [START_REF] Aki | Maximum likelihood estimate of b in the formula, log N = a-bM, and its confidence limits[END_REF][START_REF] Utsu | A method for determining the value of b in a formula logn=a-bm showing the magnitude-frequency relation for earthquakes[END_REF] is frequently used for estimating the b-value. b = log 10 e < M > -M min [START_REF] Amorese | On varying bvalues with depth: results from computer-intensive tests for southern california[END_REF] Where < M > and M min are the arithmetic mean and the minimum magnitude of the catalogue. As it uses the mean magnitude, this method provides an estimate controlled essentially by smallest events, i.e. the most numerous. Hence it is less sensitive to the occurrence of largest events than the least squares method. The essential assumption for the validity of the MLM estimate is that M min is above the completeness magnitude [START_REF] Wiemer | Minimum magnitude of completeness in earthquake catalogs: Examples from alaska, the western united states and japan[END_REF] (i.e. the magnitude above which the events are detected without sampling biases e.g. regarding the distance to sensors), and that a power-law trend is observed in the full range of magnitudes. In the case where a cut-off exists for the largest events (e.g. in Fig. 2), this method leads to an underestimation of the b-value. The MLM is often associated with the so-called Utsu's test [START_REF]A statistical significance test of the difference in b value between two earthquake groups[END_REF] for determining if the difference in b-value observed between two samples is significant. Recently authors [START_REF] Amorese | On varying bvalues with depth: results from computer-intensive tests for southern california[END_REF] used statistical tests to compare least squares and maximum likelihood methods on both natural and synthetic catalogues for estimating biases related to the sampling. They concluded that the Utsu's test gives a biased indication of the significance of a difference in b-value between two catalogues. In particular for small size samples, this test provides a too easy acceptance of the difference as significant. They considered that many previous studies associated the Utsu's test with catalogues of insufficient numbers of events, leading erroneously to consider insignificant variations of the b-value to be significant. According to their study, a sample size of 500 events is the minimum for identifying variations in b-value of about 0.1 with a good statistical reliability.

A third bias can exist when comparing cumulative (cdf ) and discrete (pdf ) event size distributions. In the case of a pdf following a pure power-law, the integration preserves the power-law shape with an exponent shifted by -1 (see equ. 1 and 2). To be correct the pdf must be calculated by dividing the number of events in each bin by the bin width. Thus the pdf is independent of the binning, whether log or linearly spaced. A potential problem may arise when calculating the pdf for magnitude frequency and referring to the Gutenberg-Richter law (equ. 6). As the Gutenberg-Richter law corresponds to a linear trend in a semilog plot, the corresponding cdf is often calculated with a linearly spaced binning of the magnitude. Because the bin width, measured in magnitude, is constant, the division by the bin width preserve the b-value (see Fig. 1). If the comparison between pdf and cdf is necessary, this problem can be avoided by converting Magnitude into Energy using equ. 7 and treating the frequency distribution as a power-law. Figure 1 presents an illustration of the difference between cumulative or discrete distribution for either Magnitude or Energy. For Magnitudes the distribution are expressed in event number (N (> M) and N (M ) respectively) respecting the usage in Earth Science (e.g. [START_REF] Schorlemmer | Variations in earthquake-size distribution across different stress regimes[END_REF]). It appears clearly that for magnitude the b-value is the same for the cumulative and discrete distributions. On the contrary, for the Energy the difference of -1 between the exponents of the cdf and the pdf is respected.

In the case where a discrepancy with a power-law appears, e.g. cut-off due to the lack of largest events, the integration induces a curvature in the whole distribution that can mask the real power-law segment. This can be seen on fig. 2 where the power-law trend of the pdf appears only between the completeness level, for smallest events, and the cut-off, for largest event (10 3 > E > 10 8 ). These two limits of the power-law trend appear less clearly on the cdf and are shifted toward smaller values. The pdf is preferable in such a case to avoid associated errors in the estimation of the power-law range.

Finite size effect

The emergence of a pure power-law (i.e. without cut-off) is theoretically possible only within a system of infinite size. In natural conditions such a hypothesis is obviously not found to be true, even if in some cases (e.g. for the whole Earth's crust) the size is large enough to be considered as infinite. In the case of a finite size system, the occurrence of the largest events is constrained by the size of the system. As a consequence the power-law distribution is affected by a exponential tail. Thus the equ. 2 describing the pdf (or equ. 1 if the cdf is considered) should be modified for taking into account the cut-off, e.g. using an exponential decay".

p(s) ∝ s -β-1 exp(-s/s 0 ) ( 9
)
Where s is the size of the event and s 0 is the cut-off size. It has been proposed that the cut-off size s 0 scales with the size of the system L obeying a power law s 0 ∝ L δ [START_REF] Alava | Statistical models of fracture[END_REF][START_REF] Girard | Failure as a critical phenomenon in a progressive damage model[END_REF][START_REF] Sornette | Critical phenomena in natural sciences, chaos, fractals, self organization and disorder: Concepts and tools[END_REF]. The value of exponent δ depends on the dimensionality and on the dynamics of the system. It can be easily determined in the case of numerical simulations by changing the system size L and analysing the size dependence of s 0 . In the case of natural systems, the identification of the finite size effect is more difficult as it can be confused with the cut-off related to divergence towards the failure (see section 2.3.3) or the one related to the finite loading rate (see section 1). Mining-induced seismicity provides a good example of cut-off probably induced by finite size effect of by finite loading rate. In particular if the data set includes all the events related to the collapse process, the effect of distance to the failure can be excluded (see the section 2.3.3). Figure 2 presents microseimic events recorded in a French coal mine (see [START_REF] Senfaute | Spatial distribution of mining tremors and the relationship to rockburst hazard[END_REF] for more details) for all events recorded during the excavation process up to and including the final rupture induced by the mine collapse. According to the critical behavior expected for this type of brittle failure process, the distribution of the whole catalogue should obey a pure power-law [START_REF] Sornette | Critical phenomena in natural sciences, chaos, fractals, self organization and disorder: Concepts and tools[END_REF]. The cut-off for largest events is more likely to be due to the finite size of the mine that can not induce events as large as determined by the extension of the power-law trend.

Alternatively, the loading induced by the cutting at the front of the excavation can be considered as fast enough for activating viscous dissipation inducing a cut-off related to finite loading rate [START_REF] Picallo | Energy dissipation statistics in the random fuse model[END_REF][START_REF] Rosso | Depinning of elastic manifolds[END_REF].

Proximity of failure

Since the work of Scholz [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF], the applied stress is the main factor evoked for explaining the change in b-value. This is based on the common laboratory observation where the increase of the applied stress is accompanied by the b-value decrease [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b-value[END_REF][START_REF] Guarino | An experimental test of the critical behaviour of fracture precursors[END_REF][START_REF] Lockner | The role of acoustic emission in the study of rock fracture[END_REF][START_REF] Meredith | Temporal variation in seismicity during quasi-static and dynamic failure[END_REF][START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF], the minimal value being reached just before the macrofailure. This dependency has been used as a proxy for estimating the stress in the Earth's crust, in particular in studies analysing the spatial variability of the b-value [START_REF] Gerstenberger | A systematic test of the hypothesis that b varies with depth in california[END_REF][START_REF] Schorlemmer | Microseismicity data forecast rupture area[END_REF][START_REF] Urbancic | Space-time correlations of b-values with stress release[END_REF][START_REF]Mapping spatial variability of the frequencymagnitude distribution of earthquakes[END_REF][START_REF] Wyss | Mapping active magma chambers by b values beneath the Off-Ito volcano, japan[END_REF][START_REF] Mori | Depth dependence of earthquake frequency-magnitude distributions in california: Implication for rupture initiation[END_REF].

Considering that the b-value decreases before the macrofailure, several authors have used this parameter for forecasting failure events in natural conditions [START_REF] Amitrano | Seismic precursory patterns before a cliff collapse and criticalpoint phenomena[END_REF][START_REF] Smith | The b-value as an earthquake precursor[END_REF][START_REF] Tsukakoshi | Decreased b-value prior to the m 6.2 northern miyagi, japan, earthquake of 26 july 2003[END_REF][START_REF] Vinciguerra | Damage mechanics preceding the September-October 1998 flank eruption at Mount Etna volcano inferred by seismic scaling exponents[END_REF][START_REF] Vinciguerra | The two behaviours of mt. etna volcano before and after a large intrusive episode: evidences from b value and fractal dimension of seismicity[END_REF]. In the case of objects for which the macrofailure can be defined clearly, e.g. landslides, cliffs or volcanoes, this approach is relevant as the collapse can be considered as the end of the progressive failure process. More problematic is the case of earthquakes that can be instead considered as part of a marginally stable process that never ends. Here the analysis needs to distinguish the macro-events, i.e. the so-called characteristic earthquakes, from the background seismic activity [START_REF] Tsukakoshi | Decreased b-value prior to the m 6.2 northern miyagi, japan, earthquake of 26 july 2003[END_REF]. The Gutenberg-Richter law, when faithfully respected, contradicts this distinction as no breakdown appears in the scaling law. So there is no reason to distinguish events of a particular size. Alternatively, if a characteristic size appears in the magnitude distribution, the former approach could make sense. The case of the California catalogue presented in fig. 1 provides a possible application of such separation between large and small events as a characteristic size appears on the magnitude distribution for M larger than 6. We calculated the b-value for successive earthquakes sequences occurring in California in the 1990-2010 period (117100 events). We used a sufficiently large number of events (1000) contained in moving windows to allow a robust estimate of the bvalue regarding the sampling bias. Figure 3 shows that significant temporal variations of the b-value exist. In order to verify if this decrease may be considered as a forerunner of large earthquakes, we identify the largest events considering different magnitude thresholds (M>7, M>6, M>5.5). One may first observe that the three M>7 events are preceded by a substantial decrease of the b-value. This could be considered as a promising result as the rate of detection is good. But several times there is a decrease of the b-value without the occurrence of such a large event so that the rate of false alarms is not negligeable. When decreasing the threshold for large events to M>6, the rate of false alarms decreases but the rate of detection too. This effect is enhanced when considering events of M>5.5. The rate of detection failure and false alarms is so high that the forecast ability of the b-value decrease becomes doubtful. An interesting point is that the magnitude distribution (Fig. 1) for events larger than M=6, displays a slight discrepancy compared to the power-law. So that these larger events could be considered as characteristic earthquakes. This could explain why a precursory behavior appears before these events, but does not explain the occurrence of false alarms. The observed decrease of the b-value before large events has been alternatively explained by epidemic models showing that foreshocks, with decreasing b-value, can result from cascades of triggered earthquakes [START_REF] Helmstetter | Foreshocks explained by cascades of triggered seismicity[END_REF]. But the characteristic size we observe is not explained by this kind of model.

The stress applied on the material has been early proposed to be an explanative parameter for the b-value decrease observed before the failure [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF]. Indeed the failure occurs when the stress is increased up to the strength of the material. So the stress has been considered as a measure of the remaining distance before the failure. The absolute value of the material strength may vary strongly depending on the considered material. Therefore for allowing comparisons between various materials, the absolute value of the stress should be normalized by the strength (i.e. the stress value at failure) of each material. The critical point theory rationalizes this statement by defining a control parameter normalizing the distance to the failure. The event size distribution can be expressed as a function of the size s and of the control parameter Δ. p(s, Δ) ∝ s -τ .exp(-s/s 0 ) [START_REF] Bossu | Complexity of intracontinental seismic faultings: The gazli, uzbekistan, sequence[END_REF] With Δ = σc-σ σc for stress controlled loading or Δ = c- c for strain controlled loading. σ c and c are the stress and strain at the critical point, i.e. the macro-failure. In this description, power-law distribution of the rupture events is associated with an exponential cut-off. According to the critical point theory, the cut-off should diverge when approaching the critical point, whereas the β exponent should remain constant. The cut-off s 0 divergence when Δ → 0 has been proposed to be well described by a power law s 0 ∝ Δ -γ [2,17,49, chap.9].

Including the exponential tail in the range used for fitting the power-law trend, as Δ decreases, could display an apparent decrease of the exponent (see e.g. [START_REF] Amitrano | Seismic precursory patterns before a cliff collapse and criticalpoint phenomena[END_REF]). This could be the origin of the numerous observations of the apparent b-value decrease before the failure. Considering the mechanism of sweeping an instability, Sornette [START_REF] Sornette | Sweeping of an instability : an alternative to self-organized criticality to get powerlaws without parameter tuning[END_REF] has analytically demonstrated that the distribution of the whole population of events, i.e. including the progressive evolution of the control parameter toward zero, obeys a power-law conforming to the Equ. 9 (see section 4 for a thorough explanation). In such a case, the cut-off is only related to finite size effect.

Summary

Summarizing the preceding review, excluding incorrect statistics methods, a few candidates arise for explaining variations of the power-law distribution of rupture events observed in natural conditions, experimental and numerical studies.

-Heterogeneity: based on theoretical considerations of the probability of failure in a heterogeneous material, Mogi [START_REF] Mogi | Magnitude frequency relations for elastic shocks accompanying fractures of various materials and some related problems in earthquakes[END_REF] has proposed that the b-value should increase as the heterogeneity increases.

-Stress: Scholz [START_REF] Scholz | The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes[END_REF] has first observed a decrease of the bvalue when increasing the stress applied on rock samples. -Proximity of failure: following the previous statement, the b-value decrease has been considered to be a forerunner of the failure. This behavior has been observed prior to some large failure events but remains controversial for large earthquakes.

-Finite size effect: the size of the largest event is constrained by the finite size of the object concerned by the failure. This has been widely and essentially investigated by numerical models of various size. The finite value of loading rate as also been proposed to induce such a cutoff.

-Divergence of the cut-off toward the failure: the exponential tail of power-law is characterized by a cut-off whose size diverges as the failure approaches, the failure is considered to be a critical point. The latter two points can lead to a biased estimate of the b-value if the range used for estimating it includes the exponential tail.

In the next section we use numerical modeling to evaluate the pertinence of some of the stated potential explanations for the observed variability of power-laws in failure phenomena. In particular we focus on the effect of mechanical properties and on the proximity of the failure.

Numerical modeling of progressive damage 3.1 Model features

The model we use here has demonstrated interesting features for investigating power-law emergence in the progressive failure process [START_REF] Amitrano | Brittle-ductile transition and associated seismicity: Experimental and numerical studies and relationship with the b-value[END_REF][START_REF]Failure by damage accumulation in rocks[END_REF][START_REF] Amitrano | From diffuse to localized damage through elastic interaction[END_REF][START_REF] Girard | Failure as a critical phenomenon in a progressive damage model[END_REF]. It includes elastic interaction within a heterogeneous material with simple damage rules. It is based on progressive isotropic elastic damage. The effective elastic modulus, E ef f , is expressed as a function of the initial modulus E 0 and damage D:

E ef f = E 0 (1 -D).
Such a relation works when the considered volume is large compared with the defect size, such as cracks, and then can be considered as a mesoscale relationship. The damage parameter D has been proposed to be related to crack density (see [START_REF] Kemeny | Effective moduli, non-linear deformation and strength of a cracked elastic solid[END_REF] for a review). The simulated material is discretized using a 2D finite element method with plane strain hypothesis. The loading consists of increasing the vertical displacement of the upper model boundary. When the stress in an element exceeds a given damage threshold, its elastic modulus is multiplied by a factor (1 -D), D being constant and small compared to 1 (here we used D=0.1). Because of the elastic interaction, the stress redistribution around a damaged element can induce an avalanche of damages that we call an event. The total number of damaged elements during a single loading step is the avalanche size s, which is comparable to the acoustic emission event amplitude observed in laboratory experiments. The Mohr-Coulomb criterion is used as a damage threshold:

σ S = μ.σ N + C ( 11
)
where σ S is the shear stress; σ N is the normal stress; C is the cohesion; and μ is the internal friction coefficient. This criterion reflects that the maximum shear stress σ S the material can sustain increases linearly with the normal stress. Thus it applies well for so-called frictional brittle materials as rocks or ice [START_REF] Jaeger | Fundamentals of rock mechanics[END_REF][START_REF] Schulson | On the initiation of shear faults during brittle compressive failure: A new mechanism[END_REF]. To obtain macroscopic behaviors differing from those of the elements, damage localization and avalanches size spreading over several magnitude orders, it is necessary to introduce heterogeneity. To simulate material strength heterogeneity, the cohesion of each element, C, is randomly drawn from a uniform distribution. Previous studies have shown that the absolute value of strength doesn't influence the behavior of the model but only changes the scale of the stress/strain relationship. On the opposite, the internal friction μ has a strong influence on the damage localization, the brittle/ductile type of macroscopic behavior and the distribution of avalanches size [START_REF]Failure by damage accumulation in rocks[END_REF]. Figure 4 shows the macroscopic strain-stress relationship obtained for simulations with internal friction μ ranging from 0 to 1. The macroscopic behavior changes from ductile to brittle, as observed in laboratory experiments for materials characterized by low/high friction. The effect of the friction on the macroscopic behavior has been discussed in depth in previous works [START_REF]Emerging complexity in a simple model of the mechanical behavior of rocks[END_REF][START_REF]Failure by damage accumulation in rocks[END_REF][START_REF] Amitrano | From diffuse to localized damage through elastic interaction[END_REF] and has been shown to be related to the interaction geometry between defects that become more anisotropic as the friction increases. This induces a more localized damage (Fig. 5). In the following we focus on the effect of μ on the avalanches distribution.

Effect of friction on the rupture events size distributions

For the purpose of the present work we carried out simulations with various internal friction values focusing essentially on their impact on the avalanches size distribution. Mechanical parameters have been chosen to correspond to rocks. The used parameters are the following: initial elastic modulus = 50 GPa, Poisson's ratio = 0.25, cohesion ranging from C min = 25 MPa to C max =50 MPa, internal friction ranging from 0 to 1. The model size was 128x256 elements. Loading was strain-controlled, i.e. the upper boundary displacement was progressively increased.

Figure 4 shows the macroscopic behavior obtained by numerical simulations with internal friction varying from 0 to 1. The corresponding damage state at the end of the simulation is shown on Figure 5, for μ=0, 0.5 and 1 respectively. Increasing the internal friction μ induces a progressive change from ductile to brittle behaviour characterised by the occurrence of a macroscopic stress drop of increasing size. In the case of internal friction very low, i.e. near zero, the stress reaches a metastable value without appearance of a particular point corresponding to the stress peak. So no macro-failure can be individuated in this case. We consider that the macro-failure corresponds to the occurrence of the largest stress drop when it is large enough for affecting significantly the stress-strain curve. In this case it corresponds to the stress peak. The macro-failure associated with the macro-scale stress drop is considered as a critical point in the following.

Regarding the events size distribution including the whole simulation (Fig. 6), one may observe that the internal friction influences both the exponent (β in Equ. 9) of the power-law trend and the tail shape. The exponent increases with the friction, whereas the tail shape progressively evolves from a cut-off, i.e. lack of large events, to an oversizing, i.e. occurrence of extreme events out of scale. Hence the macro-failure event generally does not scale with the power-law trend. It is over-or under-estimated depending on the ductile/brittle type of behaviour. In order to characterise the discrepancy between the maximum event size extrapolated from the power-law trend and that occurring during the simulation, we calculated the ratio R = max(s obs. ) max(sPL) , max(s obs. ) being the maximal observed size during the simulation and max(s P L ) being the maximal size estimated by the power-law considering the same probability of occurrence. Table 1 shows the values of β and R obtained for the different values of μ. β slightly depends on μ which could explain spatial variations of the b-value observed for earthquakes. More pronounced is the impact of μ on R that varies from values smaller than 1, corresponding to a cut-off, to values larger than 1, corresponding to extreme event compared with the power-law trend. Only in the case of μ=0.5, does the power-law trend provide a good estimate of the maximum event size, as R is close to 1. An interesting point is that the mean value of friction for the Earth's crust material is estimated to be in the range 0.5-0.6 for seismogenic conditions [START_REF] Byerlee | Friction of rocks[END_REF], including the former value. This could explain why cut-off or oversizing are rarely observed for earthquakes. To this concern, earthquakes can be considered as predictable in size.

Critical behavior toward the failure

In order to characterize how the events size distribution evolves toward the macro-failure, we calculated p(s) for successive bins of the control parameter Δ. As the simulation are realized under strain control condition, Δ = c- c , c being the strain at the macro-failure. The binning has been equally spaced in a log scale and sized for ensuring a large enough number of avalanches in each bin. For each bin the representative value of the control parameter has been calculated as the center of the bin in a log scale. For ensuring good statistical representativity we stacked 20 simulations for each configuration. Figure 7 shows the distribution obtained for different values of Δ in the case μ=0.5. One can observe that the distributions are affected by a cut-off whose size increases with decreasing control parameter. The power-law trend appears fully only for the smallest value of Δ. The τ exponent (see Equ. 10) has been estimated for the smallest value of the control parameter Δ to be 1.5±0.05. In order to characterize the evolution of the cut-off on approaching the failure, we performed a data collapse conforming to Equ. 10. The best collapse has been obtained for τ = 1.45 ± 0.05 and γ = 2 ± 0.1. We realized the same data collapse for simulations performed with different values of the internal friction. We excluded the value μ=0 as no macroscopic failure can be distinguished on the strain-stress curve for this value (see Fig. 4) so no critical point can be defined. The values of the exponents resulting from this analysis are summarised in Table 1. We observed a slight change for the τ exponent, ranging from 1.4 to 1.55. The impact is more pronounced for the γ exponent, as it varies from 2.5 to 1.3.

This demonstrates the good agreement between the simulations results and the critical point theory. Accordingly the change of events size distribution toward the failure can be fully explained by the divergence of the cut-off as the the failure is approached. But this behavior can be confused with a decrease of the τ exponent if the exponential tail is included in the estimation of the exponent. The discrepancy with the often observed decrease of the b-value, at least for natural objects, may originate in methodological biases. The analysis we developed requires a large number of events, especially near the failure. This can be difficult to obtain in natural conditions. In the case of catalogues of limited size (i.e. a few hundreds) the identification of the exponential tail can be hard to achieve (e.g. [START_REF] Amitrano | Seismic precursory patterns before a cliff collapse and criticalpoint phenomena[END_REF]). Moreover, this kind of analysis requires the definition of a well identified control parameter. In the case of natural objects, the control parameter can be difficult to define as the failure may occur under conditions mixing strain and stress (e.g. for crustal earthquakes) or under constant stress (e.g. static fatigue for landslides or rock collapse). When focusing on the short period prior to the failure, the power-law fit provides directly an estimate of the τ exponent (see fig. 7). In our analyze we observed that this exponent is slightly dependent on the mechanical properties of the material. This could explain the spatial variations of the b-value observed in numerous studies by variations of material properties, when focusing on events occurring just before the macro-failure.

Discussion and Conclusion

We first reviewed observations showing that a large variability exists in the size distribution of rupture events. Many of these reports concern natural observations of Earthquakes, at different scales and in different conditions, and acoustic emissions recorded at the laboratory. The variability concerns essentially the power-law exponent that has been observed to decrease before the failure, essentially at the laboratory but also before some large Earthquakes. Spatial variability of the exponent has been also observed. This has been interpreted as the result of variation of either the material properties, as friction and heterogeneity, or of stress. Some rare studies also stated variations in the tail shape identifying either cut-off or oversizing. Besides possible methodological biases, a few potential explanations raised from this review for explaining the observed variability of power-laws.

Theoretical considerations referring to both critical point theory [START_REF] Alava | Statistical models of fracture[END_REF] and sweeping of an instability [START_REF] Sornette | Sweeping of an instability : an alternative to self-organized criticality to get powerlaws without parameter tuning[END_REF] leads to consider that the events size distribution of the rupture event should be fully described using three exponents instead of a single one. The b-value referring to the Gutenberg-Richter law could be related to different theoretical exponents depending on how it is calculated.

Critical point theory refers to two exponents τ and γ describing the power-law distribution and the divergence of the cut-off size toward the failure, respectively. These two exponents as supposed to be constant for a given material and loading mode. The b-value could be related to the exponent τ only in the vicinity of failure when the value of the exponential tail becomes large enough to not alter the real power-law trend. The decrease of the b-value often observed before the failure should be then considered as an artifact related to the exponential tail of increasing size that is included in the power-law fit. However this apparent decrease provides precursory patterns before the failure useful for the risk assessment. The numerical simulations carried out in Section 3 have shown that τ and γ could depend on the mechanical properties of the material. In particular we displayed the effect of the internal friction that controls the brittleness of the macroscopic behaviour.

The mechanism of sweeping of an instability [START_REF] Sornette | Sweeping of an instability : an alternative to self-organized criticality to get powerlaws without parameter tuning[END_REF] allows to explain that summing the entire set of events of the progressive failure sequence, i.e. including the complete divergence of the exponential cut-off, a power-law distribution is recovered conforming to Equation 1 or 9, with the relation β = τ -1 + 1/γ (note that the term -1 comes from the fact that β is related to the cdf , contrarily to the original work of [START_REF] Sornette | Sweeping of an instability : an alternative to self-organized criticality to get powerlaws without parameter tuning[END_REF]). The former relation is well respected by the exponents estimated for numerical simulations (see Tab. 1). Thus, when considering a large set of events including a notable variation of the control parameter (e.g. a long temporal serie of earthquakes or the full set of AE events recorded during an experiment) the b-value can be related to the β exponent. In this case, the numerical simulation outcomes showed that the β exponent depends also on the mechanical properties of the considered material. This could explain the spatial variations of the b-value observed in different sites (e.g. [START_REF]Mapping spatial variability of the frequencymagnitude distribution of earthquakes[END_REF]) and also the variability of the β observed at the laboratory (e.g. [START_REF] Picallo | Energy dissipation statistics in the random fuse model[END_REF]).

Moreover, the simulations displayed a dependence between the β exponent and the brittleness of the behaviour. The change in β is associated with a shift from a cut-off to an oversizing characterized by the occurrence of extreme events. The intermediary case for which it is possible to forecast the largest event using the power-law trend of the smallest ones is obtained for μ=0.5 what corresponds roughly to the mean friction value for the Earth's crust. This suggests that earthquakes could be in the range of predictable events in size. For high internal friction corresponding to highly brittle behavior the largest event, i.e. the size of the macro-failure, is out of scale compared with the power-law trend, we characterized by a high value of the R ratio. Referring to the topic of this special issue, the ductile-brittle transition could be considered as the progressive appearance of Dragon-Kings, i.e. extreme events. Despite their unpredictability in size, considering they obey the critical point theory, these particular events remain predictable in time, at least if a temporal evolution of the control parameter is known, as they emerge from the divergence of the failure process.

The author thanks L. Girard and the two anonymous reviewers whose critical remarks helped to clarify this work, and H. Dyer for improvements of the manuscript. All (or most of) the computations presented in this paper were performed at the Service Commun de Calcul Intensif de lObservatoire de Grenoble (SCCI-CIMENT).

μ β R τ γ τ-1 + 1/γ ±0.

Fig. 3 .

 3 Fig. 3.Temporal evolution of the b-value calculated for successive windows of 1000 events with an overlap of 500 events, in California from 1990 to 2010. The 95% confidence interval is given by dotted lines. The catalog includes 117100 Earthquakes ranging from M=2 to M=7.3. The b-value is estimated for the range M=2 to M=5 by least square regression. Maximum likelihood method gives comparable results. Triangles indicates the occurrence time for Earthquakes larger than magnitude 5.5, 6 and 7.

Fig. 1 .

 1 Fig. 1. Earthquake distribution in California from 1990 to 2010 (extract from the ANSS catalog). Top: Cumulative N (> M ) and discrete N (M ) distribution of Magnitude. Bottom: Cumulative distribution cdf and probability density functions pdf of the Energy. The pdf is normalized for plotting convenience.

Fig. 2 .

 2 Fig. 2. Earthquake energy distribution of microseismicity induced by mining excavation (courtesy of G. Senfaute, Ineris).

Fig. 4 .

 4 Fig. 4. Macroscopic behavior obtained by numerical simulations for different internal friction values (H=0.5, μ varies from 0 to 1).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Damage state of the numerical samples at the end of the simulation, for μ=0, μ=0.5, μ=1

Fig. 7 .

 7 Fig. 7. Event size distribution obtained by numerical simulations for different range of control parameter Δ (indicated in the legend). Simulations have been realized with (μ=0.5). The dotted guide line corresponds to an exponent -tau-1 = -1.5.

Fig. 8 .

 8 Fig. 8. Data collapse of the probability distributions of event size s for different values of the control parameter (indicated in the legend). The best collapse is obtained with γ=1.75±0.05 and β=0.45±0.05. The insert reproduces the event size distributions showed in figure 7.

Table 1 .

 1 Exponents characterizing the event size distribution for various internal friction μ determined form numerical simulations outcomes. β is the exponent characterising the size distribution of the whole set of events during the simulation. R is the oversizing ratio as defined in section 3.2. τ and γ are the exponents characterizing the critical behavior toward the failure (see equ. 10). τ -1 + 1/γ is given for comparison with β.

		05		±0.05 ±0.1	±0.15
	0	0.85	0.004	-	-	-
	0.25	0.95	0.055	1.4	2.5	0.8
	0.5	1	1.06	1.45	1.7	0.95
	0.75	1.2	8.7	1.55	1.6	1.15
	1	1.3	20.9	1.6	1.3	1.3