Laure Pizzella
email: laure.pizzella@gmail.com

Robin Alais
email: robin.alais@mines-paristech.fr

Simon Lopez
email: slopez@brgm.fr

Xavier Freulon
email: xavier.freulon@mines-paristech.fr

Jacques Rivoirard
email: jacques.rivoirard@mines-paristech.fr

Taking Better Advantage of Fold Axis Data to Characterize Anisotropy of Complex Folded Structures in the Implicit Modeling Framework

Keywords: 3-D structural geological modeling, Implicit modeling, Potential-eld method, Non-Stationarity, Cokriging 3-D structural geological modeling, Implicit modeling, Potential-field method, Non-Stationarity, Cokriging

come   L'archive ouverte pluridisciplinaire

Introduction

Structural geological modeling aims to achieve three-dimensional models of geological objects and of their relationships. It is a mandatory and determining step for most geological studies in various fields of application: resource estimation, reservoir study, risk assessment, environmental remediation. It relies on a conceptual model proposed by geologists based on field observations. Given these preliminary data, the main goal is to choose a suitable mathematical model with parameters that give an optimal representation of the reality [START_REF] Wellmann | 3-d structural geological models: Concepts, methods, and uncertainties[END_REF].

For this purpose, explicit and implicit modeling are the two leading approaches. The former directly builds geological interfaces by solving a least-square smooth approximation problem on a discrete mesh [START_REF] Mallet | Discrete smooth interpolation in geometric modelling[END_REF]. In contrast, the latter relies on a preliminary interpolation of a three-dimensional scalar field from geological data. Only then are the surfaces of interest extracted: they are those holding the same scalar value, referred to as their potential.

Over the last decade, implicit modeling algorithms have gained popularity in geomodeling software. Two families of methods can be distinguished, depending on whether the interpolation is performed on a meshed grid [START_REF] Frank | 3d-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF][START_REF] Caumon | Three-dimensional implicit stratigraphic model building from remote sensing data on tetrahedral meshes: Theory and application to a regional model of la popa basin, ne mexico[END_REF] or not [START_REF] Houlding | 3d geoscience modeling: computer techniques for geological characterization[END_REF]. Focus here is on the meshless approach. In this context, different basis functions have been deployed in order to interpolate the scalar field: the Radial-Basis functions [START_REF] Cowan | Practical implicit geological modelling[END_REF][START_REF] Turk | Modelling with implicit surfaces that interpolate[END_REF][START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions[END_REF], the Moving Least Squares functions [Manchuk andDeutsch, 2019, Renaudeau et al., 2019] or the dual form of kriging with the use of covariance functions [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part i. modelling method coupling 3d potential-field interpolation and geological rules[END_REF]. Following in the footsteps of [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF], the approach developed in this paper is of geostatistical nature. It exploits and investigates possible improvements of the so-called potential field method.

One drawback of this seminal and efficient method is that it is not suited to the representation of complex structures such as poly-deformed areas, where the available data can be too sparse to correctly describe geometries at the chosen complexity level [START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions[END_REF]. To deal with this issue, other methods typically introduce inferred or interpreted locally varying anisotropies [START_REF] Boisvert | Kriging in the presence of locally varying anisotropy using non-euclidean distances[END_REF]. Alternatively, when facing important uncertainties, experienced users can introduce a priori geological knowledge into the model in the form of additional control points. This workaround is usually time consuming, error prone and user dependent. By its deterministic nature, it also prohibits any uncertainty analysis or quantification.

In the specific case of poly-deformed areas, which are at the core of the present study, the most challenging aspect is to understand, represent and model folded structures [START_REF] Ramsay | The techniques of modern structural geology: Folds and fractures[END_REF]. They are typically complex when formed by a succession of folding events or impacted by salt intrusions. In fact, taking into account these consecutive events can help enhance the implicit modeling framework. For instance, [START_REF] Maxelon | A workflow to facilitate three-dimensional geometrical modelling of complex poly-deformed geological units[END_REF] interpret variations in the orientation of the dominant planar fabrics. [START_REF] Laurent | Implicit modeling of folds and overprinting deformation[END_REF] suggest to integrate additional field data describ-ing folds, lineations and vergence. This approach and the inference of its parameters are discussed in [START_REF] Grose | Structural data constraints for implicit modeling of folds[END_REF]. In another vein, de la Varga and [START_REF] De La Varga | Structural geologic modeling as an inference problem: A bayesian perspective[END_REF] adopt a Bayesian point of view: field data and interpretations are taken as prior information in the model.

While all these geological measurements bring relevant information about fold geometry, one straightforward pointer of the preferential direction of continuity of the fold structure (later referred to as the anisotropy) is the hinge line (Fig. 1). In the implicit approach, this line is directly linked to the underlying scalar field, which stays constant along the hinge line [START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions[END_REF]. It turns out that working with the direction of this line has a practical asset: it allows to use isotropic covariance functions in the interpolation model, while taking into account the anisotropy of the fold structure. The hinge line is deduced in practice from geological observations and more particularly from fold axes data, which are typically readily available in most geological studies. However, they are not always fully exploited.

The main objective of this paper is to explore three possible ways of taking better advantage of these data. First, by adding a second derivative constraint along the hinge line. Under some mild assumptions on the fold structure, it is proven that this is not only an appropriate mathematical translation of the constant property of the scalar field along this line, but also an additional geometrical feature that can help to better define the model. Second, by optimizing the location and number of first order derivative data. Last, by further constraining the model using both first and second derivative data in two orthogonal directions.

We start off in Sec. 2 by summarizing the theoretical background of the potentialfield method. A new implementation of the second-order derivative values of the potential is introduced in this framework. The mathematical advantages of using first and second derivatives are given, and proofs of their ability to improve fold geometry modeling are also provided when needed. Then, a challenging synthetic case study with sparse data is presented in Sec. 3. It is used to illustrate the assets and limitations of using these mathematical objects in Sec. 4, thereby leading to the consideration of the maximum curvature criterion. These results, finally discussed in Sec. 5, should allow geomodeling software developers or users to get more insight on fold axis incorporation for the anisotropy modeling process.

Mathematical Framework

Three-dimensional implicit function construction

The key challenge of three-dimensional structural modeling in geology boils down to a reconstruction problem. It consists in approximating a three-dimensional model of the reality from sparse data sampled on the geological field under study. In the implicit modeling approach, focus is on the scalar field function Z, called potential, that maps any point p ∈ R 3 to a real value

p → Z(p) R 3 → R.
Punctual, localized and sparse measurements of this function are used to estimate it on the entire domain of interest. This is achieved by means of an interpolation technique (co-kriging, radial basis functions). The primary purpose of this scalar field is the extraction of iso-surfaces S α

S α = p ∈ R 3 , Z(p) = α .
They are the surfaces of interest, corresponding to the geological interfaces of the three-dimensional structural model. The so-called potential field method [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF] directly considers increments of the scalar field sampled on iso-surfaces. The implicit function is then defined up to an additive constant. It can be set by choosing a reference iso-surface, the potential value of which is arbitrarily set to zero [START_REF] Calcagno | Geological modelling from field data and geological knowledge: Part i. modelling method coupling 3d potential-field interpolation and geological rules[END_REF].

Typical input data characterization

In order to sample the scalar field, geologists have access to different categories of data. In this work, we focus on field measurement data (or surface data) but the following remarks can be easily generalized to borehole data, or manual or automatic surface interpretation from geophysical surveys. Field measurements may be summarized into two main categories: position point information (p ∈ R 3 ) and orientation vector data (d ∈ R 3 ). With regards to the implicit function Z to build, these field data mathematically correspond to the local evaluation of the function Z(p) or its directional derivative D d Z(p), respectively: i. The assertion that two points, p i and p j , belong to the same interface S α .

They consequently present same potential value, therefore the difference of their potential value (called increment) is null

Z(p i ) -Z(p j ) = 0 ∀ p i , p j ∈ S α .
These points are usually called contact points. Among implicit approaches, increment data constitute the peculiarity of the potential field method; other implicit methods directly integrate contact point information by itself.

Since the potential value of Z(p) at a particular point p is not geological data, increment data are geologically the most convenient to sample, but direct values of potential z i may be introduced alternatively

Z(p i ) = z i ∀ p i ∈ R 3 .
ii. By definition, the first order derivative of Z at point p in a direction d is

D d Z(p) = lim h→0 Z(p + h d) -Z(p) h d .
As the potential values have no physical meaning, there is no way to sample directional derivatives other than identifying direction lines that lie in the plane tangent to geological interfaces. Indeed, as geological interfaces are supposed to be iso-potential surfaces, all derivatives of the potential along such direction lines will vanish. Such directions will usually be sampled as unit vectors τ such that D τ Z(p) = 0, which can also be written

D τ Z(p) = τ.∇Z(p) = ∇Z(p), τ = 0,
where ∇Z(p) denotes the gradient of the scalar field, that writes

∇Z(p) =   D e x Z(p) D e y Z(p) D e z Z(p)   ,
with (e x , e y , e z ) being the standard basis of R 3 . Geologists usually sample the tangent plane of geological interfaces using two such tangent directions, respectively called strike and dip directions. Strike direction, given by the unit vector τ s , is obtained as the intersection of the tangent plane with a horizontal plane that holds the observation point. Dip direction, given by the unit vector τ d , is measured replacing the horizontal plane by a vertical plane.

Last, but not least, geological interfaces can be attributed a polarity, as formations can usually be chronologically ordered. Assuming that the potential Z is a decreasing function of the geological age of the formations, and that ε ∈ {-1, 1} is such that ετ s × τ d points toward younger formations we then know that

(ετ s × τ d ) • ∇Z > 0.
In practical applications and software packages, this inequality constraint is usually replaced by the (biased) assumption that the the gradient of the scalar field has unit norm with the consequence that one can write

∇Z(p) = ετ s × τ d ,
This trick will be denoted as the gradient norm bias in the following.

A new input data: the second-order derivative data

Just as the first-order derivative, the second-order derivative of the potential computed twice in direction d can be taken into account

D 2 d Z(p) = lim h→0 Z(p + h d) + Z(p -h d) -2 Z(p) h 2 d 2 . ( 1
)
Motivations for the addition of such data and its impact on the method are discussed in Sec. 2.5. Let us notice that the directional second-derivative of Z in a direction τ can be expressed using the Hessian matrix of Z, denoted by

D 2 Z(p) in the following D 2 τ Z(p) = t τ.D 2 Z(p)
.τ, where the left-hand side superscript t indicates the transpose of a vector or matrix.

2.4

The interpolation method: the potential field method Using a co-kriging estimator [START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF], the interpolated value of Z at point p, Z * (p), is computed by the weighted sum of different input data

[ Z(p) -Z(p 0 ) ] * = ∑ ii ′ λ ii ′ [Z(p i ) -Z(p i ′ )] + ∑ j λ j ∇Z(p j ), τ j + ∑ k dim ∑ d λ k,d D e d Z(p k ) + ∑ s λ s ( t τ s .D 2 Z(p s ).τ s ), (2)
where the λ ii ′ , λ j , λ k and λ s are coefficients to be determined and Z(p 0 ) is an arbitrary value of Z assigned at the arbitrary point p 0 (most often chosen to be null).

In the universal co-kriging context [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF], the random function, here the potential function Z, is described as a combination of a residual random function with zero mean, Y , and a deterministic drift function, m

Z(p) -Z(p 0 ) = Y (p) + m(p).
(3)

The drift function m is usually expressed as a linear combination of n basis functions

m(p) = n ∑ ℓ=1 ν ℓ f ℓ (p), (4)
where the f ℓ are basis functions and ν ℓ unknown associated coefficients. Kriging equations in this case, detailed in Appendix A, lead to a final co-kriging system K F t F 0

λ µ = K p F p , (5)
where K and K p are the covariance matrices, respectively between data among themselves, and data and interpolated points; F and F p are the universality conditions on data and interpolated points and λ and µ, the unknown variables to solve.

Motivations for second-order derivative data implementation

The main difference here, compared to the initial formulation of the method from [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF], is the addition of second-order derivative values of potential in an established direction into the co-kriging equations. Original observations, which lead to the motivations of this development, are detailed below, as well as the mathematical proof of its merits.

Empirical observations

This new development has been motivated by the observation that, in case of folds with planar hinge line, in addition to the first-order derivative values of the potential (tangent), the second-order derivative values of the potential for a point on the hinge line are (Fig. 2):

Pt Sa S 1 S 0 L 3 L 2 L 1 n H 2 H 1 DA H 3 n'
H in ge li n e i. null in the direction of the hinge line, ii. maximal in the orthogonal direction (on the tangent plane).

To illustrate this, values of different directions for second-order derivative of an arbitrary point of a hinge line (crest here) of the synthetic case (Sec. 3) are shown on Fig. 2. The same pattern can be observed for all points on the crest, with a different value for the direction angle φ .

Demonstration and application scope

Concerning the first observation, this empirical observation is supported by a theoretical demonstration developed in Appendix C for particular folds (Fig. 1). The only requirement is that the axial poles n along the hinge line are colinear, or equivalently that there exists a plane to which they are all normal, defined here as the 'tangent plane' Pt.

Here the study focuses on an upright fold, with vertical axial surface Sa. In this case, the hinge line of the studied iso-potential surface coincides with the crest of the fold. However, the demonstration that the second-order directional derivative is null along a particular curve only relies on the existence of the tangent plane along this curve. This means that the conclusions are the same for inclined folds, as long as the axial poles n are still colinear. In the case of cylindrical folds, in particular, the second-order directional derivative is null along any line parallel to the hinge line.

The demonstration provided in Appendix C shows that the property stands for both syncline and anticline folds. Moreover, the axial surface and hinge line may freely vary in direction in the tangent plane and hence characterize anisotropy of the fold, irrespective of its shape. This opens up new constraints for fold axis data, restricted in the present work to particular types of poly-deformed folds, but that could conceivably be extended to more complex structures: even in cases where the hinge line is not a planar curve, if a portion satisfies the tangent plane requirement, we can incorporate the null second-order directional derivative information along it at that location.

The second-order derivative of a three-dimensional function gives information about the local behavior of this function, and thus intrinsically, for implicit approaches, on resulting isosurface curvature. This could be a significant information, especially in case of very curvy hinge lines. Indeed, the geometric meaning of these two previous observations is then: 1. a null value of second-order derivative of potential in the hinge line direction, implies local flatness in the hinge direction 2. the maximal value of directional second-order derivative of potential in the orthogonal direction to the hinge line on the tangent plane means that the curvature of the resulting isosurface is maximal orthogonal to the hinge direction.

Interest of second-order derivative data for assessing local geometry

Null first-derivative values of potential in all directions on the tangent plane already express local flatness in the hinge direction. However, the addition of a null secondorder derivative value provides an additional information of 'smoothness'. Indeed, assuming the potential Z is sufficiently smooth, its Taylor expansion at a point p on the hinge can be written

Z(p + h) = Z(p) + ∇Z(p).h + h t D 2 Z(p)h + O( h 3 ).
For h in the tangent plane, the first-order term ∇Z(p).h is zero. If the hypotheses of Sec. 2.5.1 are satisfied and h is along the hinge, the second-order term h t D 2 Z(p)h is null as well, so that Z(p + h) -Z(p) is a third-order infinitesimal; by providing the information D 2 τ Z(p) = 0 for τ along the hinge line, one can expect the resulting model to be smoother and then better constrained.

Fold axis data sampling

These mathematical properties are especially relevant for already existing field data, such as fold axis data. This data consists in direction measurements (azimuth and plunge values). In this section, we consider fold axis data sampled on the field. These data are often largely rather interpreted than measured on geological maps but rarely integrated, however, into final models. They represent unit localized vectors and do not necessarily belong to the hinge line. Indeed, they may belong to the axis of another sub-parallel surface of the folded system, but this information is still a relevant constraint about the folded structure for the scalar field.

This direction measurement may be sampled in different manners (Fig. 1):

1. Fold axis may be measured directly on the axial surface (D A and D A ′ ), on one hinge line. This hinge zone is a fragility zone and is not often well preserved, these observations are then pretty rare, and thus valuable. 2. In fact, fold axis data come abundantly from intersection lineations measurements (D B ). These lines are the result of the intersection between bedding planes (S 0 ) and cleavage planes (S 1 ), which commonly takes place in folding processes. 3. Fold axis data may also be inferred from orientation measurements of fold limb normal vectors (gradient data), called pole vectors. Then, stereographic projection tools can be used to determine local fold axis data (D C measurement from n and n ′ pole vectors).

The first measurement is directly positioned on one hinge line. In practice, even though the latter are measures at other locations, they are commonly replaced at the hinge line location. Then, this paper studies the use of these three raw measurements on the hinge line.

However, wherever they come from, these fold axis data measurements significantly contribute to the global hinge line interpretation over the domain. This global interpretation is performed by geologists from field data but also other observations with various interpretative values, such as micro-scale structures, parasitic folds, global regional context and expert knowledge. For a complete review of fold data see [START_REF] Hudleston | Information from folds: a review[END_REF], [START_REF] Mcclay | The mapping of geological structures[END_REF]. 7). The elliptic hinge line of the fold, with semi-major and semi-minor axes, A and B, is presented on (a) (dashed line). Random data for different samplings are generated with the angle θ around this structure. Data sampling zones, A1 and A2, with their respective cut off angles,

c 1 = 2π 24 , c 2 = 4π 24 , c ′ 1 = 8π 24 , c ′ 2 = 10π
24 , are also shown, as well as the vertex of the ellipse.

Covariance function requirements

Implementation of second-order derivative value of potential leads to the need of a four-times differentiable covariance function (see Appendix B.2 for more details on successive derivations of the covariance function). Thus, the classically implemented cubic covariance model is not applicable; here the considered model is one of the 'piecewise polynomial, positive definite and compactly supported radial functions of minimal degree' introduced by [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF], denoted by Wendland C4 in the following, and which may be expressed as

C(r) = 1 3 (1 -r) 6 + (35r 2 + 18r + 3). ( 6
)
3 Case study

Presentation of the considered potential field

In this paper, an explicit theoretical potential function has been chosen as case study. Shaped to represent an overturned anticline structure, the hinge line -or crest, since they coincide here -describes an ellipse trajectory of semi-major and semi-minor axes A and B (Fig. 3). Therefore, this trajectory represents the main continuity direction of the structure, defined as anisotropy, here significantly varying in space. The explicit potential function can be expressed as

Z(x, y, z) = z -α sin(ω d) if d ∈ [l -ε ; l + ε], z otherwise, (7)
where

-d = a (x -y) 2 + b (x + y) 2 ,
α = 1.9, ω = 2, a = 0.5 and b = 0.05, with a = 1 B 2 and b = 1 A 2 , l and ε being chosen so that sin(ω d) = 0 for d = lε and d = l + ε, leading to continuous isosurfaces: ε = π ω and l = 5π 2ω . Only one iso-potential surface of this field at level Z(x, y, z) = 2.5 is shown on Fig. 3: this surface is defined as our reference of study thereafter.

Initial results on different samplings

In the implicit modeling framework, an anisotropic structure like this one requires a dense sampling, in terms of quantity and localization of data, in order to obtain an admissible model. The following section will highlight interpolation issues encountered in such a sparse data environment.

Random sampling strategy

In basic geological structural studies, sampled data are usually contact points and gradient data. In that respect, considering the initial potential field method, Fig. 4 presents two different resulting interpolated surfaces, based on two different initial samplings. These samplings have been randomly generated in different zones of the structure (denoted A1 and A2 zones on Fig. 3). Resulting models of the potential field method on each sampling is presented in Fig. 4.

-Scenario 1 (Fig. 4a): The first sampling consists in 4 random contact points (with 4 associated gradient data) generated in the A1 zone. -Scenario 2 (Fig. 4b): The second sampling is defined by 4 random contact points (with 4 associated gradient data) generated in the A1 zone, and 4 random contact points/gradient data generated in the A2 zone.

As a consequence, both sampling density and area coverage of data are higher in Sampling 2 compared to Sampling 1.

Interpolation parameters

In the following, the implemented covariance function is the Wendland C4 (Eq. 6), the range of the model is equal to the diagonal of the model ( 15) and a vertical drift is considered here

m(p) = ν z(p).
The interpolation is performed on a 20-by-20-by-20 grid (although for numerical reasons the error score V e (Sec. 3.2.3) is computed on a 20-by-20-by-80 grid). Although true values of gradient data are known in this synthetic example, gradients are defined unitary as a default choice, as it is implemented that way in software packages (Sec. 2.2).

Modeling quality indicators

This work focuses mainly on the qualitative aspect of the resulting three dimensional models, as it is generally performed in geological studies. For this purpose, the threedimensional isosurface of interest is presented for each experiment thereafter. In addition, in order to help visualization, the resulting surface of interpolation is compared to the original theoretical surface using a vertical cross section along the 'first diagonal' (x = y, profile A-B on Fig. 4a). This cross section is shown underneath each corresponding three dimensional model, even though the cross section direction is not always displayed. However, a quantitative measure, the volume error V e , is provided. This indicator estimates the volume difference between the theoretical and modeled surfaces in the region near the structure of interest; its computation is detailed in Appendix D. This measure, made possible here by the relative simplicity of the synthetic model, can be distorted by numerical errors, notably due to the grid resolution in z. An uncertainty value ∆V is computed and provided for each model, and has to be taken into account when comparing values of the indicator V e . In the end, the investigations performed in this work can be seen as attempts to minimize this error score. These different scores are reported in appendix D. To express improvements of the resulting model compared to initial results in the same scenario i (with its own volume error V e i ), the relative volumic improvement (RVI) ηV e is introduced

ηV e = V e i -V e V e i
where V e i is the volume error of scenario i; i ∈ (1, 2)

This signed indicator ηV e quantifies the reduction of error compared to the initial model It is expressed as a percentage. Negative values mean deterioration compared to initial results, 100% RVI means that the error is now null, whereas 0% RVI means no change in volume error value. More details are provided in Appendix D. These criteria, both qualitative and quantitative, have to be handled simultaneously in order to evaluate the results and then provide better insight on the scalar field variations. All resulting values for the different experiments are summarized in appendix D.4.

Critical case versus well-constrained scenario

Resulting surfaces of the potential field method with the two different samplings of Sec. 3.2.1 are shown in Fig. 4. This interpolation is performed with the choice of parametrization detailed in Sec. 3.2.2. These two starting assumptions lead to two different results:

The resulting surface with Sampling 1 (Fig. 4a) presents a 'bubbly' aspect far from the data, as it is generally observed with implicit modeling approaches in sparse data situations. With no data around, the potential value tends towards the drift (mean) value. This result is actually not satisfactory for geologists. The volumetric error is V e 1 = 104.68. Moreover, this phenomena takes place here with use of Wendland C4 covariance model, but also occurs with the more usual cubic covariance one. This example may be considered as a critical case scenario, where fold anisotropy introduction would significantly benefit the model. Proposed solutions will be evaluated in comparison with this vertical error score, which should never exceed this critical value.

On the other hand, the resulting surface obtained with the better sampling, Sampling 2, is shown in Fig. 4b. This sampling choice leads to a better result, as could be expected. The volumic error is V e 2 = 48.36, which is 56.32 better than the previous scenario. This change of scenario improves the result by 54%. Hence this scenario may be considered as a better constrained model. Moreover, despite this prevailing error, this result may be qualitatively considered as satisfying and thus, this scenario considered as a well constrained model. This can be defined as a benchmark for further comparisons with improvements made on the previous case. Still, the impact of additional anisotropy information data will also be tested on this model in order to see the impact of these new data on a satisfying well constrained model, but where there is still possible room for improvement.

These two sampling starting assumptions are mentioned as Scenario 1 and Scenario 2 in the following.

Exploratory methodology

The first scenario is not satisfying, especially considering the fact that in such situations we observe that geologists actually already do hold anisotropy information: either localized fold axis data or a final interpretation of the hinge line (Sec. 2.5.4). The parsimonious deployment of both kinds of inputs is discussed in Sec. 5.

However, the following results on the exploitation of fold axis data fit for both approaches. Indeed, the integration of the hinge line may be done by discretizing it and then use punctual fold axis data, which is the correct mathematical expression of this information. Nevertheless these conclusions raise different outstanding issues: i. Whereas in theory null second-order derivatives yield additional smoothness information compared to first order derivative values, in practice for hinge line conditioning, is it a real relevant constraint to add? ii. Although the use of tangent data was already possible in the initial method of [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF], in practice, the use of such constraints for fold geometry modeling is not necessarily always adopted by geologists, even though data of different types mentioned in Sec. 2.5.4 are available. When tangent data are taken into account, there is still uncertainty about the necessary amount of them to well constrain poly-deformed folds in sparse data environment, as well as a geographical distribution preference. These understandings would allow to define priorities for either the direct geological sampling strategy of fold axis data or the parametrization of an interpreted hinge line. iii. Considering the tangent plane, for a given data field measurement (fold axis) on the hinge line, different geometric constraints may be taken into account: at least two tangent data, a null second-order derivative value and a maximal secondorder derivative value. The relevance of each one should be investigated.

Investigations on better fold axis data integration

In the following, these three issues are addressed on the synthetic case regarding the two previous scenarios proposed in Sec. 3. The hinge line direction is known; different constraints are imposed upon it depending on the experiment. All the resulting surfaces are compared to the theoretical one (Fig. 3) through cross sections just as presented in Sec. 3.2, even though its orientation is no longer indicated. The different volumetric error indicators are listed in appendix D.4.

4.1 Use of fold axis data as tangent data, null second-order derivative data or both On the hinge line, a null first-order derivative value of potential constraint may be imposed for fold axis vector data. New observations made in Sec. 2.5.1 motivated the addition of a null second-order derivative value of potential constraint as well.

In order to characterize the contribution of each constraint in practice, Fig. 5 shows the impact of such data on the synthetic case, regarding the two different starting scenarios. 15 homogeneously distributed vectors are sampled along the hinge line 5h). For comparison, the initial results from Sec. 3.2.4 are also recalled (Figures 5a,5b). The conclusions are the following:

-Scenario 1: The introduction of null second-order derivative constraints does have an impact on the resulting model (Fig. 5e), in comparison with the initial result of this scenario (Fig. 5a). These constraints appear to flatten the unwanted bulging as expected; the resulting surface draws near the theoretical one. V e is 92 ± 4, RVI is around 13%. However, the introduction of tangent data constraints alone (Fig. 5c) visually leads to a better result: the modeled hinge line seems to be localized more accurately. In addition, V e is 85 ± 4, RVE is around 19%. Finally, the combination of both constraints (Figures 5g) leads to a qualitatively similar result that the one with tangent data constraints only. However, V e is 91 ± 4, RVI is around 13%, which is closer to the result with null second order derivatives. The differences between these three results fall within the uncertainty bounds given by ∆V : direct comparisons should be qualitative in nature. -Scenario 2: In this specific well constrained case, the RVI scores of these four comparisons indicate that the four models are equivalent (within the uncertainty range, at this resolution). Qualitatively, these models do not show perceptible differences either. As a conclusion, in this well constrained scenario, the addition of either tangent or null second-order derivative data does not seem necessary.

In the end, these different results do not bring to light the interest of the introduction of second-order derivative value made on Sec. 2.5.3. This could possibly be explained by the pre-existence of local flatness, or at least smoothness, of the model preceding its introduction. Actually this smoothness is generally inherent in geological surfaces. As a matter of fact, this local feature is the main assumption of geological modeling interpolation methods. However, second-order derivative values could bring relevant information in ambiguous surfaces in other surface modeling applications.

Quantity or spatial distribution of tangent data impact

In response to observations made on Sec. 4.1, it still appears relevant to sample fold axis data and incorporate them in terms of tangent data alone at least. The impact of the modalities of this integration, such as quantity and spatial distribution remain ill-known. For this purpose, the impact of quantity and spatial distribution of these data has been studied on the synthetic example.

Impact of the quantity of tangent data

The impact of the quantity of tangent data along the hinge line is described in Fig. 6. 29 initial tangent data regularly spaced along the hinge line have been decreased to 8, then 4 data points only.

-Scenario 1: The results show that the quality of the resulting model remains the same for these three numbers of tangent data values. The relative improvement (RVI) is about 19% for 29 tangent data (Fig. 6a), 8 data (Fig. 6c) or 4 data (Fig. 6e). This principally suggests that the integration of a limited number of tangent data already shows a positive impact on the result, thus a fine distribution of these tangent data is not necessary, at least in this particular case. -Scenario 2: In the well-constrained scenario, the number of tangent data does not seem to have a significant impact on the result. These different constraints lead to similar results (Fig. 6b, Fig. 6d, Fig. 6f) on the edge of the indicator's resolution. These statements are also similar to the ones observed in the previous experiment with derivative constraints (Sec. 4.1).

Impact of the spatial distribution of tangent data

Now considering a fixed number of tangent data, the impact of the localization of such data along the hinge line is described by Fig. 7. Three major sampling strategies have been adopted here: new tangent data regularly spaced along the hinge line, new data close to the initial inputs (contact points and gradient data) and new data far from the initial inputs. Adding anisotropy information in the last location could potentially be of some specific interest, since this zone is the ellipse's vertex area (Fig. 3) where the anisotropy direction varies strongly, in comparison to the whole structure.

-Scenario 1: With a constant number of tangent data, the results show that the use of these data close to initial contact points and gradient data (Fig. 7c), or regularly spaced (Fig. 7a), qualitatively give the best results, by flattening the budging around the hinge line. The quantitative improvement of these strategies is around 11% and 15%, respectively. However, regularly spaced data visually seem to improve the result mainly around the vertex area, whereas tangent data close to initial data improve the result around these well-constrained zones. This highlights that local information of flatness is given. On the contrary, additional tangent data placed far from contact points and gradient data (Fig. 7e) make the modeling worse by -14%. This could be explained by the contradiction between the forced local flatness given through six tangent data constraints at the vertex position here and the forced unitary norm given by gradient data which may lead to some unexpected results (Sec. ii.). -Scenario 2: The addition of tangent data does not significantly affect the resulting model whatever position of these data is (Figures 7b,7d,7f). RVI values are similar, within the limit of their uncertainties.

As a conclusion, with respect to this only case study, it seems that a parsimonious regularly spaced sampling of tangent data is the most cautious strategy to adopt. Adding anisotropy information where it is particularly needed could, counterintuitively, worsen the modeling, possibly because of the gradient norm bias.

Insights on orthogonal constraints of fold axis data on the tangent plane

The previous sections, 4.1 and 4.2, have considered only fold axis data directions, i.e. the hinge line direction. However, it is possible to take into account the orthogonal direction on the local tangent plane (Fig. 1) of these fold axis data. This consideration would be relevant for two major reasons. First, the first-order derivative value is null in this direction. In practice this second tangent data is sometimes already integrated in the model, completed with a younging direction, which indicates the polarity of the normal of the tangent plane. This is performed in order to express pole vectors (Fig. 1) on the hinge line. Indeed, this vector shall not be provided to the model as a gradient input, because of the gradient norm bias. In the end, more insights on this tangent constraint impact would benefit to the understanding of its relevance. Second, the observations made in Sec. 2.5.1 revealed that this direction is a maximum value of the second-order derivative. This new constraint could have potential interest in search of the valuation of fold axis data in the enhancement of geometric control of fold structures. For this purpose, Fig. 8 shows the consecutive impact of these four constraints, first and second-order derivatives in the two mentioned directions, for 6 fold axis data on the crest line:

-Scenario 1: The use of tangent data in the hinge line direction only (Fig. 8a) leads to the same observations previously made in this case (Sec. 4.1,4.2). However, the addition of the second tangent data (Fig. 8c) seems to conduct to a flatter hinge line, but not especially better positioned. The addition of null-second order derivative in the hinge line direction (Fig. 8e) does not add an appreciable qualitative or quantitative difference. Finally, the incorporation of second-order derivative values in the orthogonal direction (Fig. 8g) leads to noticeably better results; the resulting surface appears to be closer to the theoretical one than it has ever been so far, and the RVI score is around 68%. The volumic error V e is even better than the initial result in the well-constrained scenario S2 (33.34 instead of 48.36). -Scenario 2: In this well constrained scenario, there is no significant difference between the first three constraint influences, i.e. the addition of two tangent data (Figures 8b,8d) and null second-order derivative value (Fig. 8f). Just as in the case of Scenario 1, the final addition of second-order derivative values in the orthogonal direction (Fig. 8h) conducts to a very satisfying model. Qualitatively, the resulting surface is very close to the theoretical one, closest than it has been at any time for the different attempts in this work. RVI score is around 36%, whereas the previous results with this same sampling did not show positive improvement scores.

The observations performed on this single synthetic case lead to several conclusions. In addition to tangent data in the hinge line direction, the integration of orthogonal tangent data does not appear in this case to add sufficient information to constrain the fold structure as expected. The conclusions on the addition of null second-order derivative data are the same as the ones previously made (Sec. 4.1: these data do not affect significantly the model. However, the addition of the second-order derivative values in the direction where they are maximal, orthogonal to the hinge line, seems to affect in a very significant way the results; at least as it could be expected by users in such situation.

Unfortunately, in practice it is not possible to directly sample these values on the field. Then, the immediate integration of this constraint on real cases for now would only be possible by an empirical approach. Indeed, the value of this constraint controls the orthogonal curvature, hence the local 'pinch' of the fold. This value could potentially be estimated by a trial-and-error approach with visualization tools, for example. However, these observations are worth exploring, and should open a path for further mathematical developments on the integration of a local maximum as a criterion input for implicit methods.

Conclusion and discussion

This paper tackles a well-known issue encountered in the potential field method: when the data is sparse, classical models can lead to unsatisfactory results (e.g. presenting a 'bubbly' aspect). This is particularly the case when the geological field is made of complex folded structures. In such situations, geologists are forced to add artificial control points. The alternative proposed here is to incorporate information about the anisotropy. It can be characterized by fold axis data, which have the advantage of being readily available in most fold structural studies. This approach is motivated from a mathematical point of view in Sec. 2, using the potential method [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF] as a reference. This allows to describe poly-phased structures in a purely geometric, parsimonious manner, without having to make strong assumptions on the underlying succession of geological events. Another asset of this method is that the covariance function used for the interpolation of the potential scalar field can remain isotropic.

A challenging case study is described in Sec. 3. It serves in Sec. 4 as a founding example to gain more insight into the usage of fold axis data for characterizing the anisotropic geometries of folds. Three major conclusions are drawn from this exploration.

i. From a mathematical point of view, the second-order derivative values of the three-dimensional scalar field are shown to be of particular interest. They bring theoretical information on the curvature of the surface, notably where it is flat and might coincide with the hinge line. However, in the chosen example, adding null second-order derivative constraints to the model only slightly improves the results. When combined with tangent data, almost no improvement is observed. Obviously, no generalization can be drawn from a single example. Nonetheless, it might be conjectured that the effect of these mathematical objects are bound to be limited in other geological studies too. Indeed, one of the most important features characterizing complex iso-surfaces is their smoothness. It is at the foundation of geological interpolation methods, which try to minimize the curvature. When dealing with fold structures, curvatures turn out to be quite regular, especially in the hinge line direction. Therefore, it is our belief that tangent data suffice to capture this key characteristic. ii. Tangent data turn out to have a significant impact on the quality of the results; when ignored, undesirable effects are observed. Even with a very small number of such data, fold anisotropy integration is notably improved. It is also observed that their effect is all the more positive as fold axis vectors are sampled close to the initial gradient and contact point data. iii. The value of the second-order derivative in the orthogonal direction of the hinge line (i.e. the maximum value of second-order derivatives) appears to be particularly useful. Models incorporating this information in a synthetic example show promising results. Unfortunately, this value is hardly reachable in practice, although it could conceivably be estimated if a cross-section is available. Alternately, an arbitrary value could be chosen, which would control the local pinch of the fold. Such observations suggest to investigate the incorporation of a Gaussian curvature condition in the model; this is left for future research.

All these remarks hold for discrete fold axis data. They usually come from raw measurements, but could also be generated from the discretization of an interpreted hinge line. Indeed, the ontological boundary between these two types of model input is not immediate: the interpretation of the hinge line stems from a combination of fold axis data and other geological observations such as micro-scale structures or the regional context. In practice, the interpreted hinge line could be directly integrated as a sketch (expert interpretation) or automatically interpolated from some direction measurements (automatic interpretation). Then, it could be discretized with synthetic fold axis data. A non-negligible limitation of this technique is that the uncertainty linked to such artificial input is generally intractable. Nonetheless, interpreted data hold as much geological value as field measurements [START_REF] Frodeman | Geological reasoning: Geology as an interpretive and historical science[END_REF]. Therefore, favoring either one over the other should remain an expert choice, based on practical considerations such as the desired model complexity, the state of play, or the targeted applications. It goes beyond the methodological scope proposed here.

In any case, using a discretized version of the interpreted hinge line remains more desirable, in our opinion, than introducing purely arbitrary control points (e.g. contact points or gradient data). Indeed, the latter generate complex, non-parsimonious models of the anisotropy. In particular, contact points force the surface to pass through a particular point; this condition can be viewed as too restrictive. In addition, gradient data control the contraction of the scalar field (isosurfaces get closer to each other with the gradient norm) and thus introduce the gradient norm bias. The interpreted hinge line, on the other hand, gives some degree of freedom to the model and solves the gradient norm bias by adding tangent data. By focusing on the two leading types of information that characterize the anisotropy, namely the null and maximum curvatures, it produces a parsimonious realistic model. The number of tangent or second-derivative data to add to the model must be pondered relative to the induced computational complexity, since each constraint adds one line to the co-kriging system.

Another way of incorporating the hinge line into the model, which was not considered here, would be to use it as a drift function. This would require an explicit analytic formula describing the hinge line, instead of a collection of discrete values. In practice, it is possible to create an explicit fold potential from an initial threedimensional curve sketch using the distance to a skeleton [START_REF] Henrion | Odsim: an object-distance simulation method for conditioning complex natural structures[END_REF][START_REF] Garcia | Automatic generation of geological stories from a single sketch[END_REF]. This approach brings geological modeling and computer graphics techniques closer together, but quantifying the uncertainty of such data remains an open question. This is left for future research.

Acknowledgements The authors would like to thank Gabriel Courrioux of BRGM for being a pioneer of the geomodeler software and for initiating this research work with his colleague, Bernard Bourgine. We are particularly grateful to Christian Lajaunie, inventor of the geostatistical potential field method, who largely participated in the development of the consecutive derivatives of the covariance function. Support and advice from the Geostatistics Research Group, especially Emilie Chautru for her careful editing, were also helpful, just as the ones from the Geomodeling Team of BRGM. Authors thank Michael Hillier, Guillaume Caumon and the anonymous reviewer for their fruitful comments.

A Formulation of the potential field method in terms of kriging equations

Whereas the seminal paper by [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: principles of a method based on potential interpolation[END_REF]] adopted a dual (co-)kriging formulation, we retained here the universal (co-)kriging formulation.

The interested reader may consult the chapter 3.4 of the book by [START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF] to have more details about this derivation.

Equation 2 writes:

[ Z(p) -Z(p 0 ) ] * = ii ′ n ∑ ii ′ λ ii ′ [ Z(p i ) -Z(p i ′ ) ] + gn ∑ g λ g Z ′ g (p g ) + sn ∑ s λ s Z ′′ s (p s ) (A.8)
where the estimate [ Z(p) -Z(p 0 ) ] * is a linear combination of:

ii ′ n increment data, g n first-order directional derivative of potential data (tangent and/or components of gradient data), s n second-order directional derivative of potential data, and where for the sake of simplicity, we also index the directional derivatives Z' and Z" with g and s rather than with directions τ g and τ s :

Z ′ g (p g ) = D τg Z(p g ), Z ′′ s (p s ) = D 2 τs Z(p s ).
(A.9)

Similarly, in the following, in summations of first and second-order directional derivatives of any function, we will drop the directions and index by g and s as well.

The estimation error ε is defined by

ε = [ Z(p) -Z(p 0 ) ] * -[ Z(p) -Z(p 0 ) ] = ii ′ n ∑ ii ′ λ ii ′ [ Z(p i ) -Z(p i ′ ) ] + gn ∑ g λ g Z ′ g (p g ) + sn ∑ s λ s Z ′′ s (p s ) -[ Z(p) -Z(p 0 ) ].
(A.10)

A.1 Universality conditions

Universality conditions ensure that the estimation is unbiased in the sense that the expectation of the estimation error is null i.e. E(ε) = 0.

Noting that

E[Z(p) -Z(p 0 )] = m(p) = n ∑ ℓ ν ℓ f ℓ (p), E[Z ′ g (p g )] = m ′ g (p g ) = n ∑ ℓ ν ℓ f ℓ g ′ (p g ), E[Z ′′ s (p s )] = m ′′ s (p s ) = n ∑ ℓ ν ℓ f ℓ s ′′ (p s ), (A.11)
with the same abbreviation logic as previously and m

′ g (p g ) (resp. f ℓ g ′ (p g )) standing for D τg m(p g ) (resp. D τg f ℓ (p g )).
Then, universality conditions write:

E[ε] = E[ ii ′ n ∑ ii ′ λ ii ′ [ Z(p i ) -Z(p i ′ ) ] + gn ∑ g λ g Z ′ g (p g ) + sn ∑ s λ s Z ′′ s (p s ) -[ Z(p) -Z(p 0 ) ] ] = ii ′ n ∑ ii ′ λ ii ′ [ E[Z(p i )] -E[Z(p i ′ )] ] + gn ∑ g λ g E[Z ′ g (p g )] + sn ∑ s λ s E[Z ′′ s (p s )] -E[Z(p) -Z(p 0 )] = ii ′ n ∑ ii ′ λ ii ′ [ n ∑ ℓ ν ℓ f ℓ (p i ) - n ∑ ℓ ν ℓ f ℓ (p i ′ ) ] + gn ∑ g λ g n ∑ ℓ ν ℓ f ℓ g ′ (p g ) + sn ∑ s λ s n ∑ ℓ ν ℓ f ℓ s ′′ (p s ) - n ∑ ℓ ν ℓ f ℓ (p) = n ∑ ℓ ν ℓ ( ii ′ n ∑ ii ′ λ ii ′ [ f ℓ (p i ) -f ℓ (p i ′ )] + gn ∑ g λ g f ℓ g ′ (p g ) + sn ∑ s λ s f ℓ s ′′ (p s ) -f ℓ (p) )
= 0, (A.12) which must be true for all values of ν ℓ and implies n equations (universality constraints), as many as there are drift functions:

ii ′ n ∑ ii ′ λ ii ′ [ f ℓ (p i ) -f l (p i ′ )] + gn ∑ g λ g f ℓ g ′ (p g ) + sn ∑ s λ s f ℓ s ′′ (p s ) -f ℓ (p) = 0 ∀ℓ ∈ {1, . . . , n} (A.13) denoted T ℓ (λ ii ′ , λ g , λ s ) = 0 in the following.

A.2 Optimality conditions

Optimality conditions ensure that Var(ε) is minimal.

In this part, the following convention is adopted:

Cov(Z(p i ), Z(p j )) = C •• i, j , Cov(Z ′ g (p g ), Z(p j )) = C • ′ • g, j , Cov(Z ′′ s (p s ), Z ′ j (p j )) = C • ′′ • ′ s, j . (A.14)
When double indices are used, the notation will refer to increment values, for example C •• ii ′ , j will refer to Cov(Z(p i ) -Z(p i ′ ), Z(p j )).

The variance of the estimation error (A.10) is

Var(ε) = Var ii ′ n ∑ ii ′ λ ii ′ [ Z(p i ) -Z(p i ′ ) ] + gn ∑ g λ g Z ′ g (p g ) + sn ∑ s λ s Z ′′ s (p s ) -Z(p) = ∑ ii ′ , i 2 i 3 λ ii ′ λ i 2 i 3 C •• ii ′ , i 2 i 3 + ∑ g,g ′ λ g λ g ′ C • ′ • ′ g,g ′ + ∑ s,s ′ λ s λ s ′ C • ′′ • ′′ s,s ′ -Var[Z(p)] + 2 ( ∑ ii ′ ,g λ ii ′ λ g C •• ′ ii ′ ,g + ∑ ii ′ ,s λ ii ′ λ s C •• ′′ ii ′ ,s + ∑ g,s λ g λ s C • ′ • ′′ g,s -∑ ii ′ ,p λ ii ′ C •• ii ′ ,p -∑ g,p λ g C • ′ • g,p -∑ s,p λ s C • ′′ • s,p ) (A.15)
Lagrange multipliers are used to minimize the quadratic form Var[ε(λ ii ′ , λ g , λ s )] under the constraint that all universality conditions hold i.e. T ℓ (λ ii ′ , λ g , λ s ) = 0 for all ℓ.

This implies that Var[ε(λ .16) should have its partial derivatives with respect to each λ ii ′ , λ g , λ s and ℓ Lagrangian multipliers µ ℓ to be nul.

ii ′ , λ g , λ s )] + 2 ∑ l µ l T ℓ (λ ii ′ , λ g , λ s ) (A
The Lagrangian L is then given by .17) Derivating this Lagrangian according to λ ii ′ , λ g , λ s and µ l gives .18) which defines kriging equations and writes in matrix form:

L(λ ii ′ , λ g , λ s , µ l ) = ∑ ii ′ ,i 2 i 3 λ ii ′ λ i 2 i 3 C •• ii ′ ,i 2 i 3 + ∑ g,g ′ λ g λ g ′ C • ′ • ′ g,g ′ + ∑ s,s ′ λ s λ s ′ C • ′′ • ′′ s,s ′ + 2 ( ∑ ii ′ ,g λ ii ′ λ g C •• ′ ii ′ ,g + ∑ ii ′ ,s λ ii ′ λ s C •• ′′ ii ′ ,s + ∑ g,s λ g λ s C • ′ • ′′ g,s ) -2 ( ∑ ii ′ λ ii ′ C •• p,ii ′ + ∑ g λ g C •• ′ p,g + ∑ s λ s C •• ′′ p,s ) -2 ∑ ℓ µ ℓ T ℓ (λ ii ′ , λ g , λ s ) (A
                             ∑ i 2 i 3 λ i 2 i 3 C •• ii ′ ,i 2 i 3 + ∑ g λ g C •• ′ ii ′ ,g + ∑ s λ s C •• ′′ ii ′ ,s + ∑ ℓ µ ℓ [ f ℓ (p i ) -f ℓ (p i ′ )] = C •• p,ii ′ ∀ (i, i ′ ) ∑ g ′ λ g ′ C • ′ • ′ g,g ′ + ∑ ii ′ λ ii ′ C •• ′ ii ′ ,g + ∑ s λ s C • ′ • ′′ g,s + ∑ ℓ µ ℓ f ℓ g ′ (p g ) = C •• ′ p,g ∀ g ∑ s ′ λ s ′ C • ′′ • ′′ s,s ′ + ∑ ii ′ λ ii ′ C •• ′′ ii ′ ,s + ∑ g λ g C • ′ • ′′ g,s + ∑ ℓ µ ℓ f ℓ s ′′ (p s ) = C •• ′′ p,s ∀ s ∑ ii ′ λ ii ′ [ f ℓ (p i ) -f ℓ (p i ′ )] + ∑ g λ g f ℓ g ′ (p g ) + ∑ s λ s f ℓ s ′′ (p s ) = f ℓ (p) ∀ ℓ (A
K F t F 0 λ µ = K p F p (A.19)
Denoting N = (ii ′ n + g n + s n ), the number of data, and p, the number of points to estimate : .20) with the same abbreviation logic f ℓ g

K is N × N, K p is N × p, F is N × n, F p is n × p and 0 is a block zero matrix. K =    K •• ii ′ ,i 2 i 3 K •• ′ ii ′ ,g K •• ′′ ii ′ ,s K • ′ • g,ii ′ K • ′ • ′ g,g ′ K • ′ • ′′ g,s K •• ′′ s,ii ′ K • ′′ • ′ s,g K • ′′ • ′′ s,s ′    K p =    K •• ii ′ ,p K • ′ • g,p K • ′′ • s,p    F p =   f 1 (p 1 ) • • • f 1 (p p ) • • • f n (p 1 ) • • • f n (p p )   F =   f 1 (p i ) -f 1 (p i ′ ) • • • f n (p i ) -f n (p i ′ ) f 1 g ′ (p g ) • • • f n g ′ (p g ) f 1 s ′′ (p s ) • • • f n s ′′ (p s )   (A
′ (p g ) = D τg f ℓ (p g ) and f ℓ s ′′ (p s ) = D 2 τs f ℓ (p s ).

B Covariance function derivation developments

B.1 Definitions

Covariance function of Z between point p and point q

(p, q) → Cov(Z(p), Z(q)) R 3 × R 3 → R. (B.
r(h pq ) = h pq A = (p -q) t A (p -q) R 3 → R the anisotropic norm, (B.26)
where A is a fixed rotation/dilatation matrix of R 3 .

In the following, r(h pq ) may be denoted r for readability.

B.2 Notations

B.2.1 Differentiation with respect to variable p or q

Although by assumption the covariance function of Z between points p and q only depends on the anisotropic distance r between them, in order to compute the partial derivatives needed in our model, we still have to consider it as a function of two variables p and q in R 3 . In the following, we will need to differentiate with respect to either of them, which we will denote with a superscript indicating the variable we differentiate upon.

For instance, for a function F : R 3 × R 3 → R, we note the partial derivative with respect to p along the vector τ as D p τ F(p, q) = lim h→0 1 h (F(p + hτ, q) -F(p, q)) .

Denoting (e x , e y , e z ) the standard basis of R 3 , we can then define the gradient with respect to p as

∇ p F(p, q) =   D p ex F(p, q) D p ey F(p, q) D p ez F(p, q)   .
By differentiating this with respect to q, we can define D p,q F(p, q) and by further differentiating with respect to one or the other variable, we can define higher-order differentials, like D p,p,q F(p, q) or D p,p,q,q F(p, q).

The dimensions of these quantities must be kept in mind in the following equations: D p τ F(p, q) is a scalar, ∇ p F(p, q) is a 3D vector, D p,q F(p, q) a 3 by 3 matrix, D p,p,q F(p, q) is a 3x3x3 third-order tensor, and ∇ p,p,q,q F(p, q) a 3x3x3x3 fourth-order tensor.

Successive derivations of the covariance model K, which is a scalar function, are classically denoted by K ′ for the first derivative, and by K (n) for the n-th derivative.

B.2.2 Antisymmetry property

In the particular case where F can be written as a function of pq, i.e. there exists a function G : R 3 → R such that F(p, q) = G(pq), differentiating with respect to p is the same as differentiating with respect to q and changing the sign of the result.

In the present case, under our second-order stationary assumption, the covariance function of Z and the anisotropic distance satisfy this antisymmetry property, so we can substitute p for q in their differentials simply by changing the sign. The final sign of the substituted covariance function will then depend on the parity of the number of substitutions in the following way: (-1) n D p,p,... C Z (h pq ) = D q,q,... C Z (h pq ) where n is the number of substitutions.

B.3 First derivative of covariance function

Directional derivative at point p according to e

x Cov(D ex Z(p), Z(q)) = D p ex Cov(Z(p), Z(q)) = Cov( lim α→0 Z(p + α e x ) -Z(p) α , Z(q)) = lim α→0 1 α [Cov(Z(p + α e x ), Z(q)) -Cov(Z(p), Z(q))] = lim α→0 1 α [C Z (h pq + αe x ) -C Z (h pq )] = D p ex C Z (h pq ) = lim α→0 1 α [K(r(h pq + αe x )) -K(r(h pq ))] = D p ex K(r(h pq )).
(B.27) Thus, we obtain

D p ex Cov(Z(p), Z(q)) = D p ex C Z (h pq ) = D p ex K(r(h pq )). (B.28)
The gradient of C Z with respect to p is -K should be differentiable, -For ∇ p C Z (h pq ) to be defined where p = q, we need K ′ (0) = 0. Second-order derivative specification of C Z with respect to p then q

∇ p C Z (h pq ) =   D p ex C Z (h pq ) D p ey C Z (h pq ) D p ez C Z (h pq )   . (B.

B.3.1 Generalization: first-order derivative of the covariance function in a specific direction

Cov(D τ Z(p), Z(q)) = Cov(τ p .∇Z(p), Z(q)) = τ p .∇ p Cov(Z(p), Z(q)) = τ p . ∇ p C Z (h pq ). (
D p,q C Z (h pq ) = K (2) (r(h pq )) ∇ q r(h pq ). t ∇ p r(h pq ) + K ′ (r(h pq )) D p,q r(h pq ), (B.35) with D p,q r(h pq ) = - A r(h pq ) + A h pq t (A h pq ) r 3 (h pq ) ∈ M 3 (R). (B.36)
Let us notice again that

D p,q C Z (h pq ) = -D q,q C Z (h pq ) = -D p,p C Z (h pq ). (B.37)
Equation B.35 leads to 2 conditions -K should be twice differentiable, -For D p,q C Z (h pq ) to be defined at h pq = 0, we need K ′ (0) = 0. B.40) where D p,p,p C Z (h pq ) is a third-order tensor (3x3x3) given by

B.4.1 Second-order derivative of the covariance function in a specific direction

Cov(D 2 τ Z(p), Z(q)) = Cov( t τ .D 2 Z(p).τ, Z(q)) R 3 × R 3 → R = t τ.D p,p C Z (h pq ).τ = ∑ i ∑ j τ[i] τ[ j] D p,p C Z (h pq )[i, j]. (B.38) B.5 Third-order derivative of the covariance function Cov(D 2 τ Z(p), D µ Z(q)) = Cov( t τ .D 2 Z(p).τ , t µ .∇Z(q)) = ∑ i ∑ j ∑ k τ[i] τ[ j] µ[k] D p,p,q C Z (h pq )[i, j, k]. (B.39) Reminding that D p,p,q C Z (h pq ) = -D p,p,p C Z (h pq ), (
D p,p,p C Z (h pq )[i, j, k] = K ′ (r) D p,p,p r[i, j, k] + K (2) (r) (∇ p r[i]D p,p [ j, k] + ∇ p r[ j]D p,p [i, k] + ∇ p r[k]D p,p [i, j]) + K (3) (r) (∇ p r[i]∇ p r[ j]∇ p r[k]). (B.41)
The third-order derivative of r is also a third-order tensor (3x3x3) given by B.42) where t = Ah p,q B.6 Fourth-order derivative of the covariance function

D p,p,p r(h pq )[i, j, k] = - 1 r 3 [A i, j t k + A i,k t j + A j,k t i ] + 3 r 5 [t i t j t k ], (
Cov(D 2 τ Z(p), D 2 µ Z(q)) = Cov( t τ .D 2 Z(p).τ , t µ .D 2 Z(q).µ) = ∑ i ∑ j ∑ k ∑ l τ[i] τ[ j] µ[k] µ[l] D p,p,q,q C Z (h pq )[i, j, k, l]. (B.43)
Reminding that D p,p,q,q C Z (h pq ) = D p,p,p,p C Z (h pq ), (B.44) D p,p,p,p C Z (h pq ) is a 3x3x3x3 fourth-order tensor B.47) where t = Ah p,q .

D p,p,p,p C Z (h pq ) = K ′ (r) M 1 + K (2) (r) M 2 + K (3) (r) M 3 + K (4) (r) M 4 , (B.45) with M 1 = D p,p,p,p r M 2 [i, j, k, l] = D p,p r[i, j] D p,p r[k, l] + D p,p r[i, k] D p,p r[ j, l] + D p,p r[i, l] D p,p r[ j, k] + ∇ p r[i]D p,p,p r[ j, k, l] + ∇ p r[ j]D p,p,p r[i, k, l] + ∇ p r[k]D p,p,p r[i, j, l] + ∇ p r[l]D p,p,p r[i, j, k] M 3 [i, j, k, l] = D p,p r[i, j] ∇ p r[k] ∇ p r[l] + D p,p r[i, k] ∇ p r[ j] ∇ p r[l] + D p,p r[i, l] ∇ p r[ j] ∇ p r[k] + D p,p r[ j, k] ∇ p r[i] ∇ p r[l] + D p,p r[ j, l] ∇ p r[i] ∇ p r[k] + D p,p r[k, l] ∇ p r[i] ∇ p r[ j] M 4 [i, j, k, l] = ∇ p r[i] ∇ p r[ j] ∇ p r[k] ∇ p r[l] (B.46) and D p,p,p,p r(h pq )[i, j, k, l] = - 1 r 3 (A i j A kl + A ik A jl + A il A jk ) + 3 r 5 (A i j A kl t k t l + A ik A jl t j t l + ... + A kl A i j t i t j ) - 15 r 7 (t i t j t k t l ), (
C Demonstration of null second-order derivative along the hinge line

If we consider a twice-differentiable parametric curve γ : t

→   x(t) y(t) z(t)   ∈ R 3 and define g : t → Z(γ(t)),
where Z is our potential field (Z : R 3 → R), we have in all generality g ′ (t) = ∇Z(γ(t)).

--→ γ ′ (t) (C.48) and g ′′ (t) = ∇Z(γ(t)).

--→ γ ′′ (t) + t --→ γ ′ (t) . D 2 Z(γ(t)) .

--→ γ ′ (t).

(C.49)

Let us assume that γ is a parametrization of a planar curve included in an isopotential of Z, and that furthermore, there exists a plane that is tangent to this isopotential along the curve defined by γ. This is the case, for instance, for any parametrization of the H 3 curve of Fig. 1.

Under this assumption, g is constant, therefore in particular g ′′ = 0. Since --→ γ ′′ (t) is in the same plane as γ (plane P t in the case of Fig. 1), that is tangent to the isopotential, it is orthogonal to the gradient ∇Z(γ(t)).

Equation C.49 then reduces to t --→ γ ′ (t)∇ 2 Z(γ(t)) --→ γ ′ (t) = 0, i.e.

--→ γ ′ (t) ∈ Ker(D 2 Z(γ(t))). In other words, the directional second derivative of Z is null along --→ γ ′ (t), which is the direction of the crest/hinge line H 3 in the particular case of Fig. 1.

D Volumic error indicator

D.1 Volumic error (V e ) computation

The indicator V e provides a quantitative measure of the ability of the model to reconstruct the structure of interest of our synthetic example. It can be thought of as the volume between the two surfaces near the fold.

Let us denote by z(x, y) and ẑ(x, y) the respective elevations of the theoretical surface and the modeled surface at location (x, y). The isopotential considered in this work, described in Sec. 3 has a constant elevation of z f loor = 0.6 on a large part of the domain, but no information is provided to any model in this zone. In order not to penalize errors made far from the region of interest, instead of considering ẑ, we consider min(ẑ, z f loor ) when calculating the volume V e = zmin(ẑ, z f loor ) 1 = |z(x, y)min(ẑ(x, y), z f loor )|

D.2 Uncertainty (∆V ) on this computation

In practice, the interpolated value of the potential Z * (x, y, z) is known only on the grid points (x i , y j , z k ); at each location (x i , y j ) the elevation ẑ(x i , y j ) is estimated by taking z k such that Z * (x i , y j , z k ) is closest to 2.5, which is the value of the iso-surface considered here. This leads to an estimation uncertainty ∆V that depends on the vertical resolution (the finer the grid, the lower this uncertainty). For each model, this uncertainty is computed along with the volumic error indicator (it is almost constant but may vary a little, because this uncertainty is zero for locations where ẑ ≤ z f loor , since we use the value z f loor anyway). Several models are close in performance, in terms of this indicator; this uncertainty needs to be taken into account when comparing them. D.3 Relative volumic improvement (RVI, ηV e ) S1 : Scenario 1 ; S2 : Scenario 2 V e : Volumic error (u 3 ) V e 1 : Volumic error of initial result of Scenario 1 (u 3 ) V e 2 : Volumic error of initial result of Scenario 2 (u 3 ) where u is the arbitrary length unit used here. First experiment (Fig. 5) Second experiment (Fig. 6) Third experiment (Fig. 7) Fourth experiment (Fig. 8)

V
V
V

Figures

Figure 1 Fold Figure 2 Please see the Manuscript PDF le for the complete gure caption.

Figure 3 Please see the Manuscript PDF le for the complete gure caption.

Figure 4 Initial results of the potential eld method on two different initial samplings. Sampling 1 consists in 4 contact points (blue spheres) and 4 gradient data (black arrows). Sampling 2 consists in 8 contact points and 8 gradient data. Comparison of the impact of ( rst row) 1, (second row) 2, (third row) 3 and (last row) 4 constraints per point at the same 6 locations on the fold axis of the theoretical surface

Fig. 1 :Fig. 2 :

 12 Fig. 1: Fold geometry and planes of interest. Sa: Axial surface. n: Pole vectors belonging to the axial plane; n ′ : Other pole vectors of the same profile. Pt : Tangent plane. S 0 : Bedding. S 1 : Cleavage. H 1 , H 2 , H 3 : Hinge lines of beddings L 1 , L 2 , L 3 . The hinge line H 3 and the crest line coincide in this case. D A and D A ′ : Fold axis measurements on hinge lines H 2 and H 3 . D B : Intersection lineation between S 0 and S 1 measurement. D C : Fold axis inference from stereographic projections using n ′ pole measurements.

  Fig.3: Top (a) and side (b) view of the iso-surface Z = 2.5 of the potential of Eq. (7). The elliptic hinge line of the fold, with semi-major and semi-minor axes, A and B, is presented on (a) (dashed line). Random data for different samplings are generated with the angle θ around this structure. Data sampling zones, A1 and A2, with their respective cut off angles, c 1 = 2π 24 , c 2 = 4π 24 , c ′ 1 = 8π 24 , c ′ 2 = 10π 24 , are also shown, as well as the vertex of the ellipse.

Fig. 4 :

 4 Fig.4: Initial results of the potential field method on two different initial samplings. Sampling 1 consists in 4 contact points (blue spheres) and 4 gradient data (black arrows). Sampling 2 consists in 8 contact points and 8 gradient data.

  Fig. 5: Impact of (c)(d) 15 tangent data (green arrows), (e)(f) 15 null second-order derivative data (D2, pink arrows), (g)(h) both data on the resulting surface compared to (a)(b) the initial result with Sampling 1 or Sampling 2

Fig. 6 :Fig. 7 :

 67 Fig. 6: Impact of the number of tangent data along the hinge line in the two different scenarios

  Fig. 8: Comparison of the impact of (first row) 1, (second row) 2, (third row) 3 and (last row) 4 constraints per point at the same 6 locations on the fold axis of the theoretical surface

  21)Stationarity case of order two assumptionCov(Z(p), Z(q)) = C Z (pq) = C Z (h pq ), (B.22) with h pq → C Z (h pq ) R 3 → R. (B.23)Case of constant anisotropic distance (global anisotropy)C Z (h pq ) = K(r(h pq )),

  ηV e (RVI) : Signed relative difference of Ve compared to initial result in the same scenario (%)

  geometry and planes of interest. Sa: Axial surface. n: Pole vectors belonging to the axial plane; n0: Other pole vectors of the same pro le. Pt : Tangent plane. S0: Bedding. S1: Cleavage. H1, H2, H3: Hinge lines of beddings L1, L2, L3. The hinge line H3 and the crest line coincide in this case. DA and DA0 : Fold axis measurements on hinge lines H2 and H3. DB: Intersection lineation between S0 and S1 measurement. DC: Fold axis inference from stereographic projections using n0 pole measurements.

Figure 5 Figure 6

 56 Figure 5

Figure 7

 7 Figure 7

  Second-order derivative of the covariance functionCovariance of partial derivatives w.r.t. p according to e x and w.r.t. q according to e y Cov(D ex Z(p), D ex Z(q)) = t e x .D p,q C Z (h pq ).e y .(B.34)

	Title Suppressed Due to Excessive Length	29
	B.4
		B.33)

  This part summarizes all the intermediate computations for the results concerning the volumic error of the different investigations of Sec. 4.

		Initial V e in S1 (u 3 )	Initial V e in S2 (u 3 )
		104.68 ± 4.12	48.36 ± 3.40
			Initial values
		V e in S1 (u 3 )	V e in S2 (u 3 )	ηV e to S1 (%)	ηV e to S2 (%)
	Tangent	85.19 ±4.27	50.72 ±3.40	18.62	-4.89
	D2	91.52 ±4.17	49.76 ±3.40	12.58	-2.91
	Tangent + D2	91.39 ±4.35	49.01 ±3.40	12.70	-1.35

Pizzella et al.


Wafa Ben Hassen

Fabrice Auzanneau
email: fabrice.auzanneau@cea.fr

Marc Olivas Carrion

François Pérès

Wafa Ben Hassen
email: wafa.benhassen@cea.fr

François Peres
email: francois.peres@enit.fr

Ayeley Tchang Ani
email: ayeley.tchangani@enit.fr

A distributed diagnosis strategy using bayesian network for complex wiring networks

Keywords: complex wiring network, distributed strategy, reflectometry, Bayesian Networks, uncertainty estimation

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In recent years, we have witnessed a paradigm shift in the automotive and aeronautic industries. The appear ance of "X-by-wire" technology, replacing the traditional mechanical and hydraulic control systems with electronic ones, increases embedded systems complexity. Moreover, cumulated wire length rise up to 4 km in a modern car and 400 km in an aircraft AUZANNE AU [2012]. A part of these wiring networks have been identified as likely responsible of several tragic mishaps such as the crash of TWA Flight 800 in July 1996 and Swissair Flight 111 in September 1998. Monitoring, diagnosis and maintenance are consid ered as a nightmare when wire harness is out-of-sight or unreachable without dismantling the external structure. It is also highly expensive in terms of money, time, resources, risk and persons. Each year, US Navy spends 1.8 million person-hours for its aircraft wiring systems maintenance FURSE and HAUPT [2001]. Reflectometry is a powerful technique for detection and localisation of electric faults in wiring networks. The idea is to inject a wideband test signal down the wire. During its propagation, a part of its energy reflects back to the injection port when it crosses impedance discontinuities. Then, a received signal analysis, called "reflectogram", gives location information about the detected wire fault. In the literature, many types of reflectometry techniques are proposed. The main difference between these varieties lies on the injected test signal type and reflected signal pro-cessing techniques. Time Domain Reflectometry (TDR) uses a pulse (Gaussian, rectangular,etc.) as a test signal and measures the time delay between the injected and received signal to determine the fault location. Frequency Domain Reflectometry (FDR) uses a set of sinusoïdal waves as a test signal and analyses the phase and/or magnitude of the reflected wave FURSE et al. [2003]. Although FDR and TDR have proved their efficiency in fault detection and localisation for single wires, they are not so efficient in complex wiring networks. Indeed, both need to use a high voltage signal to detect intermittent faults. However, the test signal levels may interfere with the native ones, if applied while the target system ( aircraft, automobile, etc) is running. For those reasons, researchers have focused on innova tive methods for "on-line" diagnosis in wiring networks aiming at reducing the consequences of intermittent and arc faults. To do so, other types of reflectometry have been proposed such as Sequence Time Domain Reflectom etry (STDR) SHARMA et al. [2007] and Spread Spec trum Time Domain Reflectometry (SSTDR) that inject a pseudo-noise code TAYLOR and FAULKNER [1996], WEI and LI [2011] to diagnose wires in real time. These methods are able to detect and locate faults on live wires even when the test signal level is well below the noise margin of the signal already on the wire. In order to limit the bandwidth allocated to system signal and reduce inter-signal and Electromagnetic Compatibil ity (EMC) interference, Multiple Carrier Reflectometry


Eric Villeneuve

François Pérès

Laurent Geneste

An approach to improve risk assessment using experience feedback

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible Any correspondence concerning this service should be sent to the repository administrator: tech-oatao@listes-diff.inp-toulouse.fr This is an author's version published in: http://oatao.univ-toulouse.fr/ 28112

To cite this version:

Villeneuve, Eric and Pérès, François and Geneste, Laurent An approach to improve risk assessment using experience feedback. (2010) In: IFIP Doctoral Spring Workshop "Product and Asset Lifecycle Management, January 2010 (Huelva, Spain).

Open Archive Toulouse Archive Ouverte

!" #$ % & ' ()*+ ),-() .' / 0 1 2 3 4 ! " !" # ! $ $ % $ & & & ' $ ' $ $ $ 5 3


Khouloud Jlaiel
email: kjlaiel@enit.fr

Malik Yahiaoui

Jean-Yves Paris

Jean Denape

Acoustic signature identification of damage and wear mechanisms in a steel/glass sliding contact

Keywords: Reciprocating sliding, acoustic emission, real-time observation, Hertzian cracks, third body, accommodation mechanisms

The tribological behavior of a steel/glass ball-on-flat contact was studied by synchronizing the friction measurements with an acoustic emission device and a vision system. The results highlight two distinct friction regimes identified with low and high friction values. Their transition is characterized by a modification of acoustic emission signals. In addition, two main damage and wear mechanisms are identified: the creation and propagation of Hertzian cracks visible on the glass surface and the constitution of an interfacial layer of debris. The different accommodation mechanisms, activated successively or simultaneously, are identified for acoustic emission frequencies between 300 and 700 kHz. Eventually, this approach allows a real-time wear mechanisms identification and gives better insights about acoustic emission signals in relation to tribological systems.

Introduction

During friction, various tribological mechanisms can occur simultaneously and it is essential to better understand the evolution of these mechanisms from the contact setting to the establishment of a tribofilm, transfer film, and/or debris particles. In the common case of nontransparent frictional surfaces, this is a challenge due to the inaccessibility of the contact area during testing. The first bodies post mortem analysis gives information only on the final state of the contact area after multiple wear mechanism changes.

Over the last few decades, research has adopted several approaches to observe tribological processes. In situ techniques have been used to observe sliding dynamics, identify mechanical and/or chemical changes near the surface, and better understand the behavior of the third body for dry as well as lubricated contacts using optical microscopy, Raman spectroscopy, and camera observation. [START_REF] Yahiaoui | Forced oscillations dynamic tribometer with real-time insights of lubricated interfaces[END_REF][START_REF] Sriraman | Tribological behavior of electrodeposited Zn, Zn-Ni, Cd and Cd-Ti coatings on low carbon steel substrates[END_REF][START_REF] Eguchi | Measurement of real contact area and analysis of stick/slip region[END_REF][START_REF] Dvorak | In situ analysis of third body contributions to sliding friction of a Pb-Mo-S coating in dry and humid air[END_REF][START_REF] Teng | In situ observations of fretting wear behavior in PMMA/steel model[END_REF] The acoustic emission (AE) technique, defined as the creation of transient elastic waves generated by local micro-displacements within a material, [START_REF] Standard | Non-destructive testing--Terminology Part 9: Terms used in acoustic emission testing[END_REF] can be utilized to reveal and monitor the tribological processes involved. In other words, it consists of measuring the elastic wave signal that is released by various sources in materials under stress (e.g. fracture, plastic deformation, and plowing). The underlying premise is that each elementary contact mechanism has a specific AE signature, which could be characterized in terms of released energy and frequency.

In the literature, the AE technique has shown its sensitivity to detect unusual changes in friction systems. For example, Hase et al. [START_REF] Hase | Acoustic emission signals and wear phenomena on severe-mild wear transition[END_REF] showed that AE signals are better suited for the identification of the different wear phenomena observed on worn steel surfaces. The AE signals are affected when there is a change in the wear mode and type of wear particles caused by variations in the sliding velocity. Moreover, the same authors found that a continuous AE signal of low amplitude is detected by the generation of slip lines in the frictional interface and by the generation and transfer of wear particles. [START_REF] Hase | Acoustic emission in elementary processes of friction and wear: In-situ observation of friction surface and AE signals[END_REF] Meriaux et al. [START_REF] Meriaux | Identification of fretting fatigue crack propagation mechanisms using acoustic emission[END_REF] identified three successive stages in the global crack propagation process in a fretting fatigue test. They demonstrate that short and inclined cracks initiate and propagate first due to the shearing process (mode II). Thus, AE cumulative activity increases gradually. Second, cracks propagate in a mixed mode of shearing and tension where the AE cumulative activity stabilizes. Finally, fatigue cracks propagate by tension in a pure mode I process where the AE cumulative activity reaches its maximum value. Similarly, Rastegaev et al. [START_REF] Rastegaev | Using acoustic emission signal categorization for reconstruction of wear development timeline in tribosystems: case studies and application examples[END_REF] explained that an increasing severity of wear is accompanied by a concomitant increase in the AE energy, which occurs in parallel with the decrease in the median frequency of the AE power-spectral density. Recently, Haas et al. [START_REF] Haas | Examination of undisturbed acoustic emission generated by experimentally modelled two-body abrasive wear events[END_REF] identified distinctive AE frequency ranges in plowing and cutting events. Also, Jlaiel et al. [START_REF] Jlaiel | Tribolumen: A tribometer for a correlation between AE signals and observation of tribological process in real-time-application to a dry steel/ glass reciprocating sliding contact[END_REF] observed a change in the AE signals associated with low friction, Hertzian cracks formation, and third body formation. Regarding the frequency band, a correlation map for the AE frequency band of several tribological phenomena has been established by Hase et al. [START_REF] Hase | Correlation between features of acoustic emission signals and mechanical wear mechanisms[END_REF] based on a literature review. This map allocates frequency bands and amplitudes to various interfacial accommodation modes and sources and even operating conditions. For instance, abrasive wear is identified between 0.2 and 1 MHz. On the same map, the authors represent a sliding friction frequency band under 200 kHz and a crack propagation one between 100 and 700 kHz. More surprisingly, large frequency bands related to fatigue and tensile testing are observed on the map at 600 kHz. These large frequency bands obviously enclose several elementary accommodation mechanisms of the tribosystem which are difficult to distinguish during the test and in post mortem analysis of the worn surfaces.

Eventually, most of the past studies propose correlations between AE signals and interfacial mechanisms without a direct and synchronous observation of these mechanisms. The main goal of this study is then to clearly correlate the AE signatures of the elementary wear damage and the accommodation mechanisms occurring within the contact using a real-time observation of the contact area between a steel ball and glass flat.

Materials and methods

The Tribolumen tribometer used in this study allows simulating the fretting phenomena and reciprocating sliding close to fretting conditions on transparent materials. [START_REF] Jlaiel | Tribolumen: A tribometer for a correlation between AE signals and observation of tribological process in real-time-application to a dry steel/ glass reciprocating sliding contact[END_REF] This apparatus ensures the contact between a soda-lime glass sample and a 100Cr6 steel ball (AISI 52100, 6 mm in diameter, 0.015 ± 0.005 µm of average roughness) in a ball-on-flat configuration. The dry friction test was carried out under pure sliding conditions with an imposed stroke of 100 µm (i.e. an amplitude of ±50 µm), an excitation frequency of 10 Hz, and a normal load of 5 N, which correspond to a maximum contact pressure around 700 MPa (according to Hertz' theory). The maximum velocity relative to the glass sample displacement is around 9.36 mm/s. A high-speed camera is positioned on the soda-lime glass sample side. The camera is used at a frame rate of 500 fps to observe the contact area in real-time testing. The AE acquisitions were performed using two piezoelectric sensors fixed close to the contact. The sensors were set using a water-based adhesive containing a styrene acrylic copolymer. This adhesive allows a good transmission of acoustic signals and dries in a few seconds and the sensors are easily detached after the experiments. The sensors (Micro-80 sensor from Mistras) have a large frequency band with a maximum sensitivity between 100 and 1000 kHz. The coupling of the sensors was then controlled by the Hsu-Nielsen source method described by the standard NF EN 1330-9. A threshold value of 27 dB was chosen to get finer detection without being disturbed by environmental conditions. The device synchronizes AE signals with friction data acquired by the tribometer. Before testing, the samples were cleaned for 15 min using an industrial detergent, then rinsed with demineralized water, and finally ultrasonically cleaned with ethanol for 10 min before being dried in an oven at 60 °C for 15 min.

The measurements of displacement and tangential and actuator forces were acquired using an HBM Quantum X system. In parallel, the measured acoustic signals that include various characteristic parameters to be considered were analyzed. They were directly processed using Mistras NOESIS™ software. The contact area image sequences were obtained through a high-speed camera and were analyzed using the ImageJ software. All the experiments were repeated three times and similar results were obtained.

Results and discussions

Frictional results

The experiments display classical fretting logs (Figure 1(a)). As the displacement is imposed, the fretting logs display a quasi-square section shape. The friction is relatively low at the beginning of a test and increases after 200 cycles to reach a steady value until the end.

Figure 1(b) presents the coefficient of friction µ as a function of the cycle number. At the very beginning, the coefficient of friction is around 0.12. After about 200 cycles, this value increases abruptly to reach 0.8. Finally, after 420 cycles, the coefficient of friction stabilizes around 0.68.

Wear results

The worn surfaces of the soda-lime glass and the 100Cr6 steel ball surfaces were observed systematically by optical microscopy without cleaning after 1200 cycles (Figure 2). The wear scar on the flat glass sample shows an expected elliptical shape with the smallest axis corresponding to the diameter of the wear scar on the steel ball (around 178 ± 1 µm) and the largest axis oriented in the sliding direction. The central zones of the two worn surfaces have abrasion scratches. In the third body approach, this zone corresponds to the interface generating debris associated with the source flow and the internal flow. Debris is preferentially agglomerated on the soda-lime glass and ejected outside of the contact area, in the displacement direction, representing the wear flow.

Figure 3 shows a typical scanning electron microscopy observation of the third body after separating the ball and the flat. It is noted that debris is compacted together and adheres to the surface of the soda-lime glass sample to form an adherent interfacial layer. When the contact is opened and the ball is removed, the breakdown of this interfacial layer inevitably occurs as a result of the separation of tensile stress following the end of the tribological test (Figure 3(a)). Outside the wear track, ejected debris is found in fine particles with a diameter of 1 ± 0.5 µm and much larger particles with an average diameter of 8.5 ± 0.5 µm (Figure 3(b)).

The energy dispersive spectroscopy maps show that iron elements were transferred from the ball to the flat (Figure 4). These elements are more numerous in the central wear track than in the external wear track.

AE results

Several characteristic parameters of AE signals were calculated and extracted from the recorded waveforms. For this study, only AE amplitude (A AE ), absolute energy (E abs ), and centroid frequency ( f c ) will be discussed as the most representative parameters. The AE amplitude represents the maximum amplitude of the signal during the test. The absolute energy represents the integral of the squared signal over the burst duration. The centroid frequency represents the frequency center of mass of the AE signals, that is, it characterizes the overall frequency content of an AE signal.

Figure 5 presents the variations of the amplitude, the centroid frequency, and the cumulative absolute energy during the friction test. Three distinct stages are identified synchronously with the evolution of the coefficient of friction: • In the final stage (420 cycles to the end), the coefficient of friction is again relatively stable but high. During this third phase, the AE signals are low in energy (about 10% of the total AE energy emitted) and are characterized also by low amplitudes (29-35 dB), but high frequencies (500-700 kHz), compared with the second stage.

•

Correlation between tribological behavior and AE signals

In this section, to better understand the phenomena occurring during the test, a representation of the coefficient of friction (Figure 6) is accompanied by images extracted from the video showing the evolution of damage processes (Figure 7) at specific cycle numbers noted with letters from a to j. Figure 7(a) shows the initial contact between the ball and the flat under a normal load of 5 N. The central, circular dark area is the contact area. Following Hertz's theory, the radius of this contact area is determined and is about 58 µm, as can be observed. It is formed by elastic deformation. The white rings around this area also called Newton's rings, corresponding to the optical interference fringes. They are caused by the divergence of the ball away from the glass sample surface.

In the first stage (until 200 cycles), the coefficient of friction remains very low (about 0.12) and no degradation in steel or glass surfaces is observed (Figure 7 After 200 cycles (beginning of the intermediate stage), natural screens have been destroyed by friction. Direct interactions between the surfaces of steel and glass samples occur, resulting in the creation/breakdown of adhesive junctions and therefore in a rise of coefficient of friction. This mechanism generates AE signals characterized by an acoustic amplitude close to 49 dB and a centroid frequency of 366 kHz. These results are similar to those observed by Kolubaev et al. [START_REF] Kolubaev | Analysis of acoustic emission during sliding friction of manganese steel[END_REF] They show that when the real contact occurs between the conjugated surfaces and when the contact area increases, both the coefficient of friction and AE signal increase. The acoustic frequency reached in this case is around 340 kHz. Furthermore, as reported by Hase et al., [START_REF] Hase | Correlation between features of acoustic emission signals and mechanical wear mechanisms[END_REF] the acoustic Moreover, according to Sun et al., [START_REF] Sun | Wear monitoring of bearing steel using electrostatic and acoustic emission techniques[END_REF] as soon as the sliding contact is effective, a lot of debris can easily be generated which leads to an increase of the AE signals.

At 233 cycles, a first Hertzian crack is observed on the glass sample (Figure 7(e)). The crack has a semicircular shape in which the concave side indicates the sliding direction when it occurred. This type of crack is usually observed on brittle materials such as glasses and ceramics. It originates from the stress distribution in the flat sample, which is induced by the ball displacement. For a sliding ball-on-flat contact, tensile stress appears behind the contact and leads to a fracture. In this stage, the centroid frequency is around 501 kHz. However, the AE amplitude increases to 37 dB, because the strain energy released by cracking is higher than that released by the formation of wear debris.

When the ball moves, the number of cracks increases. Figure 7(f) shows a series of three consecutive cracks at the back of the contact. In addition, with the reciprocal sliding of the ball, cracks can also appear completely or partially in the opposite sliding direction. The cracks are equally spaced (around 36 ± 1 µm). This inter-crack spacing is generally related to the coefficient of friction and the fracture toughness of the glass. [START_REF] Hills | Brittle fracture from a sliding Hertzian contact[END_REF][START_REF] Geng | Mechanical behavior of stiff coating on glass under sliding contact[END_REF] Friction also produces compressive stress in the front of the sliding surface that closes the cracks generated during the previous part of the cycle when the displacement followed the opposite direction, so they are not more visible on the video. Such cracks are visible again when they are re-opened by tension behind the contact when the sliding is once more in the appropriate direction.

Figure 7(g) shows that detached particles continue to progress in terms of width. Their accumulation produces an adherent interfacial bed of powder that reduces the first bodies' surface interactions. A part of the detached particles is also ejected outside the contact zone and creates the wear flow Q w .

At 345 cycles (Figure 7(h)), the abrupt increase of the coefficient of friction up to 0.8 is accompanied by the increase of the AE amplitude to 55 dB. This is due to the appearance of the largest Hertzian crack on the glass surface. In this intermediate stage of testing, cracks length and opening are measured from the in situ observation. The maximum value of the crack opening reaches 10.8 ± 0.2 µm and follows the coefficient of friction evolution (Figure 8(a)). In contrast, the crack length increases monotonously from 44 to 155 ± 5 µm and follows the evolution of the acoustic cumulative E abs (Figure 8(b)).

The detection of crack initiation and propagation during sliding contact has been investigated by Cadario and Alfredsson [START_REF] Cadario | Fretting fatigue crack growth for a spherical indenter with constant and cyclic bulk load[END_REF] They showed that an increase in acoustic events and a variation of the AE signal slope are related to crack growth. Moreover, Ramadan et al. [START_REF] Ramadan | Detection of stress corrosion cracking of high-strength steel used in prestressed concrete structures by acoustic emission technique[END_REF] studied the AE originating from crack initiation and propagation. The frequency band associated with these phenomena is from 0.05 to 0.7 MHz. The AE amplitude attributed to the crack initiation was in the range between 23 and 30 dB. After that, the AE amplitude that is between 30 and 40 dB is associated with crack growth. The AE amplitude that is >55 dB is associated with material failure. All these results are close to the results presented above.

In the final stage, starting at 420 cycles (Figure 7(i) and (j)), the coefficient of friction starts to stabilize at 0.68. The adherent interfacial layer continues to grow until it becomes extended to the entire contact area. Newton's rings and also cracks are no more visible. The AE signals are low compared to the second stage. This is similar to Ito et al. [START_REF] Ito | The relationship between AE and dissipation energy for fretting wear[END_REF] results. They showed that in a steady-state stage the accumulation of wear debris between the contacting surfaces leads to relatively low values of AE signals despite a high coefficient of friction.

The shearing accommodation in this layer corresponds to a centroid frequency around 650 ± 1 kHz. In the same way, the crack sizes (length and opening) start to stabilize, respectively, around 160 ± 5 and 6 ± 0.2 µm. The accommodation mechanism is then made by crack opening/closing characterized by a centroid frequency at around 550 ± 1 kHz.

Conclusions

Ball-on-flat dry friction behavior of soda-lime glass against 100Cr6 steel was studied using an AE technique associated with direct observation by video through the contact. The study highlights the importance of synchronizing friction measurements with real-time observation and AE signals continuously recorded during friction testing. The AE technique shows high sensitivity to tribological phenomena and has made it possible to detect several sources of AE generated during testing. These sources can be associated with damage or wear mechanisms: crack initiation and propagation, creation/breakdown of adhesive junctions (particle detachment), and shear accommodation mechanism in the interfacial layer. This approach allows a real-time wear mechanism identification and gives better insights into the AE signals in correlation with the tribological system.

The results show that:

  In the initial stage (until 200 cycles), the coefficient of friction was low and no AE signals were recorded over the experimental AE threshold detection level. • In the intermediate stage, from 200 cycles until 420 cycles, the coefficient of friction increases abruptly. During this second phase, AE signals (hit) are recorded. They represent about 90% of the total AE energy emitted and are characterized by a wide range of frequencies (360-700 kHz) and amplitude (29-55 dB).

Figure 1 .

 1 Figure 1. (a) Fretting log: tangential force versus number of cycles and (b) coefficient of friction sliding average for soda-lime glass against 100Cr6 ball steel (F N = 5 N, δ = ± 50 µm, f = 10 Hz).

Figure 2 .

 2 Figure 2. Optical microscopy of worn surfaces of (a) soda-lime glass and (b) 100Cr6 steel ball after 1200 cycles (F N = 5 N, δ = ± 50 µm, f = 10 Hz). ↔ Sliding direction.

Figure 3 .

 3 Figure 3. Scanning electron microscopy (SEM) observation of the third body (a) in the central wear track and (b) in the external wear track.

Figure 4 .

 4 Figure 4. SEM observation and EDS analysis (Si and Fe) of the third body. SEM, scanning electron microscopy; EDS, energy dispersive spectroscopy.

  (b)). So, the accommodation has taken place by shearing in residual natural screens of the first bodies' surfaces. No AE signals were recorded over the experimental AE threshold.

Figure 5 .

 5 Figure 5. Evolution of (a) acoustic amplitude and (b) cumulative absolute energy during testing.

Figure 6 .

 6 Figure 6. The evolution of the coefficient of friction.

Figure 7 .

 7 Figure 7. Images of the contact area extracted from the video showing the evolution of the damaging process during sliding friction. (a) and (b) in the first stage, from (c) to (h) in the intermediate stage, (i) and (j) in the final stage.

Figure 8 .

 8 Figure 8. (a) Synchronized evolution of the coefficient of friction and crack opening during the test and (b) synchronized evolution of cumulative energy and crack length during testing.

•

  The accommodation mechanisms taking place by shearing in natural screens of the first bodies do not create AE signals over the experimental threshold. • The direct interactions between the first bodies' surfaces (creation/breakdown of adhesive junctions) produce AE signals with a centroid frequency around 366 ± 1 kHz. • The initiation of the first Hertzian cracks is characterized by a centroid frequency around 501 ± 1 kHz. • In the stabilization period, the accommodation mechanism takes place by shearing in the interfacial layer at a centroid frequency around 650 ± 1 kHz. • In the stabilization period, the accommodation mechanism also takes place by cracks opening and closing at a centroid frequency around 550 ± 1 kHz. • Regarding crack sizes, the largest crack is characterized by the highest AE features. • The crack opening follows the coefficient of friction evolution. • The crack length follows the acoustic cumulative E abs .

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.