Ayoub Belhachmi
email: abelhachmi@slb.com

Azeddine Benabbou
email: abenabbou@slb.com

Bernard Mourrain
email: bernard.mourrain@inria.fr

A spline-based regularized method for the reconstruction of complex geological models

Keywords: Implicit modeling, Structural modeling, High thickness variations, Splines, Diffusion PDE, Interpolation

The study and exploration of the subsurface requires the construction of geological models. This task can be difficult, especially in complex geological settings, with various unconformities. These models are constructed from seismic or well data, which can be sparse and noisy.

In this paper, we propose a new method to compute a stratigraphic function that represents geological layers in arbitrary settings. This function interpolates the data using piecewise quadratic C 1 Powell-Sabin splines, and is regularized via a selfadaptive diffusion scheme. For the discretization, we use Powell-Sabin splines on triangular meshes.

.

The regularization is considered as the most challenging part of any implicit modeling approach. Often, existing regularization methods produce inconsistent geological models, in particular for data with high thickness variations. To handle this kind of data, we propose a new scheme in which a diffusion term is introduced and iteratively adapted to the shapes and variations in the data, while minimizing the interpolation error.

Introduction

Subsurface structures are usually described by geological models. Often, these structures can be very complex, with several geological faults and many stratigraphic layers. The geological model must faithfully represent the complexity observed in the data. To build these models, several approaches have been proposed in the past [START_REF] Mallet | Discrete modeling for natural objects[END_REF][START_REF] Sprague | Interpretive tools for 3-d structural geological modelling part ii: Surface design from sparse spatial data[END_REF][START_REF] Caumon | Surface-based 3d modeling of geological structures[END_REF]). The foundational approach in structural modeling, known as explicit structural modeling, consists of modeling the contact surfaces between geological units typically with parametric and/or polygonal surfaces (Mallet 1992, Sprague and[START_REF] Sprague | Interpretive tools for 3-d structural geological modelling part ii: Surface design from sparse spatial data[END_REF]). An alternative approach, referred to as implicit structural modeling, consists of representing the stratigraphic layers by an implicit stratigraphic function. More specifically, the data points are along interfaces delimiting stratigraphic layers, represented as a set of equipotential surfaces called horizons. Implicit structural modeling has been widely used to build geomodels from stratigraphic data [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Evans | 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][START_REF] Hillier | Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions[END_REF]). Among the implicit methods, we can distinguish two main classes of interpolation methods. The first class includes methods based on dual kriging or radial basis functions, where the interpolation is based on the data points locations [START_REF] Lajaunie | Foliation fields and 3d cartography in geology: Principles of a method based on potential interpolation[END_REF][START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3d from structural data: The potential-field method[END_REF]. The resulting linear system to solve is dense, and its size depends on the number of data points. The second class includes mesh-based methods, which rely on the discretization of the domain of study using a mesh and then the function is approximated using piecewise continuous polynomial basis functions supported around the nodes of the mesh, leading to a sparse linear system to solve whose size is depending on the number of mesh nodes. For example the discrete smooth interpolation methods (DSI) (See e.g. [START_REF] Mallet | Discrete modeling for natural objects[END_REF][START_REF] Evans | 3d-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF] and [START_REF] Irakarama | Finite element implicit 3d subsurface structural modeling[END_REF]). They discretize the stratigraphic function on a volumetric/surfacic mesh using finite element basis functions. Several modeling softwares have been developed in the oil and gas industry based on these methods, for example volume based modeling (VBM) by [START_REF] Souche | Construction of structurally and stratigraphically consistent structural models using the volume-based modelling technology: Applications to an australian dataset[END_REF].

In this paper, we present a new method to construct geological models from stratigraphic data. In our case, the stratigraphic layers are represented by an implicit spline function computed on a triangular mesh conforming to the geological faults, using a Powell-Sabin finite element scheme [START_REF] Speleers | Isogeometric analysis with powell-sabin splines for advection-diffusion-reaction problems[END_REF].

The stratigraphic data can be sparse and irregularly distributed over the domain. The locally supported basis functions, used in the finite element method, can lead to an underdetermined interpolation problem. Hence, a regularization term is introduced to enforce a regular behaviour of the implicit function all over the domain. In our specific context, we define a regular function as a function satisfying the following regularity properties: 1. Mean value property and extrema at the boundaries principle. A regular function must only contain extrema at the boundaries (e.g. interfaces of the geological layers), whilst mean values elsewhere, ensuring no local extrema. 2. Function must be smooth, through minimizing a conventional smoothness energy. 3. Function should adhere to input data particularly in the proximity of faults or boundary data, with only minimal regularization influence driving the function in these areas.

4. Consistent shape of the function far from data.

Smoothing techniques based on the minimization of second-order derivatives are extensively used as a regularization in the implicit modeling problem, see [START_REF] Irakarama | Finite element implicit 3d subsurface structural modeling[END_REF]. The Hessian smoothing energy is used by [START_REF] Renaudeau | Formulation continue du problème de modélisation implicite de structures géologiques discrétisée avec des méthodes de réduction de maillage[END_REF] as a regularization energy; Their approach involves constructing the implicit function using locally defined moving least square interpolants and explicitly minimizing the Hessian energy. The DSI methods consider many smoothing operators according to the choice of the basis functions used for interpolation. This includes the Hessian operator [START_REF] Irakarama | Finite element implicit 3d subsurface structural modeling[END_REF], smooth gradient [START_REF] Souche | Construction of structurally and stratigraphically consistent structural models using the volume-based modelling technology: Applications to an australian dataset[END_REF]), constant gradient (Frank et al. 2007) and Laplacian operator [START_REF] Lévy | Constrained discrete fairing for arbitrary meshes[END_REF]. In his pioneering work, [START_REF] Mallet | Discrete smooth interpolation[END_REF] introduced as a global roughness, a functional expressed using finite differences scheme involving the minimization of the first order partial derivatives as well as the second partial derivatives. [START_REF] Smith | Gridding with continuous curvature splines in tension[END_REF] introduced as regularization a partial differential equation combining the Laplacian and bilaplacian operators. For seismic data regularization, Fomel and Claerbout (2001) suggested the use of a diffusion tensor allowing for anisotropic smoothing in some predefined directions.

Existing regularization techniques can fail to reproduce geologically coherent solutions when dealing with models presenting high thickness variations, as it will be shown in this paper and also discussed in other works [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF][START_REF] Renaudeau | Formulation continue du problème de modélisation implicite de structures géologiques discrétisée avec des méthodes de réduction de maillage[END_REF][START_REF] Irakarama | Finite Difference Implicit Structural Modeling of Geological Structures[END_REF]). Geological models with strong thickness variations are models that represent geological structures with high variations in the shape and thickness of the layers. These variations are often anisotropic, meaning they vary differently along different directions. These strong variations have long been a concern and a challenge for different implicit modeling methods [START_REF] Laurent | Iterative thickness regularization of stratigraphic layers in discrete implicit modeling[END_REF][START_REF] Renaudeau | Formulation continue du problème de modélisation implicite de structures géologiques discrétisée avec des méthodes de réduction de maillage[END_REF][START_REF] Irakarama | Finite Difference Implicit Structural Modeling of Geological Structures[END_REF].

In this paper, we present a new regularization method based on diffusion energy. To better handle models presenting high thickness variations, we introduce a conductivity term in the diffusion equation that is iteratively adapted to the shape and the variation in the data. The diffusion equation provides solutions that respects the mean value property and the maxima at the boundaries. In this set of solutions, we seek the smoothest solutions, minimizing the Hessian energy. Furthermore, in order to ensure that the computed function is both smooth, and of a natural behaviour in the proximity of domain boundaries and faults, we use an approximation of the Neumann boundary condition. Further details can be found in Sect. 2.1.3. For the discretization, we use Powell-Sabin splines, as an alternative to the standard linear finite elements. They are quadratic splines with global C 1 continuity [START_REF] Dierckx | On calculating normalized powell-sabin b-splines[END_REF], defined on triangulation and their use for the discretization of second-order smoothness energies (e.g. Hessian energy) is straightforward. Our method produces solutions with higher regularity and smoothness, free from local extrema and boundary artifacts.

This paper is structured as follows. In Sect. 2, we present the mathematical framework of the geological modeling problem, and introduce the components of the regularization method we propose. In Sect. 3, we describe our method and provide details of the discretization using Powell-Sabin splines on a triangular mesh. In Sect. 4.5, we provide examples to demonstrate the limitations of classical smoothing techniques as regularization methods, particularly with the Hessian smoothing energy, when dealing with models that exhibit high thickness variations. We also provide examples that illustrate the effectiveness of our method in handling complex geological models that involve thickness variations and discontinuities. Finally, we conclude with a discussion on the shortcomings of the method, as well as the proposed solutions, and future perspectives.

Fitting and regularization

We consider a domain Ω ⊂ R 2 with N p stratigraphic data points and an implicit function u interpolating the data: (x p , y p , u(x p , y p )) for p = 1 : N p . Each subset of points with equal value represents a horizon. The values assigned to the horizons, have a strong impact on the resulting implicit functions [START_REF] Renaudeau | Formulation continue du problème de modélisation implicite de structures géologiques discrétisée avec des méthodes de réduction de maillage[END_REF]. Furthermore, it is common to conduct a prepossessing to determine optimal values that are adapted to the thickness variations of the layers [START_REF] Collon | 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France)[END_REF]. In our approach, we choose to assign arbitrary values, that increase in ascending order following the stratigraphic order of the layers. We also maintain a constant gap between the values assigned to two consecutive horizons. The domain Ω is discretized with a triangular mesh conformal to the discontinuities (geological faults). On the triangulation, we define basis functions in the space F of regular functions. The implicit function u is defined as

u(x) = N s ∑ i=1 B i (x)u i = B(x)U, ∀x ∈ Ω , (1)
where the row vector

B(x) = (B i (x)) N s i=1 is a basis of F, U = (u i ) N s i=1
is a column vector of unknowns and N s is the dimension of F. To obtain a solution of the interpolation problem of the geological data points over the domain Ω , we minimize the fitting energy E f it , with η p fitting weights associated to each of the data points.

E f it (u) = N p ∑ p=1 η p (u(x p , y p ) -z p ) 2 .
(2)

In our context, the data can be sparse and noisy, and basis functions are locally supported. This leads to an underdetermined system with an infinite number of possible solutions. To restrict the space of solutions a regularization term is introduced. The implicit modeling problem is posed as a minimization of a sum of energies

min u E f it (u) + λ E reg (u), (3)
where E f it is the fitting energy associated to the data constraints, E reg is the regularization energy associated to the regularity constraints, λ is the regularization weight controlling the tradeoff between the fitting and regularization constraints.

Regularization

The regularization of the implicit function is the main ingredient of any mesh-based method. Regularization techniques are used for selecting a specific set of functions, in the solution space, which satisfy the regularity constraints introduced beforehand.

Hessian smoothing and biharmonic equation

The Hessian energy also called squared Hessian [START_REF] Stein | Natural boundary conditions for smoothing in geometry processing[END_REF], roughness, bending or thin-plate energy, depending on the field of application, is a well-known regularization technique in data smoothing, image processing and many other engineering fields.

It is a smoothing energy based on the minimization over all the domain Ω of the squared partial second derivatives entries of the Hessian matrix,

H u ∈ R d×d . E H 2 (u) = Ω ∥H u ∥ 2 F = Ω (u 2 xx + 2u 2 xy + u 2 yy )dΩ . (4)
Minimizers of the Hessian energy are solutions of the biharmonic equation, see [START_REF] Stein | Natural boundary conditions for smoothing in geometry processing[END_REF])

∆ 2 u = 0 in Ω , (5)
where ∆ 2 is the bilaplacian operator. Furthermore, [START_REF] Stein | Natural boundary conditions for smoothing in geometry processing[END_REF] showed that the natural boundary conditions are

n T H u n = 0 on ∂ Ω , ∇∆ u • n + ∇(t T H u n) • t = 0 on ∂ Ω , (6)
n is the unit normal vector on the boundaries ∂ Ω and t the unit tangential vector along the boundaries ∂ Ω . They present higher smoothness and satisfy a specific natural boundary condition, studied and analyzed by [START_REF] Stein | Natural boundary conditions for smoothing in geometry processing[END_REF], forcing the function to be only linear in the normal direction to the boundaries. This behaviour is interesting in our application, since no explicit boundary condition is available.

Harmonic functions and diffusion equation

Harmonic functions are solutions of the standard Laplace equation

∆ u = 0 in Ω , ∇u • n = g on ∂ Ω . (7)
which describes the stationary state of the heat diffusion phenomenon, where g(x) is a known flux function on the boundaries. These functions exhibit properties that share a commonality with our approach, including maxima at the boundaries and mean value. The values of the implicit function are fixed for each horizon data points. A harmonic fitting funtion, will vary between consecutive horizons without any local extrema. However, the use of the diffusion equation to describe the stratigraphy can be a difficult task because the gradient of the implicit function is unknown on the boundaries of the domain. H 1 (Ω ) denote the Sobolev space of functions defined on Ω , function space that consists of functions with square-integrable derivatives up to the first order. Then, the weak form of Laplace equation using integration by parts is

Ω ∆ uv ds = 0 ⇔ Ω ∇u • ∇v ds - ∂ Ω (∇u • n)v dx = 0, ∀v ∈ H 1 (Ω ), (8)
⇔ Ω ∇u • ∇v ds - ∂ Ω gv dx = 0, ∀v ∈ H 1 (Ω ).
Solving Eq. ( 8) subject to ∇u • n = 0 is equivalent to the minimization of the Dirichlet energy [START_REF] Evans | 3d-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF])

E ∇ 2 (u) = 1 2 Ω ∥∇u∥ 2 dΩ . (9)

Handling the boundaries

When minimizing the Dirichlet energy Eq. ( 9) or solving standard Laplace equation Eq. ( 8), the natural boundary condition that emerges is the vanishing Neumann condition. It implicitly imposes that the normal derivative on the domain boundaries of the implicit function vanishes, which results in contours of the implicit function ending perpendicular to the domain boundaries to fulfill this condition. This is an undesired effect that makes the use of the diffusion equation unsuitable for our application. [START_REF] Irakarama | Finite element implicit 3d subsurface structural modeling[END_REF] proposed a free boundary discretization for the Laplacian avoiding this bias at the boundary. They use linear Lagrange elements as basis functions, and Crouzeix-Raviart nonconforming linear elements as test functions. Their method is based on a specific choice of basis functions and test functions spaces. In our work we avoid restricting ourselves to a specific choice of discretization. We propose a simple technique to avoid the implicit vanishing Neumann boundary condition by constructing an approximation of Neumann boundary condition instead of neglecting the integral of the normal derivative on the boundaries. The normal derivative on the boundaries, is the rate of change of the implicit function along the normal direction to the boundaries. It describes how the implicit function contours finish on the boundaries. To obtain an estimation of Neumann boundary condition, we propose constructing a first solution u 0 of the minimization problem Eq. (3) using the Hessian smoothing energy as regularization with a high regularization weight λ , in order to obtain a regular solution on the boundaries. We then use the normal derivative on the boundaries of the solution u 0 as a Neumann condition such that

Ω ∇u • ∇v ds = ∂ Ω (∇u 0 • n)v dx, ∀v ∈ H 1 (Ω ). (10)
Solving Eq. ( 10) is equivalent to minimizing the following diffusion energy [START_REF] Evans | 3d-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF])

E Di f (u) = 1 2 Ω ∥∇u∥ 2 dΩ - ∂ Ω gu dx with g = ∇u 0 • n on ∂ Ω . ( 11
)
Solution of the problem Eq. ( 3) with the diffusion energy Eq. ( 11) as regularization, inherits good behaviour near the domain boundaries and presents the regularity properties of the diffusion equation.

Isotropic and anisotropic diffusion

Anisotropic nonlinear diffusion is a powerful technique in image and geometry processing. It was introduced by [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] to smooth and denoise images while preserving sharp edges, based on the equation div(c∇u) = 0, with diffusion coefficient c = g(∇u) such as g(∇u) = 1 1+∥∇u∥ 2 or other functions inverting the norm of the gradient. This method is only isotropic, since it uses a scalar diffusion coefficient, and may not be sufficient to accurately describe models presenting anisotropic thickness variations.

The anisotropic diffusion process can be described using a variant of the heat diffusion partial differential equation (PDE)

∆ C = div(C∇u) = 0. ( 12
)
The anisotropic Laplacian ∆ C = div(C∇u) has spatially varying coefficients that weight the derivatives along different directions, specified by a diffusion tensor field C, such that C = (c i j ) with c ji = c i j for i, j = 1 : d. [START_REF] Andreux | Anisotropic laplace-beltrami operators for shape analysis[END_REF] introduced an anisotropic Laplace-Beltrami operator, for shape analysis, using a diagonal tensor encoding some geometrical extrinsic quantities. Expressed in the orthonormal basis of principal curvature directions, the diagonal entries of their tensor chosen as function of the principal curvatures, privileging diffusion in one of the principal curvature directions. The diffusion tensor C can be also defined

C = c ∥ b ⊗ b + c ⊥ (I d -b ⊗ b), (13)
where the vector b ∈ R d represents the anisotropy direction, c ∥ and c ⊥ are the parallel and perpendicular diffusion coefficients, respectively, and I d is the identity matrix. This representation is used in the context of magnetized plasmas in Tokamak [START_REF] Giorgiani | A high-order non field-aligned approach for the discretization of strongly anisotropic diffusion operators in magnetic fusion[END_REF], where the parallel diffusion c ∥ is chosen several orders of magnitude higher (up to 10 9 ) than the perpendicular diffusion coefficient k ⊥ , favoring the diffusion in the anisotropy direction. [START_REF] Wang | Fast quasi-harmonic weights for geometric data interpolation[END_REF] proposed optimization on the diffusion tensor to construct an optimal weights for geometric data interpolation, based on a positive semidefinite tensor representation where the coefficients of this tensor are unknowns of the problem.

In our context, models can exhibit high levels of anisotropy in which thickness properties can vary significantly along different directions. Moreover, the anisotropy directions are unknowns unlike the other fields of applications of the anisotropic diffusion equation. Multiple techniques exist for estimating global directions of anisotropy, however identifying local anisotropy directions remains a very challenging task. Inspired by the success of the anisotropic diffusion in all these application fields, we propose an iterative scheme in which the diffusion tensor is iteratively adapted to the anisotropy present in data. We consider the anisotropic Laplacian ∆ C defined by the symmetric diffusion tensor field C(x) : Ω -→ R d×d , where d is the dimension of Ω . For d = 2, we use c 1 , c 2 , c 3 scalar functions to denote the diffusion coefficients, respectively, in the directions: xx, xy, yy.

∆ C = div(C∇u) with C = c 1 c 2 c 2 c 3 . ( 14
)
The tensor representation of the diffusion generalizes the standard Laplacian for C = I d and the isotropic nonhomogeneous diffusion for C = c(x)I d , with c(x) a scalarvalue diffusion function. Solving the anistropic diffusion equation is equivalent to minimizing the following energy [START_REF] Evans | 3d-reconstruction of complex geological interfaces from irregularly distributed and noisy point data[END_REF])

E Di f (C) (u, g) = 1 2 Ω ∇u T (x)C(x)∇u(x) dΩ - ∂ Ω g(x)u(x)dx, (15)
where g(x) corresponds to the Neumann condition g(x) = ∇u 0 (x) T Cn for x ∈ ∂ Ω and n is the interior normal to ∂ Ω with u 0 a regular solution on the domain boundaries constructed using the technique described in Sect. 2.1.3.

A fitting method with diffusion regularization

In this section, we investigate the use of the anisotropic diffusion PDEs as regularization.

Then we propose a new formulation based on the nonlinear anisotropic diffusion with an iterative scheme in which the diffusion tensor is iteratively adapted to the variations and the anisotropy in data.

New formulation of regularization based on diffusion PDE

We consider the following nonlinear minimization problem with respect to the unknowns u and C

min u,C E f it (u) + λ E H 2 (u) + µE Di f (C) (u, g). (16)
The regularization term in our new formulation is a sum of two energies, weighted by two scalars µ and λ , controlling respectively the tradeoff between the diffusion and the smoothing. The diffusion energy is based on the anisotropic diffusion PDE, providing anisotropic solutions that respect the mean value property and the maximum principle, ensuring no local extrema. The Hessian energy term is added to ensure smoothness of the solution. In the set of solutions of the anisotropic diffusion equation with the diffusion tensor C, the smoothest solutions are selected by minimizing the Hessian energy. The diffusion energy defined with the anisotropic diffusion tensor C, generates solutions that are adapted to the anisotropy present in data, providing extra degrees of freedom for the modeling problem. Thus, the formulation of our problem is nonlinear regarding the two variables: the implicit function u and the diffusion tensor C.

To solve this nonlinear minimization problem, we propose an iterative scheme that involves the two problem unknowns u and C. In our scheme we solve first for u, then for C and we repeat until the difference between the fitting energies of consecutive iterations is less than a chosen threshold ε. The diffusion energy is quadratic in u and is convex if the tensor C is positive semi-definite. In this case, the solution is unique. However, the diffusion energy with respect to the coefficients of the diffusion tensor C is linear and can be either positive or negative. Thus, minimizing this energy with respect to C can be challenging as the solution is not well-defined. To address this issue, we introduce a proxy problem that replaces the original diffusion energy when optimizing for C. We propose to minimize instead a quadratic convex energy with respect to C. To do so, we first optimize the following objective function with respect to C for a given u min

C ( Ω ∥C∇u∥ 2 dΩ + ω 1 E H 2 (C) + ω 2 Ω ∥C -I d ∥ 2 F dΩ ). (17)
We then inject C solution of Eq. ( 17) and we optimize with respect to u

min u (E f it (u) + λ E H 2 (u) + µE Di f (C) (u, g)). (18)
Minimizers of the quadratic energy Ω ∥C∇u * ∥ 2 dΩ with respect to the components of the diffusion tensor C for a given u * , give minimal quantity C∇u * . Thus, the diffusion tensor C is adapted based on the intensity and the direction of the gradient ∇u * . In regions where the gradient ∇u * is strong, we minimize strongly the tensor C in the direction of the gradient ∇u * , allowing less diffusion regularization to permit strong variations of the implicit function. Conversely, where ∇u * is small, we minimize less the tensor C and therefore regularize more reducing the variations of the implicit function in the direction of ∇u * . The Hessian smoothing term ω 1 E H 2 (C) on the components of the diffusion tensors is added to obtain smooth diffusion coefficients with minimal oscillations. A third term ω 2 Ω ∥C -I d ∥ 2 F dΩ of a distance penalization from the standard diffusion where the diffusion tensor is the identity matrix I d , to avoid the trivial null tensor as solution for this problem.

To construct an initial solution u 0 , we only use the Hessian energy as regularization associated with a high regularization weight λ 0 , favoring regularity over a good fitting of data. C is then obtained via optimizing Eq. ( 17) with u = u 0 . Using C and g(x) = ∇u 0 (x) T Cn as Neumann boundary condition we obtain the next u. At step (i), C i is obtained via minimizing the Eq. ( 17) based on u i-1 . In the same way u i is obtained via minimizing Eq. ( 18) based on C i-1 and g(x) = ∇u i-1 (x) T C i n on the boundaries as Neumann boundary condition. We keep iterating until the difference between the fitting energies of the current and the previous iteration is less than a chosen threshold ε.

Algorithm 1 Iterative algorithm for nonlinear anisotropic diffusion energy minimization.

Require:

Initial solution u 0 ← min u E f it (u) + λ 0 E H 2 (u) for a high regularization weight λ 0 . while dif ≥ ε do C i ← min C E Di f (C, u i-1 ) + ω 1 E H 2 (C) + ω 2 Ω ∥C -I d ∥ 2 F dΩ ; g = ∇u T i-1 C i n; u i ← min u E f it (u) + λ E H 2 (u) + µE Di f (C i ) (u, g); di f = E f it (u i-1 ) -E f it (u i ); end while

Tangential diffusion along faults

Discontinuities such as faults are often encountered in geological structural modeling. The generated mesh is conformal to the faults in discrete implicit approaches. The nodes on the faults are duplicated from either sides allowing a jump in the implicit function. Faults are considered as inner boundaries of the domain Ω . We note ∂ Ω in ⊂ ∂ Ω the set of faults. Each fault is composed of a set of edges. In cases where the stratigraphic layers exhibit a large thickness variation along the faults, solutions of our method might not respect the mean value property and no local extrema along the faults. For this reason, we add an extra constraint to our regularization formulation Eq. ( 15) only on faults. Along the faults ∂ Ω in , in the tangential direction to the edges composing the faults, we apply the anisotropic diffusion PDE

(c t u t ) ′ = 0, (19)
where c t = t T Ct and u t is the derivative of the implicit function u in the direction t. The weak formulation is derived by multiplying Eq. ( 19) by a test function v and integrating over ∂ Ω in

∂ Ω in c t u t v t ds = [c t (x)u t (x)v(x)] x∈E , ∀v ∈ H 1 (Ω ), (20)
E is the set of nodes on the extremities of the faults. As illustrated in Fig. 1, we integrate Eq. ( 20) over the fault first from e 1 towards e 2 and then in the opposite direction from the other side. The resulting implicit boundary condition is a vanishing derivative in the direction t at the extremities of the faults. In order to avoid u t (x) = 0 for x ∈ E, we set u(x) = u i-1 (x) and ∇u(x) = ∇u i-1 (x) for x ∈ E as Dirichlet and Neumann boundary condition, leading to solving over all Ω the following PDE system

       div(C∇u) = 0 in Ω , ∇u T Cn = g on ∂ Ω , (c t u t ) ′ = 0 on ∂ Ω in , u(x) = u i-1 (x) , ∇u(x) = ∇u i-1 (x) f or x ∈ E. (21)
Equivalently we minimize the following modified diffusion energy

E Di f (C) (u, g) = 1 2 Ω ∇u T (x)C(x)∇u(x) dΩ - ∂ Ω g(x)u(x)dx + 1 2 ∂ Ω in c t (u t (x)) 2 dx (22) s.t u(x) = u i-1 (x), ∇u(x) = ∇u i-1 (x) f or x ∈ E.

Discretization

In this section, we discretize the energies used in the iterative scheme Eqs. ( 17)-( 18) as well as the fitting energy Eq. ( 2). Our formulation is generic (i.e. it works for any choice of basis functions (B i ) i=1:N s of F where one is able to discretize the diffusion and squared Hessian energy). The implicit function u is defined as

∀x ∈ Ω , u(x) = ∑ N s i=1 B i (x)u i = B(x)U, where B(x) = (B 1 (x), ..., B N s (x)
) is a basis of F and U = (u 1 , ..., u N s ) T is the vector of unknowns. Minimizing the fitting energy Eq. ( 2) is equivalent to solving the following linear system

N p ∑ p=1 η p B(x p , y p ) T B(x p , y p )U = N p ∑ p=1 η p B(x p , y p ) T z p . ( 23
)
We write this system as

A fit U = Z (24)
In order to discretize E H 2 (u) over the domain Ω , we consider the bilinear form h over

F h(u, v) = Ω (u xx v xx + 2u xy v xy + u yy v yy )dΩ . (25)
We define a sparse matrix

H ∈ R (N s ×N s ) E H 2 (u) = Ω ∥H u ∥ 2 F = U T HU s.t H i j = h(B i , B j ) for i, j = 1 : N s . ( 26
)
Minimizing E H 2 (u) gives the linear system HU = 0. We also express the diffusion energy using the basis functions as

E Di f (C) (u, g) = 1 2 Ω U T ∇B(x) T C(x)∇B(x)U dΩ - ∂ Ω g(x)B(x)Udx (27) + 1 2 ∂ Ω in U T (t T Ct)(B t (x) T (B t (x))Udx, where B t (x) = (∂ t (B i ))
is the derivative of the basis functions in the tangential direction to the faults. Therefore, minimizing Eq. ( 27) with respect to u, gives a linear system of the form ∆ C U = N C , where

∆ C ji = Ω C(x)∇B i • ∇B j dΩ + ∂ Ω in (t T C(x)t)∂ t (B i )∂ t (B j ) dx, i, j = 1 : N s , (28)
N C j = ∂ Ω g(x)B j dx. ( 29
)
To discretize the anisotropy symmetric tensor

C = c 1 c 2 c 2 c 3 , each diffusion coefficient c 1 , c 2 , c 3 is discretized using the basis functions (B i (x)) N s i=1 in F such that ∀x ∈ Ω c j (x) = N s ∑ i=1 B i (x)c j,i = B(x)C j , j = 1 : 3. ( 30
)
We use the flat vector

C =   C 1 C 2 C 3 
 to store the degrees of freedom of the three diffusion coefficients. The terms of the minimization problem (ii) are discretized as follows ∥C∇u∥ 2 = (C∇u) T (C∇u) for a given u.

(31)

C∇u = c 1 u x + c 2 u y c 2 u x + c 3 u y = u x u y 0 0 u x u y   c 1 c 2 c 3   = u x u y 0 0 u x u y ⊗ B(x) C. ( 32
)
Using the formulation above, we get

Ω ∥C∇u∥ 2 dΩ = Ω C T     u 2 x u x u y 0 u x u y u 2 x + u 2 y u x u y 0 u x u y u 2 y   ⊗ [B(x) T B(x)]   C dΩ . ( 33
)
The Hessian energy for the diffusion coefficients is simply

E H 2 (C) = E H 2 (c 1 ) + E H 2 (c 2 ) + E H 2 (c 3 ) = Ω C T (I 3 ⊗ H)C dΩ . (34)
The distance from the identity matrix explained in Sect. 3.1 can be expressed as

Ω ∥C -Id∥ 2 F dΩ = Ω (c 1 -1) 2 + (c 3 -1) 2 + 2c 2 2 dΩ . ( 35
)
Then we minimize it by solving the following linear system

    1 0 0 0 2 0 0 0 1   ⊗ M   C =   1 0 1   ⊗ Ω B(x) T dΩ , with M = Ω B(x) T B(x)dΩ . (36)

Powell-Sabin splines basis functions

We introduced the regularization formulation Eq. ( 15) and the fitting energy Eq. ( 2) without specifying a choice of basis functions for the space F. A variety of basis functions can be considered, without any changes in the regularization formulation.

In this paper, we use Powell-Sabin spline basis functions, referred to as PS-splines. We consider the use of PS-splines basis functions for the discretization of the iterative minimization problem Eqs. ( 17)-( 18) formulated in the previous section. We recall the main properties of these spline basis functions and we refer to [START_REF] Dierckx | On calculating normalized powell-sabin b-splines[END_REF] for the details of the construction and the proof of the properties. The considered PS-splines are piecewise quadratic splines with global C 1 continuity, defined on triangulation.

They are constructed on a specific PS-refinement (Fig. 2), have a local support, are convex and form a partition of unity. We use the normalized B-spline representation introduced by [START_REF] Dierckx | On calculating normalized powell-sabin b-splines[END_REF]. Additionally, PS-splines are based on a PS-6 split method that involves dividing each triangle of the input mesh into six quadratic Bezier triangles, while increasing the number of unknowns by only a factor of three (Fig. 2). The use of PS-splines presents many advantages. First, they are based on triangulations, which offer an advantage, over grids, when an adaptive refinement is necessary. The discretization of second order derivatives used in the Hessian energy minimization is straightforward compared to the classic linear discretization.

Linear system weighting

The minimization problem Eq. ( 18) gives rise to a linear system of the form

AU = B, ( 37
)
where the matrix A and the column vector B are expressed as

A = A f it + λ H + µ∆ C and B = Z + µN C . ( 38
)
The weights λ and µ associated, respectively, with the diffusion and Hessian energy, control the tradeoff between fitting the data and regularizing the solution, as well as the tradeoff between the diffusion properties and smoothness.

We use the weighting function m that gives the average of the absolute values of the matrix elements, defined from

R N s ×N s - → R as m(M) = 1 N 2 s N s ∑ i, j=1 |M i j |. ( 39
)
We apply this function to the matrices A f it , H, ∆ C resulting, respectively, from the fitting, Hessian and diffusion energies. The quantities

m(A f it ) m(H) and m(A f it )
m(∆ C ) estimate, respectively, the average relative weight of A f it to H and to ∆ C . We introduce two coefficients 0 ≤ α ≤ 1 and β ≥ 0 such that

   λ = (1 -α)β m(A f it ) m(H) , µ = αβ m(A f it ) m(∆ C ) . ( 40
)
The values of α and β define the solutions of Eq. ( 17): β is controlling the tradeoff between fitting and regularization, while α the balance between the diffusion and the Hessian smoothing effect. For α = 0.5 and β = 2, the three components have an equal contribution. A small value of β is privileging the fitting but can compromise the regularity of the surface mostly for models with high thickness variations. Taking 0 ≤ α ≤ 1 ensures staying in between the two extrema: the regularization is complete smoothing for α = 0 and only a diffusion for α = 1.

Applications

Classical smoothing approaches based on the minimization of the second derivatives are widely used as regularization for the implicit modeling problem Eq. ( 3). However, these methods might fail to produce regular solutions when used on models exhibiting high thickness variation. In this section, we consider the Hessian smoothing energy as regularization and we show resulting solutions on models presenting high thickness variation. Figures 3a,b demonstrates the interpolation of sparse data in one dimension. The data points are irregularly distributed and exhibit high variations in the implicit function values assigned to each data point. The resulting implicit functions are shown in Fig. 3 as a multicolored line, with the first derivative of the implicit function as a blue dashed line and the second derivative as a red line. In one dimension, the Hessian is simply Ω (u xx ) 2 dx. In Fig. 3, the univariate normalized PS-splines [START_REF] Speleers | Multivariate normalized powell-sabin b-splines and quasi-interpolants[END_REF]) are used for interpolation. They are piecewise quadratic basis functions, which explains why the second derivative of the implicit function is piecewise constant.

When minimizing the second derivatives, the resulting function tends to have a smooth first derivative, which means that it cannot have strong changes. As a consequence of this, the function is more exposed to local extrema, which are points where the first derivative changes sign (see Fig. 3a). Our method minimizes globally the first and the second derivatives but tolerates a high jump in the second and first derivatives in the region between the two plateaus; see Fig. 3b. To benchmark the method in two dimensions, we use synthetic models with high thickness variation introduced in the data. We refer to these models as bell, rings, extracted layers and faulted synthetic, respectively, as shown in Figs. 4a,b,c,d. The thickness variation is present mainly in the y-direction in models Figs. 4c,d, while, in Fig. 4b along the radial direction and Fig. 4a along x andy. Data points along the rings have the same isovalue. The three rings have, respectively, 1, 2, 10 as values.

The bell model has horizons with values of -1, 0, 1, 2, and 3 from bottom to top and is a benchmark model of [START_REF] Renaudeau | Formulation continue du problème de modélisation implicite de structures géologiques discrétisée avec des méthodes de réduction de maillage[END_REF]. The faulted synthetic model has three horizons with values of 0, 1, and 2 from bottom to top, while the extracted layers model has horizons with values of 6, 7, 8, and 9 from bottom to top. Hessian smoothing on these models, as shown in Fig. 4, produces irregular solutions, where the mean value property and the maximum principle are violated with the presence of bubbles that are oscillations around local extrema in the resulting implicit function. These bubbles represent an undesired artifact as they are geologically incoherent, and occur in the regions where the variation of thickness in the data is the strongest. In Fig. 4a, bubbles are present on the boundaries and around the middle curved horizon, therefore the extraction in red of the isovalue corresponding to this horizon is geologically inconsistent. In Figs. 4b andc, the isovalues 1 and 9 appear twice and are highlighted in red. Solutions obtained via the minimization of the Hessian energy are really smooth, with desired behaviour near the domain boundaries, except the boundaries where a high thickness variation is present. These solutions can be seen to follow the data trend, with minimal bias near boundaries (Fig. 4d left side). However, in the presence of high thickness variations, they are highly oscillative (Fig. 4d right side). Despite measuring minimal Hessian energy, these solutions are geologically inconsistent.

Handling the boundaries

The use of the diffusion equation implies implicitly a vanishing Neumann boundary condition, leading to contours finishing perpendicular to the boundaries, as shown in Fig. 5a. This strong artifact is undesired for our application due to the presence of discontinuities inside the domain in Fig. 5a. To overcome this limitation, as described in Sect. 2.1.3, we first construct a function obtained using the Hessian smoothing with a high weight in Fig. 5b, giving a regular function inside the domain and naturally finishing on the boundaries. We use its normal component of the gradient on the boundaries as Neumann initial boundary condition and to compute the first diffusion tensor by minimizing Eq. ( 17). For the next iterations, we use the normal component of the gradient of the previous solution to update the Neumann boundary condition.

Our solution finally respects all the regularity criteria and fitting constraints with a natural behaviour near boundaries, free from the artifact as shown in Fig. 5c. In this section, we present a comparison between three diffusion schemes. Standard diffusion with the diffusion tensor set as the identity matrix, isotropic diffusion as introduced in Sect. 2.2 and anisotropic diffusion as introduced in Eq. ( 15). We consider the iterative scheme on u and C, respectively, Eqs. ( 18) and ( 17 The three diffusion schemes produce regular and globally smooth solutions, respecting the mean value property and the maximum principle at the boundaries. In the first model, where the thickness variation is only along the y-direction, almost no difference is visually observed between the three different solutions. However, significant differences are observed in the regions with anisotropic high thickness variations: the flat region under the curved horizon in the bell model and between the rings with isovalue 1 and 2 in the rings model. In these regions the standard diffusion fails to produce smooth curves and some undesired oscillations are observed. The isotropic diffusion produces fewer oscillations and the curves are smoother than the standard diffusion. The anisotropic diffusion using the diffusion tensor is visually producing solutions with higher smoothness and respecting the regularity criteria, with a perfect fitting as well of the curved horizon in the middle highlighted in red; see Fig. 6i. Decreasing the regularization weight improves the fitting of this horizon in the standard and isotropic case, as shown in Figs. 6g,h, respectively; however, it introduces more oscillations.

Standard

Isotropic Anisotropic

(a) (b) (c) (d) (e) (f) (g) (h) (i)
Fig. 6 Resulting implicit functions using different diffusion schemes in high thickness variations models

Mesh refinement effect

The level of refinement of the mesh determines the quality of fitting and the regularity of the solutions of the minimization problem (Eq. 16). Our method, shown in Figs. 6c, f, i, produces successfully smooth solution that respects all the regularity criteria and accurately fits the data. This is achieved on meshes with a coarse level of refinement as shown in the previous comparison section. The three examples shown in Fig. 7 demonstrate that refining the mesh leads to an enhancement in the quality of data fitting for the three diffusion schemes (Fig. 7). Nevertheless, for the standard diffusion (Fig. 7a) and isotropic diffusion (Fig. 7b), some oscillations can be seen in the flat region beneath the top of the curved horizon. Meanwhile, the anisotropic diffusion (Fig. 7c) provides a solution with higher smoothness and regularity and the no local extrema property is holding. The main difference between the solutions of the three diffusion schemes lies in the smoothness of the curves in the flat region, where the thickness variation is strong and the implicit function varies weakly. In the case of interest for this region, a high refinement is necessary to maintain the good properties of the solutions of our method. Advanced refinement techniques, such as adaptative mesh refinement, can be an effective way to improve the quality of solutions. Adaptive refinement techniques allow the mesh to be more refined in areas where the high thickness variations are present and in the anisotropy directions. We refer to the model in Fig. 7 as bell refined. In this section, we show the effectiveness of our method on complex geological models that contain faults. We first evaluate our method on a faulted model presenting no significant thickness variations (Fig. 8), and another presenting strong thickness variations (Fig. 9). The faulted model in Fig. 8 contains seven horizons. The values assigned to the horizons are in the ascending order from the bottom to the top, ranging from one to seven. The thickness variation in this model is not significant/pronounced, since the thickness of the layers between these horizons are relatively uniform. Figure 8 shows resulting implicit functions using the Hessian smoothing and our method as regularization; Figs. 8a andb. Our method produces a regular solution throughout the entire domain, displaying geologically consistent behaviour near the boundaries and faults. Similarly, the Hessian smoothing solution respects all the regularity criteria.

Visually, the two solutions are quite similar, except near the fault inside the domain due to the use of the tangential diffusion constraint introduced on the faults for our method. The model in Fig. 9 has been obtained by removing horizons 4, 5 and 7. To further complicate the problem, the value associated with the horizon 6 has been changed to 4. These modifications induce a very strong thickness variations, particularly along the fault inside the domain. We refer to the models Fig. 8 and Fig. 9 respectively as Geo and Geo refined.

The Hessian smoothing fails to produce a regular solution (Fig. 9a). The maximum principle and the mean value property are violated throughout the domain. The resulting implicit function increases until it reaches the value 4, and then decreases, which induces a wrong extraction of the horizon 4, being extracted twice.

The method Eq. ( 15) with no tangential diffusion along the faults, generates a regular solution within the domain. However, along the faults it is not maintaining the mean value and the maximum principle, resulting in the formation of bubbles. To overcome this limitation, we introduced the tangential diffusion along the faults in the diffusion energy (Eq. 22). It ensures the respect of the regularity criteria along the faults, thereby eliminates the bubbles as shown in Fig. 9c. The average fitting error in Table . (1) is computed by dividing the fitting energy Eq. ( 2) by the number of data points, and all fitting weights are equal to 1. For all the models, we converge after three or four iterations. In our experiments, we opt for a parameter α = 0.5 for the models exhibiting moderate thickness variations and 0.9 for models with strong thickness variations. The parameter β is typically set equal to one, while we choose small values for the models with a sparse density of data points and a very refined mesh.

To determine the high regularization weight denoted λ 0 in Algorithm. (1), which is associated to the Hessian energy to construct the initial solution u 0 , we set

λ 0 = p m(A f it ) m(H) , (41)
with p a positive scalar. For all the models, a p = 10, gives a regular solution. For Geo and Geo refined models, we choose a p = 100. While for Bell and Bell refined models, we choose p = 1000. Finally, we choose p = 50 for the Extracted layers model.

Model Number of iterations Average fitting error α β

Extracted layers (Fig. 6d) 4 6.567e -3 0.5 1 Faulted synthetic (Fig. 5d) 4 0.02385 0.5 1 Bell refined (Fig. 7) 3 2.3116e -4 0.9 1 Bell (Fig. 6i) 4 3.3e -3 0.9 0.1 Rings (Fig. 6f) 3 4e -3 0.9 1 Geo (Fig. 8b) 4 1.8e -3 0.5 1 Geo refined (Fig. 9c) 4 6.391e -6 0.9 0.001

Table 1 Number of iterations, the fitting error and the parametrs (α,β ) per model.

Discussion

Our method depends on the initial solution u 0 obtained using the Hessian smoothing energy. In cases where strong curvatures are present in data, starting from a very penalized solution u 0 (i.e. obtained using a high smoothing weight) can result in a final solution where some regions are flat. To address this, a less penalized solution u 0 should be privileged.

Our approach is based on approximating the Neumann boundary condition to prevent the issue of vanishing implicit Neumann boundary condition. We first consider as an approximation the normal component on the boundaries of the gradient of the initial solution u 0 . During the iterative process, we use the normal component of the gradient of the previous solution as Neumann boundary condition. In areas near faults and domain boundaries where data are unavailable, the final solution will tend to have the same behaviour of the initial solution.

The minimization of Eq. ( 17) for tensorC aims to minimize the diffusion coefficients (c i ) i=1:3 , while keeping c 1 and c 3 within the range of 0 to 1. However, this approach does not guarantee that the resulting tensor C will be positive semidefinite. Further investigations are needed to ensure that the tensors remain semidefinite positive, and to study the undesired consequences (e.g. oscillations or numerical instability) of any violations of this condition on the resulting solutions.

The solution U of Eq. 37 is used to evaluate the implicit function all over the domain. In models without faults, we discretize the domain using a regular grid, and the implicit function is evaluated on the grid's nodes. The isovalues are then extracted linearly on each grid element using marching squares. However, in presence of discontinuities, we refine the mesh generated conformal to the faults for interpolation, and we use a marching triangles method for the extraction of the isovalues. It is worth nothing, that the evaluation mesh used for the extraction of the isovalues, is refined up to three times compared to the mesh used for interpolation.

Conclusion and perspectives

In this paper, we introduce a new regularization method based on autoadaptive anisotropic diffusion. Our method is an iterative scheme, where a diffusion term is adapted to the variations present in the data. Contrary to existing methods, our approach is data driven and no preprocessing is needed to assign values to geological horizons. Instead, our iterative scheme adapt the diffusion tensor to the anisotropy and thickness variations present in the data to obtain a regular solution. Furthermore, in our iterative scheme, the estimation of the Neumann boundary conditions is performed simply by integrating the normal derivative of the previous solution on the boundary.

Extending our iterative scheme to three dimensions using Powell-Sabin splines is feasible. [START_REF] Speleers | Multivariate normalized powell-sabin b-splines and quasi-interpolants[END_REF] presented a method for constructing a normalized basis for the multivariate quadratic spline space defined over a generalized Powell-Sabin refinement of a triangulation in R S , s ≥ 1. The extension of our regularization formulation to three dimensions with the anisotropic diffusion tensor is possible, as well as to other discretizations as long as one is able to discretize the diffusion and the Hessian energy. However, it is worth noting that we did not conduct experiments in the volumetric context. [START_REF] Irakarama | Finite element implicit 3d subsurface structural modeling[END_REF] proposed a free boundary discretization for the Laplacian, by using a specific choice of the space of basis functions and test function space. A similar approach can be applied to obtain a free boundary discretization for the anisotropic Laplacian.

Wang and Solomon (2021) introduced a Cholesky factorization of the diffusion tensor and considered an optimization on the elements of this factorized matrix instead of the diffusion tensor elements. The same parameterization can be used in our method to ensure that tensors stay positive semidefinite.

When constructing geological models, it is common to encounter situations where new data become available after the initial model has been computed. In such cases, a local update technique can be employed to enrich specific regions in the model using the new data, without making significant changes to the rest of the model. In future works, we will explore the use of local update techniques in such cases.

Fig. 1

 1 Fig. 1 Fault inside the mesh with extremities e 1 , e 2 . Red line represents the fault, blue arrows represent the direction of the fault edges in both sides.

Fig. 2

 2 Fig. 2 PS-refinement of a triangular mesh, (a) input, (b) refined mesh

Fig. 3

 3 Fig. 3 One dimensional interpolation using (a) Hessian energy (b) our method. Black dots are data points to interpolate, red line in the graph corresponds to the second derivative of the implicit function, while the blue dashed line represents the first derivative.

Fig. 4

 4 Fig. 4 Resulting implicit functions using the Hessian smoothing on models presenting high thickness variations and discontinuities. (a) Bell. (b) Rings. (c) Extracted layers. (d) Faulted synthetic. Black dots are data points ordered along horizons. Red contours represent the extracted isovalues that are geologically inconsistent, while green contours represent extracted isovalues that are geologically valid .

Fig. 5

 5 Fig. 5 (a) Resulting implicit function with implicit vanishing Neumann condition on the boundaries. (b) Resulting implicit function with a high smoothing weight using the Hessian energy. (c) Resulting implicit function free from boundary artifact using our method

  ) for the isotropic and the anisotropic diffusion. The final solutions obtained are shown for comparison. For comparison, we create models with thickness variations introduced in the data along one direction (y direction) in Figs.6a, b, c first model and along two directions in the second and third, respectively, inFigs. 6d, e, f and Figs. 6g, h, i. In Figs. 6d, e, f data points along the rings have the same isovalue, and the three rings have 1, 2, 10 as values. In Figs.6g, h, i the values on the horizons are: -1, 0, 1, 2, 3, respectively, from the bottom to the top.

Fig. 7

 7 Fig. 7 Resulting implicit function on a refined mesh using (a) standard diffusion, (b) isotropic diffusion, (c) anisotropic diffusion

Fig. 8

 8 Fig. 8 Geological benchmark model of Renaudeau (2019). (a) Resulting implicit function using the Hessian energy, (b) resulting implicit function using our method. (c) A zoomed-in view of the area within the black frame in Figure (a), (d) A zoomed-in view of the area within the black frame in Figure (b).

Fig. 9

 9 Fig.9(a) Resulting implicit function using the Hessian energy, (b) resulting implicit function using our method without the tangential diffusion constraint, (c) resulting implicit function using our method with the tangential diffusion constraint.

Acknowledgments The authors would like to thank Boniface Nkonga (Université Côte d'Azur and INRIA) and Angelos Mantzaflaris (Université Côte d'Azur and INRIA) for their insightful discussions on the anisotropic diffusion equation. The

authors also thank Thomas Laverne, Emmanuel Malvesin and Julien Renaudeau (all from SLB) for valuable discussions on implicit structural modeling, as well as for providing the synthetic geological model used in this work.

The Ph.D. thesis of Ayoub Belhachmi is funded by both SLB and the Association nationale recherche technologie.

Conflict of interest The authors declare that they have no conflict of interest.