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 self-adaptive diffusion scheme. For the discretization, we use Powell–Sabin splines on triangular meshes. Compared to classical interpolation methods, the use of piecewise quadratic splines has two major advantages. First, they have the ability to produce surfaces of higher smoothness and regularity. Second, it is straightforward to discretize high-order smoothness energies like the squared Hessian energy. 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.
