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

When too few field measurements are available for the geological modeling of complex folded structures, the results of implicit methods typically exhibit an unsatisfactory bubbly aspect. However, in such cases, anisotropy data are often readily available but not fully exploited. Among them, fold axis data are a straightforward indicator of this local anisotropy direction. Focusing on the so-called potential field method, this work aims to evaluate the effect of the incorporation of such data into the modeling process. Given locally sampled fold axis data, this paper proposes to use the second-order derivatives of the scalar field in addition to the existing first-order ones. The mathematical foundation of the approach is developed, and the respective efficiencies of both kinds of constraints are tested. Their integration and impact are discussed based on a synthetic case study, thereby providing practical guidelines to geomodeling tool users on the parsimonious use of data for the geological modeling of complex folded structures.

geomodeling tool users on the parsimonious use of data for the geological modeling

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 (Wellmann and Caumon 2018).
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 (Mallet 1992). 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 past 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 (Frank et al. 2007;Caumon et al. 2012) or not (Houlding 1994). 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 (Cowan et al. 2003;Turk and O'brien 2005;Hillier et al. 2014), the moving least squares functions Manchuk and Deutsch 2019;Renaudeau et al. 2019 or the dual form of kriging with the use of covariance functions (Lajaunie et al. 1997;Calcagno et al. 2008). Following in the footsteps of Lajaunie et al. (1997), the approach developed in this paper is of a 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 (Hillier et al. 2014). To deal with this issue, other methods typically introduce inferred or interpreted locally varying anisotropies (Boisvert et al. 2009). 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 (Ramsay and Huber 1987). 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, Maxelon et al. (2009) interpret variations in the orientation of the dominant planar fabrics. Laurent et al. (2016) suggest integrating additional field data describing folds, lineations and vergence. This approach and the inference of its parameters are discussed in Grose et al. (2017). In another vein, de la Varga and Wellmann (2016) 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 (Hillier et al. 2014). It turns out that working with the direction of this line has a practical asset: it allows us 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 axis 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 Sect. 2 by summarizing the theoretical background of the potential field 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 Sect. 3. It is used to illustrate the assets and limitations of using these mathematical objects in Sect. 4, thereby leading to the consideration of the maximum curvature criterion. These results, finally discussed in Sect. 5, should allow geomodeling software developers or users to get more insight on fold axis incorporation for the anisotropy modeling process.

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 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 isosurfaces S α They are the surfaces of interest, corresponding to the geological interfaces of the three-dimensional structural model. The so-called potential field method (Lajaunie et al. 1997) directly considers increments of the scalar field sampled on isosurfaces. The implicit function is then defined up to an additive constant. It can be set by choosing a reference isosurface, the potential value of which is arbitrarily set to zero (Calcagno et al. 2008).

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 can be grouped into two main categories: position point information (p ∈ R 3 ) and orientation vector data (d ∈ R 3 ). With regard 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 the same potential value, and therefore the difference of their potential value (called increment) is null 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 ii. By definition, the first-order derivative of Z at point p in a direction d is 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 isopotential 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 where ∇ Z (p) denotes the gradient of the scalar field, that writes 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 In practical applications and software packages, this inequality constraint is usually replaced by the (biased) assumption that the the gradient of the scalar field is of unit norm, with the consequence that one can write 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 Motivations for the addition of such data and its impact on the method are discussed in Sect. 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 where the left-hand side superscript t indicates the transpose of a vector or matrix.

The Interpolation Method: The Potential Field Method
Using a co-kriging estimator (Wackernagel 2013), the interpolated value of Z at point p, Z * (p), is computed by the weighted sum of different input data 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 (Chiles and Delfiner 2009), 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 The drift function m is usually expressed as a linear combination of n basis functions 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 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 Lajaunie et al. (1997), 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): 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 the second-order derivative of an arbitrary point of a hinge line (crest here) of the synthetic case (Sect. 3) are shown in 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 collinear, 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 isopotential surface coincides with the crest of the fold. However, the demonstration that the second-order directional derivative is null along a particular curve relies only 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 collinear. 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, 3) in all directions on the tangent plane. φ and φ + π are the directions along the hinge line for this given point. φ + π 2 and φ + 3π 2 are the orthogonal directions of the hinge line on the tangent plane on the resulting isosurface curvature. This could be significant information, especially in the case of very curvy hinge lines. Indeed, the geometric meaning of these two previous observations is then: 1. a null value of the second-order derivative of potential in the hinge line direction implies local flatness in the hinge direction; 2. the maximal value of the 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 the potential in all directions on the tangent plane already express local flatness in the hinge direction. However, the addition of a null second-order derivative value provides additional information on smoothness. Indeed, assuming that the potential Z is sufficiently smooth, its Taylor expansion at a point p on the hinge can be written as For h in the tangent plane, the first-order term ∇ Z (p).h is zero. If the hypotheses of Sect. 2.5.1 are satisfied and h is along the hinge, the second-order term 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. These data consist of 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 subparallel 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. The 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, and these observations are then quite rare, and thus valuable. 2. In fact, fold axis data come abundantly from intersection lineation 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. Stereographic projection tools can then 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 contribute significantly 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 microscale structures, parasitic folds, global regional context and expert knowledge. For a complete review of fold data see Hudleston and Treagus (2010) and McClay (2013).

Covariance Function Requirements
The implementation of the second-order derivative value of the potential gives rise to the need for 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 Wendland (1995), denoted by Wendland C4 in the following and which may be expressed as 3 Case Study

Presentation of the Considered Potential Field
In this paper, an explicit theoretical potential function has been chosen as a 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 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 isopotential surface of this field at level Z (x, y, z) = 2.5 is shown in Fig. 3: this surface is hereafter defined as our reference of study.

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 acceptable model. The following section will highlight interpolation issues encountered in such a sparse data environment.  7). The elliptic hinge line of the fold, with semi-major and semi-minor axes A and B, is presented in 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 cutoff 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

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). The resulting models of the potential field method on each sampling are presented in Fig. 4.
-Scenario 1 (Fig. 4a): The first sampling consists of four random contact points (with four associated gradient data points) generated in the A1 zone. -Scenario 2 (Fig. 4b): The second sampling is defined by four random contact points (with four associated gradient data points) generated in the A1 zone, and four random contact points/gradient data points generated in the A2 zone.
As a consequence, both sampling density and area coverage of the 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 The interpolation is performed on a 20 × 20 × 20 grid (although for numerical reasons the error score V e (Sect. 3.2.3) is computed on a 20 × 20 × 80 grid).

Fig. 4
Initial results of the potential field method on two different initial samplings. Sampling 1 consists of 4 contact points (blue spheres) and 4 gradient data points (black arrows). Sampling 2 consists of 8 contact points and 8 gradient data points 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 (Sect. 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 three-dimensional isosurface of interest is presented for each experiment hereafter. In addition, in order to aid 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 the improvements in the resulting model compared to initial results in the same scenario i (with its own volume error V e i ), the relative volume improvement (RVI) ηV e is introduced where V e i is the volume error of scenario i; i ∈ (1, 2).
This signed indicator ηV e quantifies the reduction in error compared to the initial model. It is expressed as a percentage. Negative values indicate deterioration compared to initial results, and 100% RVI means that the error is now null, whereas 0% RVI means no change in volume error value. Additional 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 into the scalar field variations. All resulting values for the different experiments are summarized in Appendix D.4.

Critical Case Versus Well-Constrained Scenario
The resulting surfaces of the potential field method with the two different samplings of Sect. 3.2.1 are shown in Fig. 4. This interpolation is performed with the choice of parametrization detailed in Sect. 3.2.2. These two starting assumptions lead to two different results, as follows.
The resulting surface with Sampling 1 (Fig. 4a) presents a "bubbly" aspect far from the data, as 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 volume error is V e 1 = 104.68. Moreover, this phenomena takes place here with the use of the Wendland C4 covariance model, but also occurs with the more common 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 volume 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 a better constrained model. Moreover, despite this prevailing error, the result may be qualitatively considered satisfying, and thus this scenario considered 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 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 (Sect. 2.5.4). The parsimonious deployment of both kinds of inputs is discussed in Sect. 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 using 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 (Lajaunie et al. 1997), 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 Sect. 2.5.4 are available. When tangent data are taken into account, there is still uncertainty about the necessary amount of data to constrain well poly-deformed folds in a sparse data environment, as well as a geographical distribution preference. These understandings would allow us 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 points, a null second-order derivative value and a maximal second-order 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 Sect. 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 Sect. 3.2, even though its orientation is no longer indicated. The different volume error indicators are listed in Appendix D.4.

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 Sect. 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. Fifteen homogeneously distributed vectors are sampled along the hinge line and expressed in terms of tangent (Fig. 5c, d), null second-order derivative (D2) (Fig. 5e, f) or both (Fig. 5g, h). For comparison, the initial results from Sect. 3.2.4 are also recalled (Fig. 5a, b). 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, and 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, and RVE is around 19%. Finally, the combination of both constraints (Fig. 5g) leads to a qualitatively similar result as the one with tangent data constraints only. However, V e is 91 ± 4, and 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 a second-order derivative value made in Sect. 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 on ambiguous surfaces in other surface modeling applications.

Quantity or Spatial Distribution of Tangent Data Impact
In response to observations made in Sect. 4.1, it still appears relevant to sample fold axis data and incorporate them at least in terms of tangent data alone. 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. Twenty-nine initial tangent data points regularly spaced along the hinge line have been decreased to eight, and then four data points only.  (Fig. 6a), eight data points (Fig. 6c) or four data points (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 points does not seem to have a significant impact on the result. These different constraints lead to similar results (Fig. 6b, d, f) on the edge of the indicator's resolution. These statements are also similar to the ones observed in the previous experiment with derivative constraints (Sect. 4.1).

Impact of the Spatial Distribution of Tangent Data
Now considering a fixed number of tangent data points, 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 points, 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 bulging around the hinge line. The quantitative improvement with 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 (Sect. 2.2). -Scenario 2: The addition of tangent data does not significantly affect the resulting model regardless of the position of these data (Fig. 7b, d, f). RVI values are similar, within the limit of their uncertainties.
As a conclusion, with respect to only this 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 Sects. 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 Fig. 8 Comparison of the impact of 1 (first row), 2 (second row), 3 (third row) and 4 (last row) constraints per point at the same six locations on the fold axis of the theoretical surface 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 point 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 an understanding of its relevance. Second, the observations made in Sect. 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 sixfold 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 (Sects. 4.1, 4.2). However, the addition of the second tangent data (Fig. 8c) seems to lead to a flatter hinge line, but not especially better positioned. The addition of a 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 volume 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 points (Fig. 8b, d) 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) leads to a very satisfying model. Qualitatively, the resulting surface is very close to the theoretical one, closer than it has been at any time for the different attempts in this work. The 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 those previously made (Sect. 4.1: these data do not significantly affect 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 the results in a very significant way, at least as could be expected by users in such a situation.
Unfortunately, in practice it is not possible to directly sample these values on the field. The immediate integration of this constraint on real cases for now would there-fore 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 are 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 Sect. 2, using the potential method (Lajaunie et al. 1997) as a reference. This allows us 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 advantage 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 Sect. 3. It serves in Sect. 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 improves the results only slightly. 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 is bound to be limited in other geological studies as well. Indeed, one of the most important features characterizing complex isosurfaces 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 investigating 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 microscale 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 (Frodeman 1995). 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 three-dimensional curve sketch using the distance to a skeleton (Henrion et al. 2010;Garcia et al. 2018). 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.
We are particularly grateful to Christian Lajaunie, inventor of the geostatistical potential field method, who participated extensively 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. The authors thank Michael Hillier, Guillaume Caumon and the anonymous reviewer fo their fruitful comments.

Appendix A-Formulation of the Potential Field Method in Terms of Kriging Equations
Whereas the seminal paper by (Lajaunie et al. 1997) adopted a dual (co-)kriging formulation, we retained here the universal (co-)kriging formulation. The interested reader may consult chapter 3.4 of the book by (Chiles and Delfiner 2009) for greater detail about this derivation.
Equation 2 writes 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 ), 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 (A.10)

Appendix 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 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 which must be true for all values of ν and implies n equations (universality constraints), as many as there are drift functions denoted T (λ ii , λ g , λ s ) = 0 in the following.

Appendix A.2 Optimality Conditions
Optimality conditions ensure that Var( ) is minimal. In this part, the following convention is adopted. (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 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 should have its partial derivatives with respect to each λ ii , λ g , λ s and Lagrangian multipliers μ to be null. The Lagrangian L is then given by (A.17) Deriving this Lagrangian according to λ ii , λ g , λ s and μ gives which defines kriging equations and writes in matrix form Denoting N = (ii n + g n + s n ) as the number of data points and p, the number of points to estimate 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.
with the same abbreviation logic f g (p g ) = D τ g f (p g ) and f s (p s ) = D 2 τ s f (p s ).

Appendix B Appendix B-Covariance Function Derivation Developments Appendix B.1 Definitions
Covariance function of Z between point p and point q Stationarity case of order two assumption Case of constant anisotropic distance (global anisotropy) where A is a fixed rotation/dilatation matrix of R 3 . In the following, r (h pq ) may be denoted as r for readability.

Appendix B.3 First Derivative of the Covariance Function
Directional derivative at point p according to e x is The gradient of C Z with respect to p is First derivative specification of C Z with respect to p (B.33)

Appendix B.4 Second-Order Derivative of the Covariance Function
Covariance of partial derivatives w.r.t. p according to e x and w.r.t. q according to e y is Cov(D e x Z (p), D e x Z (q)) = t e x .D p,q C Z (h pq ).e y . (B.34) Second-order derivative specification of C Z with respect to p then q .35 leads to two 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.

Appendix 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 and 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)).

Appendix D Volume Error Indicator Appendix D.1 Volume 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 Sect. 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 = z − min(ẑ, z f loor ) 1 = |z(x, y) − min(ẑ(x, y), z f loor )|

Appendix 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 isosurface 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 volume 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. Appendix D.3 Relative Volume Improvement (RVI, ηV e ) S1 : Scenario 1; S2 : Scenario 2 V e : Volume error (u 3 ) V e 1 : Volume error of initial result of Scenario 1 (u 3 ) V e 2 : Volume error of initial result of Scenario 2 (u 3 ) where u is the arbitrary length unit used here. ηV e (RVI): Signed relative difference of V e compared to initial result in the same scenario (%) ηV e = V e i − V e V e i i ∈ (1, 2)