AdaHRBF v1.0: Gradient-Adaptive Hermite-Birkhoff Radial Basis Function Interpolants for Three-dimensional Stratigraphic Implicit Modeling

. Three-dimensional (3D) stratigraphic modelling is capable of modeling the shape, topology, and other properties of strata in a digitalized manner. The implicit modeling approach is becoming the mainstream approach for 3D stratigraphic modelling, which incorporates both the off-contact attitudes and the on-contact occurrence information of stratigraphic interface to estimate the stratigraphic potential field (SPF) to represent the 3D architectures of strata. However, the magnitudes 15 of SPF gradient controlling variation trend of SPF values cannot be directly derived from the known stratigraphic attribute or attitude data. In this paper, we propose an Hermite-Birkhoff radial basis function (HRBF) formulation, AdaHRBF, with an adaptive gradient magnitude for continuous 3D SPF modeling of multiple stratigraphic interfaces. In the linear system of HRBF interpolant constrained by the scattered on-contact attribute points and off-contact attitude points of a set of strata in 3D space, we add a novel optimizing term to iteratively obtain the true gradient magnitude. The case study shows that the 20 HRBF interpolants can consistently establish accurate multiple stratigraphic interfaces and fully express the internal stratigraphic attribute and attitude. To ensure harmony of the variation of stratigraphic thickness, we adopt the relative burial depth of stratigraphic interface to the Quaternary as the SPF attribute value and propose a new stratigraphical thickness index (STI) to represent the variation trend of stratigraphic thickness in SPF. In addition, the proposed stratigraphic potential field modeling by HRBF interpolants can provide a suitable basic model for subsequent geosciences numerical simulation

stratigraphic potential field modeling using the AdaHRBF method. The SPF attribute value is set to the relative burial depth of strata, i.e., mean distance from a given stratigraphic surface to the top surface of the Quaternary, meanwhile, we propose a new stratigraphic thickness index (STI) to express the variation trend of stratigraphic thickness. The distributions of attribute, thickness, and attitude of strata in 3D space can be fully expressed by the SPF and its gradient vector field.

65
The key of implicit modeling methods is to interpolate a 3D scalar field function whose equipotential surfaces indicate the boundaries of geological bodies. These surfaces can represent ore grade boundaries or stratigraphic interfaces. This scalar field is interpolated from stratigraphic interface points and attitude data with either discrete interpolation schemes or continuous interpolation schemes.
For discrete interpolation schemes of implicit modeling with a special mesh, Mallet (1989); Mallet (1992) proposed a discrete 70 smooth interpolation (DSI) method of producing values only at the mesh points on the stratigraphic interface instead of explicitly computing a function defined everywhere. The GoCAD (www.pdgm.com/products/skua-gocad/) software was developed based on the DSI method to meet the needs of geological, geophysical, and petroleum reservoir engineering modeling (Collon et al., 2015;Zhang et al., 2020). Caumon et al. (2013) proposed a discretizing finite-element method (FEM) to generate 3D models of horizons on a tetrahedral mesh, using stratigraphic interface traces of unknown attribute 75 values and attitude measurements from 2D geological maps, remote sensing images, and digital elevation models. Hillier et al. (2013) presented a structural field interpolation (SFI) algorithm using an anisotropic inverse distance weighted (IDW) interpolation scheme derived from eigen analysis of strike/dip measurements. Gonçalves et al. (2017) proposed a vector potential-field solution from a machine learning perspective, recasting the problem as multivariate classification in a compositional data framework, which alleviates some of the assumptions of the cokriging method.

80
Since the continuous interpolation schemes does not depend on a mesh for its definition, the stratigraphic interfaces can be extracted at any desired resolution in the specific volume of interest. There is already a dual kriging or cokriging formulation for continuous potential field modeling of multiple stratigraphic interfaces. Lajaunie et al. (1997) proposed an implicit potential field modeling method using the dual formulation of kriging interpolation that considers known points on a geological interface and plane attitude data such as stratification or foliation planes. Calcagno et al. (2008) cokriged the location of 85 geological interfaces and attitude data from a structural field to interpolate a continuous 3D potential-field scalar function describing the geometry of geological bodies. Geomodeller 3D (www.geomodeller.com), an implicit geological modeling application, utilizes the implicit potential field method by cokriging or the dual formulation of kriging (Lindsay et al., 2012; implementation, to generate 3D geological models based on an implicit potential-field cokriging interpolation approach and to enable stochastic geological modeling and inversions of gravity and topology in machine-learning and Bayesian inference frameworks. Renaudeau et al. (2019) proposed an implicit structural modeling method using locally defined moving least squares shape functions and solved a sparse sampling problem without relying on a complex mesh. Irakarama et al. (2020) introduced a new method for implicit structural modeling by regularization operators using finite differences. To reduce the impact of regularly occurring modeling artifacts that result from data configuration and uncertainty, Von Harten et al. (2021)  95 proposed an approach that is a combination of an implicit interpolation algorithm with a local smoothing method based on the concepts of nugget effect and filtered kriging known from conventional geostatistics.
For continuous radial basis function (RBF) or HRBF interpolation schemes of implicit modeling without a mesh, Cowan et al. (2003) constructed an implicit model of the orebody or stratigraphic interface using a volumetric RBF interpolation function with an equipotential surface that includes the interface points, and conventionally assigned an attribute value of zero and a 100 "±" sign to indicate the inside and outside of the interface. Hillier et al. (2014) presented a generalized interpolation framework using RBF to implicitly model 3D continuous geological interfaces from on-surface points with gradient constraints as defined by strike-dip data with assigned polarity. Leapfrog Geo (www.leapfrog3d.com) is an implicit geological modeling software package that models scattered data for interfaces using fast RBF interpolation methods (Vollgger et al., 2015;Stoch et al., 2020;Creus et al., 2018;Basson et al., 2016;Basson et al., 2017). Martin and Boisvert (2017) developed a RBF-based 105 implicit modeling framework using domain decomposition to locally vary orientations and magnitudes of anisotropy for geological boundary models. Zhong et al. (2019); Zhong et al. (2021) introduced combination constraints for modeling ore bodies based on multiple implicit fields interpolation through RBF methods, in which a multiply labeled implicit function was defined that combines different implicit sub-fields by the combination operations to construct constraints honoring geological relationships more flexibly. Guo et al. (2016); Guo et al. (2018); Guo et al. (2020); Guo et al. (2021) proposed an explicit-110 implicit-integrated 3D geological modelling approach for the geometric fusion of different types of complex geological structure models; therein, the HRBF-based implicit method was used to model general strata, faults, and folds, and the skinning method and the free-form surface were used to model local detailed structures. Wang et al. (2018) proposed an implicit modeling approach to automatically build a 3D model for orebodies by means of spatial HRBF interpolation directly based on geological borehole data.

115
However, the above RBF or HRBF interpolants, which use only the on-contact point datasets for each geological interface or assign an approximate gradient vector for each on-contact point according to its nearest attitude measurements, cannot be accurately consistent with actual attitude survey data. Moreover, RBF/HRBF-based methods construct implicit field functions separately for each geological interface and extract the zero value equipotential surfaces to locate the geological interface.
Therefore, it is difficult to maintain topological consistency between geological bodies, let alone to represent their internal 120 attributes and structural attitudes. Our AdaHRBF interpolation scheme yields an HRBF linear system that is analogous in form https://doi.org/10.5194/egusphere-2022-1304 Preprint. Discussion started: 23 January 2023 c Author(s) 2023. CC BY 4.0 License.
to the previously developed implicit potential field interpolation method based on cokriging of contact increments using parametric isotropic covariance functions.

125
The geological boundaries and structural attitudes on geological maps and cross-sections are the most common data used for 3D geological modeling. Besides the geological boundaries extracted from boreholes, cross-sections, and geological maps, structural attitude (including strike direction, dip direction, and dip angle) data from geological maps play important roles in characterizing the shape and distribution of geological bodies. As shown in Fig. 1a, stratum S1 is between its bottom surface f1 and top surface f2; a fault interface F divides the 3D space into two sub-domains D1 and D2. We can extract the on-contact 130 boundary points and off-contact attitude points of the strata and fault from the cross-section AA' (Fig. 1b) and geological map ( Fig. 1c). The SPF modeling method can jointly reconstruct a 3D geological model using these data extracted from geological maps and cross-sections. A field in a spatial domain n defines the function f=f(p) at a point p∈ n in domain n , and f(p) is also called field function.
The SPF defines the 3D space as a scalar function f(p) at any point p, meanwhile, the stratigraphic interfaces are simulated and expressed as specific equipotential surfaces satisfying f(p)=fk (i = 1, ..., K) in the SPF. In practice, this specific function value fk may correspond to the age of the stratigraphic interface or a relative distance from a reference interface (Mallet, 2004). 140 Therefore, a stratum occupies the space between its bottom surface fk and top surface fk+1, while there are countless disjoint equipotential surfaces in each stratum (Mallet, 2004). A well-known problem is how to interpolate unknown points by a function f(·) using known points of the space n . The key problem of SPF modelling is to obtain surfaces that are consistent with known on-contact points on the stratigraphic interfaces and the off-contact attitudes of the strata. The stratigraphic interface points define the distribution of reference equipotential surfaces, while the attitude points define the gradient vectors 145 of the scalar field.
The SPF modeling by the HRBF interpolant satisfies both the on-contact attribute constraint and off-contact attitude constraint.

HRBF Interpolant
With the joint constraints of ( ) = and ( ) = , the optional solution is to obtain equipotential surfaces that are as smooth as possible, that is, to ensure the energy function of SPF, which represents the degree of equipotential surface smoothness and unevenness, as small as possible (Bachman and Narici, 2000). Therefore, the energy function (E) of the SPF 155 is defined by the second-order derivative of ( ) as: where, ‖ − ‖ denotes the Euclidean distance between locations p and ; ( ) is the radial basis function, herein, for which the cubic function ( ) = 3 was used in this study; is the Hamiltonian operator; is the Hessian operator, in Furthermore, the type of basic functions (e.g., Gaussian, multi-quadric, and thin plate spline, as shown in Table 1) affects the result of spatial interpolation. We adopt the cubic function as the basis function in this study, i.e., φ(r) = 3 .  According to the joint constraints, the weight coefficients , , and of the interpolant are determined by the following linear system: , whose elements = [1 0 0] , = [0 1 0] , and = [0 0 1] , respectively.
Once we have the weight coefficients α i , , and the polynomial coefficients (c1, c2, c3, c4) by solving the above HRBF

190
linear system, we can substitute the weight coefficients and polynomial coefficients into the HRBF equations, then the interpolant function f(p) and its gradient function ( ) can be easily obtained.

Determination of Gradient Direction
The gradient of the SPF is an important feature of stratum shape, because it indicates the attitude of a stratum. For construction 195 of a scalar field ( ), the gradient constraints ( ) = can also be added into modeling process. As shown in Fig  The gradient is a vector with magnitude and direction (which is the same as the normal direction of the stratigraphic interface). The X-axis, Y-axis, and Z-axis components of the normal direction, , , and , in the 3D Cartesian coordinate system can be derived from the strike, dip and angle of dip of the stratigraphic interface as following:

Optimization of Gradient Magnitude
However, it is difficult to obtain the gradient magnitude through any geological observation. The exact definition of gradient magnitude (‖ ‖) is the change of an attribute value over unit distance along the gradient direction. The gradient magnitude reflects the rate of change of the scalar field values, which is caused by the difference of stratum thickness at different locations.
As shown in Fig. 3a, a larger gradient magnitude 1 indicates that the stratum becomes thinner, whereas a smaller gradient 215 magnitude 3 indicates that the stratum tends to become thicker. We assume that the gradient magnitude changes gradually everywhere in the scalar field; therefore, every equipotential surface inside of the stratum changes uniformly. However, as shown in Fig. 3b, if we force all the gradient magnitudes to be equal, it may cause the inconsistent SPF changes with neighbors; that is, the trends of some equipotential surfaces inside of the stratum change suddenly compared to other equipotential surfaces. constraint. When →0, the solution is close to interpolation. [ where the diagonal coefficient matrix is given by = ( ]. Given the 230 gradient magnitudes = |l 1 l 2 ⋯ l M | T , then the gradient = l × is the product of lk and normal vector n k : We initially set (t=0) to a nonzero constant vector and ( =0) = 1. After solving the HRBF system, we can get the function of scalar field ( ), then the gradient vector on the attitude observed point is easily obtained according to = ∇ ( ).
We record the HRBF coefficients calculated at the t-th time as α and , and record the gradient magnitude at the attitude 235 observed point as . After solution of the linear system in Equation (6), we estimate the gradient magnitudes and the gradient constraint at next iteration step as t = × . Accordingly, we shrink the coefficient +1 to fit more closely to the update gradient constraint. Our idea is that when gradient magnitudes converge, the resulting implicit function interpolates the converged .
In this study, we calculate the increment of λ from: where UQ() is the upper quartile of differences of all gradient magnitudes. Given the updated +1 and , we substitute them into the (t+1)-th HRBF system (Eq. 6) and solve for the updated coefficient of implicit function. This iterative process continues until the stopping criteria is satisfied.
We use two stopping criteria to finish the iterations. Firstly, for all observed attitude points, if the sum of differences of gradient 245 magnitudes between two consecutive iterations is less than or equal to a small enough threshold ε, we stop the iterations on convergency. Secondly, the number of iterations reaches a given number Niterate, we also obtain the final results of α , , and where | | represents the absolute value of a real number and M is the number of observed attitude points. The basic steps of the iterative calculation of gradient magnitude are given in the pseudo code (Fig. 4).

4.
Calculate known points by , and .

Verification Experiments
Two experimental fields in 2D space, with gradient changing in direction or magnitude, were designed to verify the AdaHRBF 255 method. The experimental results show that the different gradient magnitude settings apparently affect the modeled fields, moreover, the AdaHRBF method is effective to iteratively obtain the true gradient magnitude of the fields. We modeled an analytic field of 1 ( ) = (( − 300) 2 + ( ) 2 ) 3 2 with the changing gradient direction and magnitude as show in Fig. 5a.
Then we sampled attribute and attitude points from the analytic field with different locations as shown in Fig. 5b. Hence, we can retrieve the coefficients and of the HRBF formula and the polynomial coefficients, respectively. We compared two different experimental settings: (1) Assuming that gradient is a unit vector and each gradient magnitude is 1, we used the HRBF interpolant to reconstruct the field as shown in Fig. 5c. Although the field values at the sampling points are equal to the given attribute values, the retrieved field values change irregularly, thus we obtained a large number of exceptional values in the reconstructed field.
(2) The true gradient magnitude was obtained via the iterative AdaHRBF method introduced above. In this condition, we more accurately restored the field (as shown in Fig. 5d) and also got the optimized gradient magnitude after 265 the iterations, which was close to the true value.

270
We also modeled a potential field of 2 ( ) = ( ) 3 with the changing gradient magnitude as show in Fig. 6a. It is known that each direction of gradient points is the positive Y-axis direction. We sampled attribute points and attitude points as shown in Fig. 6b. We also compared two different experimental conditions: (1) Assuming that each fixed gradient magnitude is 1, we used the HRBF interpolant to reconstruct the field as shown in Fig. 6c. (2)

280
We overlaid above-mentioned two fields to generate a new potential field of 3 ( ) = 1 ( ) + 2 ( ) as show in Fig. 7a, and the sampled attribute points and attitude points are shown in Fig. 7b. We also compared two different experimental conditions: (1) Assuming that each fixed gradient magnitude is 1, we used the HRBF interpolant to reconstruct the field as shown in Fig.   7c. (2) The true gradient magnitude was obtained via the iterative AdaHRBF method, and we more accurately restored the potential field (as shown in Fig. 7d).

290
The study area is located in the Lingnian-Ningping manganese ore zone, in Debao County, southwestern Guangxi Zhuang Autonomous Region, China (Fig. 8). The study area mainly consists of strata from the late Paleozoic to the late Triassic-Pliocene (T3-N2). The middle Permian (P2) strata are in para-unconformity contact with early Triassic (T1) strata; the middle Triassic (T2) strata are in angular unconformity contact with Quaternary. There is a left strike-slip inverse fault, the Nacha Fault, in the middle of the study area. It dips to the southeast, with a NE strike direction of 45°, a dip angle of about 70°, and 295 a total length of about 12 km, extending outside the study area. The footwall slid to the west relative to the hanging wall, and the slip distance is about 600 m. There are two synclines (I and Ⅲ ) and an anticline (II) in the study area. Syncline III is located in the middle of the study area with a high symmetry. The axis of syncline III strikes nearly northeast and its south limb is cut by the Nacha Fault. Anticline II is located in the northwest of the study area with a good symmetry, the fold axis striking about 30° northeast.  Faults, unconformable strata, and intrusive rocks all cause discontinuities in a SPF (Calcagno et al., 2008). We used the fault surface samplings to interpolate the potential field of the Nacha Fault (Fig. 9a). We extracted the zero equipotential surface of 305 the fault potential field to reconstruct the surface model of the Nacha Fault which divides the study area into two sub-domains ( Fig. 9b). In each sub-domain, the coefficients of the HRBF linear system were separately solved according to the joint samplings of the SPF and its gradient. According to the comprehensive stratigraphic column, the burial depth of each stratigraphic interface relative to the top surface of the Quaternary was used as the attribute value of the SPF (Fig. 10) for implicit function interpolation. The SPF defines the 3D space as a scalar function f(p) at any point p, where f is defined as the relative burial depth in this study. Burial depth decreases as geological time progresses; therefore, earlier deposited strata are assigned a relatively larger burial depth, while 315 later deposited strata are assigned a relatively smaller burial depth. Each point inside of a stratum has its own burial depth relative to the top surface of the Quaternary, therefore, the depth values in the field decrease gradually from bottom to top in strata. In this context, the SPF is fitted by a scalar function of the relative burial depth. When the relative burial depth is used as the attribute value of the SPF, we can set the initial gradient magnitude ‖ ‖ ≅ 1 if the strata underwent heterogenous deformation. However, if we use geological age as the attribute value of the SPF, ‖ ‖ can no longer be initially assumed to 320 be 1 because the stratigraphic age and distance along the gradient direction are from different measured variables.

330
The attribute points and attitude points of each stratigraphic interface and fault plane extracted from the geological map and cross-sections were used as the original dataset for 3D SPF modeling. The 3D points of stratigraphic interfaces extracted from the geological map and cross-sections were regarded as samplings of the SPF. The gradient vectors which are transformed from the off-contact stratigraphic attitude points were regarded as the samplings of the gradient of SPF.

335
There are 1410 known on-contact attribute points and 34 off-contact attitude points scattered throughout the study area (Fig.   12a). The known attitude sampling points are scattered on the south limb of fold I, the north and south limbs of fold II, and the north and south limbs of fold III. There are 17 attitude sampling points in the north side of the Nacha fault and 17 attitude sampling points on the south side. The distribution of the dip directions and dip angles is shown in Fig.12b. On a specific grid resolution, we modeled the scalar field of gradient magnitude before and after optimization for each attitude point (Fig. 14). Along the north side of the Nacha Fault in Fig. 14a, the gradient magnitudes obtained by interpolation in area 355 B exceed the maximum values. Compared with the scalar field of gradient magnitude before optimization, the scalar field of gradient magnitude after optimization (Fig. 14b) more smoothly represents changes in the strata. The Carboniferous strata have the largest true gradient magnitude, while he true gradient magnitudes of the Devonian strata are smallest. Furthermore, we cut four cross-sections of the gradient magnitude scalar field, as shown in Fig. 15.

Stratigraphic Potential Field (SPF)
After the optimized gradient magnitude for each attitude point was obtained, all scatted attribute points and attitude points were finally substituted into HRBF linear system to respectively solve the HRBF coefficients (α , ) and the polynomial coefficients ( 1 , 2 , 3 , 4 ) for each side of the Nacha Fault. On a specific grid resolution, we generated the regular discrete grids as interpolated points in 3D space. Then the points above the digital elevation model (DEM) were removed from the 370 interpolated points. Finally, we reconstruct the SPFs in 3D space before and after optimization of the gradient magnitude according to the respective HRBF interpolant of each sub-domain (Fig. 15). In this study, the SPF represents the relative burial depth in 3D space. The larger field value represents earlier deposited strata with larger relative burial depth, and vice versa.
The same stratigraphic interfaces in different sub-domains share the same field value. The field values change abruptly at the Nacha Fault because the conformable strata were cut by the fault plane.

375
The SPFs are both constrained so that the interpolated SPFs values at the attribute points are equal to the initial relative burial depths, but the SPFs values may abruptly change or produce outliers at some locations. Obviously, the SPF values change nonuniformly with gradient magnitude before optimization, which caused the SPF values that originally belonged to the Carboniferous strata to be interpolated as those of other strata and sequentially resulted in incorrectly extraction of the stratigraphic interfaces. The abnormal SPF values (areas A, B and C in Fig. 16a), are not continuously distributed along 380 stratigraphic interfaces but appear at irregular intervals. This abnormal stratigraphic potential field causes separated, discontinuous, and dispersed stratigraphic interfaces to be extracted through equipotential surface tracking. However, reconstructing the SPF through optimization of gradient magnitude for each attitude point (Fig. 16b) avoids the generation of either abnormal field values or of the wrong equipotential surfaces. This geologically plausible SPF can be appropriately constrained by the known gradient direction and the optimized gradient magnitude at the attitude sampling points. We cut the SPF along four section lines, and the SPF value also changes more uniformly from older to younger strata after gradient magnitude optimization than using a fixed gradient magnitude of 1, as shown in Fig. 17.

Three-Dimensional Models of Strata
Once the field was interpolated in 3D space, the specific equipotential surfaces were extracted from the implicit volumetric function as stratigraphic interfaces within each main structure bounded sub-domain. We used the marching cube method to 395 extract the equipotential surfaces with a specific relative burial depth from the stratigraphic interfaces by connecting all the points with the same field value in the stratigraphic potential field (Fig. 18). The 3D surface model extracted from the potential field shows that the geometrical shape of each equipotential (iso-depth) surface is smooth, and the topology is consistent. The interface model on both sides of the Nacha Fault restores the location of the fault in the south limb of syncline III. Sequentially, according to the range of relative burial depth of stratigraphic top and bottom, two stratigraphic solid models were reconstructed from these equipotential surfaces before and after optimization of gradient magnitude for each attitude point, respectively, combined with sub-domain boundaries and DEM (Fig. 19). Many abnormal potential field values and additional unreasonable geological bodies were extracted from the model before optimization, especially in areas A and C as shown in Fig. 19a. These abnormal potential field values lead to the occurrences of additional strata fragments that do not conform to the rule of sediments. Therefore, the HRBF interpolation with the initial fixed gradient magnitude of 1 roughly reflects stratigraphic on-contact information and captures the structure of syncline I in the north. However, several details are different from the stratigraphic structure on the geological map. Where the Nacha Fault passes through syncline III, the strata 410 on the south side of the fault plane should correspond to the same strata on the north side. However, the Devonian strata corresponded to the Permian strata in area B as shown in Fig. 18a, which is inconsistent with the geological structure. The geological model extracted using the optimized gradient magnitude for each attitude point is shown in Fig. 19b. Overall, the obtained geometries follow more closely the shape of the folds and stratigraphic on-contact lines. From north to south in the study area, anticline II and syncline III were successfully modeled with the Nacha Fault correctly represented as an inverse Four cross-sections through the solid models (see the geological map for cross-section lines) were cut, and the cross-sections of the solid model are more consistent with the original structural relationships on the geological map after gradient magnitude optimization than using a fixed gradient magnitude of 1, as shown in Fig. 20. The highest stratum and section coincidence percentages on cross-sections are 74.50% (T2) and 78.03% (Section 16) before optimization, respectively, as shown in Table 2. However, the highest stratum and section coincidence percentages on crosssections are 98.99% (D1) and 98.01% (Sections 13 and 15) using the optimized gradient magnitude for each attitude point, respectively, as shown in Table 3. The total coincidence percentage on cross-sections increases from 67.03% to 98.27% after optimizing gradient magnitude. In order to represent the variability of stratigraphic thickness, we introduce a definition of stratigraphic thickness index (STI) at any location p in a stratum: where ftop(p) and fbottom(p) are the potential field values of the top and bottom surfaces, respectively, of a stratum where the point p is located, and l(p) is the gradient magnitude obtained by the previously described iterations at the location p. The

440
normalized STI represents the true thickness of the stratum passing through the location p; therefore, we can analyze stratigraphic thinning and thickening effects by comparing STI at different locations in the stratum. For each stratum in the study area, the STI of each stratum may be different everywhere.
The STI, which is the ratio of the difference between the potential field values of a stratum's top and bottom surfaces to the gradient magnitude at each point, is a normalized indicator to represent the thickness variation of strata. We restored the STI 445 scalar field using the optimized gradient magnitude and found that the STI values depend on the strata in the study area (     We cut the STI scalar field along four section lines, as shown in Fig. 23. The STI at the core of syncline III gradually decreases from the northeast to the southwest, which reflects the southwest plunge of this syncline. At the core of anticline II, the STI is lower in the northeast, but is higher in the southwest where the thicker Devonian strata occur, which reflects the northeast plunge of this anticline. 465 Figure 23. Cross-sections of the stratigraphic thickness index (STI) field.

Discussions
Due to their robustness and stability, RBF and HRBF interpolants are a good choice for modeling when geological sampling data is relatively sparse and uneven (Cowan et al., 2003;Guo et al., 2016). However, existing RBF and HRBF interpolants 470 implicitly reconstruct a single geological interface and extract it as the zero-value equipotential surface. Moreover, existing RBF and HRBF interpolants need several independent scalar fields to simulate geological interfaces, which makes it difficult to ensure topological consistency among different geological interfaces. The AdaHRBF proposed in this study can create only one scalar field and extract multiple continuous stratigraphic interfaces for a set of conformable strata.
The AdaHRBF proposed in this study improves the use of altitude data in SPF modeling by optimization of gradient magnitudes.

475
In additional to use of attitude information as the gradient directions of SPFs, we use the gradient magnitude as a new constraint to control the rate of change of SPF values. The gradient of a SPF is a vector with certain direction and magnitude, in which https://doi.org/10.5194/egusphere-2022-1304 Preprint. Discussion started: 23 January 2023 c Author(s) 2023. CC BY 4.0 License.
the gradient magnitude provides constraints on the thickness of deformed strata. Therefore, it is extremely important to construct HRBF linear systems with accurate gradient magnitudes in 3D SPF modeling. As a "chicken-and-egg" problem, it is difficult to determine the exact gradient magnitude through the geological measurements or prior structural knowledge. We

480
proposed an iterative optimization method which alternates between estimation of SPF and gradient magnitudes so that the gradient magnitudes progressively converge towards the values being adaptive to the stratigraphic architecture. The optimized gradient magnitudes more accurately simulate the variations of the SPF between the top and bottom surfaces. Jessell et al. (2014) highlighted two limitations of current implicit modeling schemes: (1) they are incapable of interpolating or extrapolating a fold series within a continuous structural style; (2) the shape of fold hinges they produce is not controlled 485 and may yield inconsistent geometries. To overcome these two limitations, we adopted two strategies: (1) a 3D stratigraphic potential field modeling method based on HRBF interpolant was used to interpolate or extrapolate a fold series within a structurally continuous domain; (2) a number of structural attitude points were sampled on both limbs of the folds to control the geometries of fold hinges.
There are several choices for the value of the potential field, e.g., the sorted serial number or depositional time for each 490 stratigraphic interface (Mallet, 2004). However, the thickness of the stratum is not necessarily proportional to the sorted serial number and deposition time. We chose the burial depth of each stratigraphic interface relative to the top surface of the Quaternary as the SPF value so that stratigraphic thickness is considered as a constraint in the modeling. Compared with using the sorted serial number or depositional age of stratigraphic interfaces as the potential field value, our solution is more in line with 3D SPF modeling. We derived the gradient direction from the attitude points; moreover, we used the gradient magnitude 495 as a constraint to control the rate of change of the SPF.

Conclusions
The purpose of this study is to establish a framework for 3D SPF modeling by using the HRBF interpolant with adaptive gradient optimization constrained by on-contact attribute points and off-contact structural attitude points. We applied this method to a study site in the Lingnian-Ningping area, and a geological map, 4 cross-sections, and a DEM were used as original 500 data to model a SPF whose field value was taken from the relative burial depth of the stratigraphic interfaces. The results show that the implicit modeling of the SPF by HRBF interpolant and optimization of gradient magnitude can be effectively adapted to 3D geological modeling using the sampling points from a geological map and cross-sections. A SPF can express the parameters of a stratum such as property, shape and topology in 3D space. Because 3D stratigraphic potential fields can be coupled with various geoscience numerical simulation methods, they have a broad prospect for application in related fields 505 such as metallogenic prediction.
However, the modeling process is complicated because the sub-domains are required to be divided manually. In actual geological surveys, the geological structure may be more complex and include a large number of faults, unconformable strata https://doi.org/10.5194/egusphere-2022-1304 Preprint. Discussion started: 23 January 2023 c Author(s) 2023. CC BY 4.0 License.
(Central South University) for their kind assistance with data collection and Professor Jeffrey Dick (Central South University) for revising scientific English writing of this manuscript.