A Triangulated Irregular Network Constrained Ordinary Kriging Method for Three-Dimensional Modeling of Faulted Geological Surfaces

Three-dimensional (3D) modeling of geological surfaces, such as coal seams and strata horizons, from sparsely sampled data collected in the field, is a crucial task in geological modeling. Interpolation is a common approach for this task to construct continuous geological surface models. However, this problem becomes challenging considering the impact of the faults on geological surfaces. Existing methods tend to solve this problem through three steps, including interpolating stratum and fault surface, applying a fault modeling method to modify the geological surface, and optimizing the modified surface to pass sample points fallen into the fault displacement zone. This paper presents a more concise method to generate a faulted geological surface, in which 1) a constrained Delaunay triangulated irregular network (CD-TIN) is constructed to facilitate the neighborhood search process of the ordinary kriging (OK) interpolation, 2) the CD-TIN is also directly constrained by horizon cut-off lines formed from theoretical fault displacement profiles, and 3) subsequently, neighbors of the location to be estimated are selected effectively in the CD-TIN considering the fault topology. The proposed method significantly improves the time efficiency of the OK interpolation by utilizing the CD-TIN and incorporates fault effects directly into the interpolation process by inserting fault horizontal cut-off lines into CD-TIN. Moreover, by integrating the fault effects directly into the interpolation process, the surface modeling process is accomplished in a single stage instead of two separate stages of interpolation first and then modifying the surface in the fault area. By this strategy, the proposed method significantly improves the time efficiency of the OK interpolation algorithm and achieves more accurate modeling of the faulted geological surface. Experiments were designed to compare the performance of our method with several commonly used approaches, and the results indicate that the proposed TIN-constrained OK method achieves better accuracy and efficiency in modeling faulted geological surfaces than other methods. This method could also be used in geospatial interpolation studies, such as meteorological data interpolation.


I. INTRODUCTION
In subsurface studies, three-dimensional (3D) structural modeling helps geologists understand the physical world [1]- [3] and has been widely utilized in applications such as resource The associate editor coordinating the review of this manuscript and approving it for publication was X. Huang . management, disaster control, and tectonic evolution. Faults lead to the distributed damage of geological layers, and constructing faults in structural models play a vital role in mitigating risks in mining design. However, many works concentrate on the modeling of the physical properties of faults [4], [5], whereas the distributed damage of the surrounding geological layers draws less attention. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Constructing a structural model with layers and faults typically involves two steps. First, the surfaces of geological layers, i.e. the surfaces that represent the boundary between two rock layers, is modeled from the collected data. Second, the surface is cut and displaced by faults using a fault modeling method.
Traditionally, the modeling of surfaces relies on defining geological boundaries and linking cross sections to generate 3D surfaces manually [6], [7]. Automatic modeling methods based on extensive interpolators are then applied. These methods interpolate sample points extracted from horizontal observations of the geological surface to reconstruct a surface model, and they can be categorized into explicit and implicit approaches. Explicit methods, such as the inverse distance weighting interpolation method, and geostatistical methods (mainly kriging techniques; [8]- [11]) are utilized to explicitly calculate a continuous 3D surface by a function z = f (x) that interpolates properties Z = {z i , i = 0, 1, . . . , n} of 2D sample points P = {x i ( x i , y i ), i = 0, 1, . . . , n} to get 3D points (x, y, z) on the surface. Alternatively, in the implicit surface reconstruction method, a geological surface is considered as an isovalue surface of a 3D scalar field f (x, y, z) interpolated from the collected data. Commonly used implicit methods include discrete smooth interpolation [12], [13], and a set of radial basis functions (RBF) [14], [15]. RBF techniques, which have an extensive application in computer-aided design, interpolate a scalar field f = f (x, y, z) based on the observations and then extract and visualize the zero isovalue surface (f = 0) of the scalar field. The Hermite RBF (HRBF) method [16] leveraging the data to be interpolated with second-order information (normal vectors). A compactly supported RBF method is given to decrease the interpolation time by constraining the neighborhood radius. Reference [17] and [18] developed a multiscale CSRBF approach using an Octree-based spatial index and normal vectors, advances in both time efficiency and interpolation accuracy.
Faults are then introduced to intersect and adjust the surface. Traditional fault modeling methods insert topological discontinuities by cutting and moving the smooth surface in a 3D Cartesian coordinate system [18], [19]. Then, implicit approaches are introduced to interpolate the fault surface into an implicit surface of a scalar field, then the regional model is modified in this scalar field [20]- [22]. In most of these methods, a fault is presented as topological discontinuity in a geological layer. First, fault geometries are computed by interpolating fault observations or estimated from seismic images based on semblance (e.g., [24], [25]). Second, fault surfaces are intersected with geological surfaces. Third, to keep the kinematic consistency of the surface, local movements of the intersection lines, namely, the fault cutoff lines, needs to be performed and propagated to the whole surface by manual [2], [26]. This fault modeling process increases the time required to build a valid model [27], so mathematical fault displacement modeling methods are suggested [28], [29], in which the slips on and around the fault geometry are modeled by interpolation [30] or 3D parametric fault representation techniques [20], [31].
Although these methods are widely used in geological applications, an interpolated surface, which is initially constrained by all sample points and then cut and moved by faults, may no longer pass some of these controlled points. Fitting observations to get a geological surface that honors both horizontal and fault observations become difficult [32]. Optimization methods, such as numerical optimization and particle swarm optimization, can be used to move the faulted surface to honor the constraint of observations. The former finds an ''optimal'' set of parameters for trishear algorithms of fault modeling [33], [34], and the latter chooses kinematic parameters by minimizing a cost function consisting of the misfit between the data observed and the built surface model [32]. However, the complexity of the modeling workflow grows, as numerical calculation and kinematic interpretation is introduced by these methods, and, therefore, they should be implemented cautiously.
In contrast to the two-phase modeling process of these previous methods, we proposed a new workflow to reconstruct a geological surface with the presence of faults in a single processing step. As previous researches, this method also presents fault by topological discontinuity. However, it combines the fault modeling and interpolation processes and produces a geological surface honoring the constraint of faults and original sampling points. Geological modeling tasks are set to benefit from this method at two points: 1) no calculation required for the displacement of geological layer, and 2) the optimization methods do not need to be introduced, which will significantly improve the time efficiency of geological modeling. This method is achieved by leveraging a triangulated irregular network (TIN)-constrained ordinary kriging (OK) method. The OK method is used as the main interpolation method because of its effectiveness in geological applications and the simplicity of its mathematical form. The constrained Delaunay triangulations (CDT) algorithm [35] is performed to construct a TIN using the sample points and fault cutoff lines to constrain its neighborhood search process for locating neighbor points used in the interpolation. A workflow is then developed to construct the 3D models of faulted geological surfaces. The experimental results show that the proposed interpolation method improves both the accuracy and efficiency of the faulted geological surface model compared to the OK, HRBF, and multiscale CSRBF methods.
The remainder of this paper proceeds as follows. Section 2 highlights the related works in the literature. In Section 3, the framework of the 3D modeling of the faulted surface is described in detail. Section 4 describes the experimental results. Section 5 concludes the paper and discusses further research directions.

II. RELATED WORK
Kriging techniques provide the best linear unbiased estimation for un-sampled points as well as the ability to quantify the uncertainty of estimates by error variance [36], [37], and they are frequently adopted to build geological layers [10].
To implement kriging, the neighborhood search method is needed to find nearby points of location for estimation. The distance-based neighborhood method, namely, identifying sample points based on their distance to the target location, is commonly implemented to accelerate the interpolation process [38], [39]. To improve the spatial balance of selected data, references [39] and [8] divided the space into quadrants around the target location and selected data within each (sub)quadrant based on distance. Other criterion for defining a neighborhood include K-nearest neighbors (KNN), which selects k closest points to a target point using a spatial data structure constructed from the whole sample dataset. The introduced spatial data structures in kriging studies include the k-dimensional tree [40]- [42] and the fast KNN search [43]. These structures provide storage and a fast index for the collected sample points and can avoid brute force computation of distances between the target point and points in the whole dataset.
However, these neighborhood search methods may not be optimal insofar as (i) they ignore the spatial topology among data, and only a geometric standpoint is applied to measure the closeness and redundancy of data [44]. The spatial topology, which defines containment, connectivity, and adjacency for spatial features, can be helpful in supporting better data selection. For example, a sample point (say point A in Fig. 1) that is geographically close to the target point for which the value is to be interpolated (star point in Fig. 1) but is located on the opposite side of the fault surface should be not be considered as a neighbor and should be excluded from the subsequent interpolation calculation. (ii) The brute force computation of distances between interpolated points and all sample points in the collected dataset is a very timeconsuming procedure, particularly when the interpolation is conducted using big data [40].
Since a fault in a geological layer can be represented as a linear feature (for an unbounded fault terminating at the boundary of the modeling domain) or a polygonal feature (for a finite fault with a limited length in the modeling domain), another strategy, based on line intersection, is adopted according to [19]. For each sample point, a line segment is drawn to connect it with the location for which its value is to be estimated, and it is verified whether this line segment intersects with the fault surface. Those verified sample points are then filtered from the result of the neighborhood search in each estimation. However, efficiency is a major concern in this approach, as the intersection operation in each estimation can be time consuming because the size of sample dataset, fault numbers, and estimation locations increases.
To address the aforementioned problems in the neighborhood search of the kriging method, our proposed TINconstrained neighborhood search method uses a TIN to restrain the neighborhood search process. TIN comprises a series of vertices, which are the locations of sampled observations of the geological layer here, with 3D associated coordinates connected by edges to form a triangular tessellation. The CDT algorithm used in triangulation can force certain required segments (linear features) into the triangulation. Thus, by inserting a polygon consisting of fault cutoff lines into a TIN using the CDT algorithm and removing the internal triangles of this polygon, a fault can be represented as a hole bounded by fault cutoff lines in the TIN.
With the constraints of neighborhood defined in the TIN, the fault effect can be automatically incorporated during the interpolation process. Moreover, the spatial structure of the TIN stores adjacent triangles for every triangle in the TIN. As all vertices of triangles are sample points, the neighbors in triangles adjacent to the triangle containing the target point to be estimated can be reached in O(1) time, which significantly reduces the amount of computation needed for neighborhood search and accelerates the entire interpolation process.

A. METHOD OVERVIEW
There are four steps involved in the proposed 3D modeling of faulted geological surface workflow (Fig. 2). These are (i) preprocessing the collected data, i.e. the extraction of sample points, fault data, and boundary, (ii) computing VOLUME 8, 2020 fault geometries to obtain horizontal cutoff lines for each fault, (iii) interpolating sampled points with the constraint of fault horizontal cutoff lines, and (iv) building the geological surface within the geological boundary. Each step is listed below and described in detail later in this paper:

1) PREPROCESS THE COLLECTED DATA
Sample points of the geological surface can be extracted from archival geological data, such as drilling data, and geological maps, such as cross sections. Fault observations, including locations and displacement parameters, and boundary of the geological layer are also extracted from geological maps in this step.

2) GET HORIZONTAL CUTOFF LINES FOR EACH FAULT
The fault horizontal cutoff lines, defined as the hanging wall line and the footwall line for a given fault and its intersecting layer, are generated in this step by constructing a fault local coordinate system and calculating displacement along the fault surface.

3) IMPLEMENT INTERPOLATION
The TIN-constrained OK method is implemented in this step, as detailed in Section 3.3, with an input of the extracted geological sample points and computed fault horizontal cutoff lines and an output of a set of points located on the geological surface.

4) BUILD THE SURFACE WITHIN THE GEOLOGICAL BOUNDARY
After combining boundary and fault horizontal cutoff lines with estimated points within the boundary, faulted geological surface are then constructed by using the CDT algorithm based on these data.

B. GET FAULT HORIZONTAL CUTOFF LINES
In this work, we adopt three steps to calculate the horizontal cutoff lines for each fault: (i) compute a fault surface, (ii) build a local coordinate system for each fault, and (iii) calculate the horizontal cutoff lines.
In addition to the observed points on the fault surface, the parameters listed in Table 1 are also used to establish a fault surface and a local coordinate system, including the   [32]. By creating a fault medium plane (Fig. 3a), which is a reference plane aligned to v s and v d , signed distances can be defined as the distances between points and the plane (Fig. 3b). The geometry of a fault can then be computed (Fig. 3c) by interpolating these signed distances, and the fault local coordinate system G F can be established. The coordinates of a point X = {x, y, z} within To acquire the fault cutoff lines, profiles indicating the relationship between the fault extent and amplitude are expected, which can be utilized to calculate the displacement for each point on the fault cutoff lines with fault characteristics. Basically, the amplitude of the fault displacement should be at its maximum at point O and decreases to zero over a distance. Some other profiles are also available [19]. A parametric throw profile [20], using Hermite splines to present the attenuation followed by along-dip D s and along-strike D d (Fig.  4a), also can be used. Other profiles could be used instead, and more complex profiles could be obtained by mixing several simple profiles [45].
Fault horizontal cutoff lines can then be calculated along the fault surface in the local system and transformed to the coordinate system (usually Cartesian) where the sample points are represented (Fig. 5a). Fault cutoff lines and horizontal observation points are then interpolated to reconstruct the geological surface. Fig. 5b illustrates these data and a possible geological surface implementation.
where λ i is the weight parameters corresponding to each sample point determined by variogram, which defines the variance of the difference between values at two locations, such as an exponential model in which h is the distance between two sample points, a is the range, C 0 denote the nugget effect and C denote the partial sill. The last three parameters are fitted from an empirical variogram: where γ ij denotes the semivariance value between x i and x j , and k is the number of sample points used in the estimation.
After weights are calculated,ŷ (x 0 ) can be calculated from Eq. (1). However, as the volume of sample in the dataset grows, the dimensions of the matrices constructed by the dataset in Eq. (2) increases, which leads to a rapid increase in the time cost. The aforementioned distance-based or quadrantdistance-based neighborhood search methods (hereinafter, distance-based methods) are employed to limit the count k of sample points that participated in the estimation by selecting nearby neighbors.
The neighborhood search process follows Tobler's first law of geography, which states, ''Everything is related to everything else, but near things are more related than distant things'' [46]. When the neighbor count is limited, the time required to solve Eq. (2) can be theoretically constant in each estimation process, which helps when applying OK interpolation to a large dataset.

2) TIN-BASED NEIGHBORHOOD SEARCH METHOD
We explored the TIN-based neighborhood search method to tackle two challenges-the utilizing of spatial topology among data and the time-consuming nature of distance computations mentioned in Section 2-by constraining and accelerating the neighborhood search process of the kriging method. In our method, the neighborhood search process would benefit from the structure of the TIN, not only to applying fault's effect on neighborhood selection but also to accelerate this process by accessing adjacent triangles in the TIN.

a: CONSTRUCT TIN-BASED NEIGHBORHOOD SEARCH NETWORK
A TIN-based neighborhood search network, constructed by a CDT algorithm using sample points and fault cutoff lines, is defined to store the spatial relationship and topology among these data. A Poly2Tri library (http://code.google.com/p/poly2tri/) is used to the implement CDT algorithm, and the steps to establish the search network are as follows: -First, project sample points and fault cutoff lines of the surface to a plane. This is because the Poly2Tri library was developed for a 2D data set.
-Second, construct a triangulation network within the geological boundary by applying the CDT algorithm to these data.
The hanging wall and footwall of a reverse fault overlaps when the data is projected to 2D space (Fig. 6a). However, the CDT algorithm cannot generate overlapping areas-for example, Fig. 6b illustrates a TIN constructed by the CDT algorithm using the sample points and cutoff lines of the fault shown in Fig. 6a. Here, a non-overlapping discontinuous area, bounded by fault cutoff lines, is created, and this appears to be not in line with Fig. 6a. So, to model a reverse fault in the TIN and keep a correct topology for sample points, lines linked to a point in the hanging wall line (brown) and a corresponding point in the footwall line (green) are reconnected (Fig. 6c).

b: NEIGHBORHOOD SEARCH
Since each vertex in the TIN-based search network denotes a sample point or a node on the fault cutoff lines, the neighborhood search process can be executed by utilizing the structure of TIN. Fig. 7 shows the basic geometry components and their relationship in a constrained Delaunay TIN. All elements are categorized into three classes-Vertex, Edge, and Triangle. An Edge object serves as one of the three edges comprised in a Triangle object and connects two sample points storing their locations in Vertex objects. Each Triangle object also keeps its adjacent triangles, which can be accessed by the member function of GetNeighbors().
Then, to choose neighbor points in the search network to estimate the value at target location x, two steps are expected. The first step is to locate the triangle containing x, and this is achieved by a two-phase procedure of selecting triangles whose enclosing rectangle contains x, and locating the triangle T 0 containing x from these triangles. The order of the triangle, which measures the connectivity from x to each triangle in the TIN, is defined in the following step. The 0 th order triangle denotes T 0 confirmed in the first step, and the t th order triangles (t > 0) are given the definition of triangles that shares vertices or edges with (t-1) th order triangles. Then, by specifying the order of triangle t, all the vertices in 0∼t order triangles will be selected and used in the estimation.
To define the difference between the TIN-based neighborhood method and the distance-based neighborhood method, we should consider the following scenario. Suppose we select neighbor points from dataset X = {X i (x i , Z i (x i )), i = 1, 2, . . . , n} to estimate the value Z(x 0 ) at point X 0 . Using the distance-based neighborhood search, all sample points are sorted according to their distance to X 0 , and the six nearest points, i.e. X 4 ∼X 9 , were chosen (Fig. 8a). Even when the local surface is divided into two sides by fault f , x 8 and x 9 are apparently located at the opposite side of X 0 of the fault.
As a comparison, in our TIN-based neighborhood search method, all the vertices of 0∼t order triangles were selected in this process and subsequently interpolated. In Fig. 8, suppose t = 1, vertices in 0 th order (sample points X 4 , X 6 , and X 7 ) and 1 st order triangles (sample points X 1 , X 2 , X 3 , and X 5 and fault cutoff line points X f 1 ∼X f 4 ) were selected (Fig. 8b). All of these neighborhood points are located in a continuously changing region where X 0 is set and are obviously more meaningful for estimating the value at X 0 .
Compared with the other KNN methods mentioned in Section 2, only the TIN-based neighborhood search concerns the constraints of faults. By inserting faults in the TIN and using the TIN as a search network, this method can select sample points that are not only spatially close but also highly topologically correlated to the target point of estimation, which ensures that the prediction of the OK method is produced from a continuous domain.

A. MODELING OF SIMULATED DATASET
To evaluate the performance of the proposed TIN-constrained OK method, a surface model of a geological surface S was simulated from a process z = s (x, y), where x, y [−50, 50] m and z = 0 m. A finite normal fault f was then simulated to intersect, cut, and deform S, with parameters listed in Table 2.
A data set P = {p i (x i , y i , z i ), i = 1, 2, . . . , 100} including 100 points was randomly sampled from the faulted surface, and then normal vectors of the surface corresponding to each sample point were calculated. The simulated surface, fault, and sampled data set are shown in Fig. 8a.
In our experiments, OK interpolation methods using two distance-based neighborhood search methods were employed to compare the performance in both accuracy and efficiency with TIN-constrained OK. Two implicit methods using second-order information, including HRBF, and multiscale CSRBF methods, were also added to the comparison. By leveraging the gradient constraints, the latter two methods provide better interpolation results than the original RBF method.
For readability, the used methods are abbreviated as TIN-OK, D-OK (OK method using distance-based neighborhood), Q-OK (OK method using quadrant-distance-based neighborhood), HRBF, and multi-CSRBF (multiscale CSRBF).

1) INTERPOLATION RESULTS OF SIMULATED DATASET
In the first scenario, four interpolators were applied to the simulated dataset to estimate z values for unsampled locations to form surface models drawn in Fig. 9 b, c, d, and e, respectively. All surfaces are rendered within the same color range from light purple to red (corresponding to the value range [−5 √ 3, 5 √ 3] m of z), and the dark areas shown in the figures indicate that the estimated values exceed this value range.
From Fig. 9, the values predicted from the D-OK and Q-OK methods in the fault displacement zone are narrowed because they chose sample points from both sides of fault when estimating values of locations near fault. The surfaces modeled from HRBF and multi-CSRBF change sharply near the fault when additional normal data are used (Fig. 9e). In line with expectations, our proposed TIN-OK performs better than the other methods with a surface model (Fig. 9b) that very similar to the simulated ''ground truth'' surface (Fig. 9a).
A further experiment was conducted to compare the performance of the TIN-OK and D-OK method on fault modeling. For our proposed TIN-OK method, the fault can be constructed with the constraint of fault horizontal cutoff lines in TIN; for the D-OK method, to construct fault with OK interpolated surface, a fault modeling method using a fault  local frame was implemented based on the work presented in reference [32] and detailed in section III. B. Fig. 10 depicts the fault modeling result using an additional fault modeling method after the D-OK interpolation. Compare to the result illustrated in Fig. 9(b), the surface in Fig. 10 appears some dark areas within the fault displacement zone. Because the deformation of the surface in the fault modeling process caused a misfit of sample points. Although it can be corrected by implementing optimization methods mentioned in Section I, the increasing time consumption is another concern. The same conclusion can be derived from the Q-OK, HRBF, and multi-CSRBF methods, since they also have no built-in methods to handle fault structures.

2) COMPARATION OF INTERPOLATION ACCURACY
A K-fold cross-validation (with K equal to ten here) was subsequently implemented to compare the accuracy of these interpolators. In each cross-validation, sample points were divided into K subsets, of which nine were interpolated, and the remaining one was used to acquire errors to validate. This process ends after all subsets have served as a validation set. The interpolation error at point (x i , y i ) is defined as an unsigned value of difference between the estimated value z i and the z i . This K-fold cross-validation operated ten times for each method in our experiment. The average errors from ten cross-validations were then interpolated and are illustrated in Fig. 11. The darker the color on the map, the greater the error in cross-validation.
In all four maps, the error increases when the sample point approaches the center of fault. The results from D-OK and Q-OK interpolation have larger dark areas (bigger errors), and the errors reach a maximum value of 4.55 m for the D-OK method (Fig. 11b) and 4.53 m for the Q-OK method (Fig. 11c). However, errors in most areas in the results from the proposed TIN-OK interpolation were below 0.40, with a maximum value of 1.52 (Fig. 11a). HRBF has better accuracy than D-OK and Q-OK, but it exceeded others in its error area, the same as a maximum value of 9.17 m (Fig. 11d). The multi-CSRBF approach also produced large errors in the center and two of four corners (Fig. 11e), maximally 7.19 m, which may be caused by a sharp change in the normal vectors near the fault. Thus, the proposed TIN-OK showed the highest interpolation accuracy compared to other approaches.  One of ten K-fold cross-validation results is illustrated in the boxplot in Fig. 12. The result shows the same trend illustrated in Fig. 11, i.e. the introducing of TIN improves the accuracy of OK method when interpolating faulted surface.

3) COMPARISON OF INTERPOLATION TIME
The time efficiency of the four methods was compared using simulated datasets with different sizes (50, 100, 500, 1000, 5000, 10000, 20000, 30000, 40000, and 50000) on a laptop with an Intel R Core TM i5-8250U processor clocked at 1.60 GHz and 8 GB of random-access memory. Fig. 13a shows the time consumption trend of the four methods when the size of the simulated dataset increases. Limiting the count of neighbors, selected by the different neighborhood methods, that participate in each estimation is expected to linearly increase time in all kriging methods. The slope of the line produced by the TIN-OK method is smaller than the other two OK methods with a reduction in the time used of 96.1% (D-OK) and 91.3% (Q-OK) when the count of simulated points reaches fifty thousand. In terms of the other methods, the HRBF method uses all points in each estimation, which means that the computational cost required to solve the system of RBF equations increases cubically and the laptop used in this experiment cannot implement interpolation with more than one thousand points; the multi-CSRBF method, in the opposite situation, runs even more efficiently than all of the TIN-OK method benefits from the Octree indexing and multiscale interpolations (Fig. 13b). Table 3 summarizes the time complexity of the four comparison methods. The whole interpolation process consists of three parts in our time consumption analysis: constructing a sample point index for the neighborhood search, interpolating, and visualizing. In the first part, the TIN-OK and multi-CSRBF methods have a worst-case time complexity of O(nlogn) for operating the CDT and Octree construction algorithms, respectively, for the whole dataset whereas D-OK, Q-OK, and HRBF cost no time in this stage.
In the interpolating stage, for the three OK methods (TIN-OK, D-OK, and Q-OK), each estimation in interpolation is calculated after indexing k neighbors and making the prediction. TIN-OK has an O(n) time complexity to locate a triangle, as we mentioned in Section 3.3.2.2; the D-OK and Q-OK methods sort sample points based on distance (in each sector for the Q-OK method) and create an index dynamically to find the k nearest neighbors in each estimation, which produces a worst-case time complexity of O(n) and O(nlogn), respectively. All of these methods index neighbors with O(1) time complexity, except the HRBF. Then, to solve the linear equation constructed in the OK methods and get an estimate from selected neighbors, the O(k 3 ) time is taken; the HRBF reaches O((n + 3n) 3 ) complexity to calculate the isovalue surface in this process; whereas the multi-CSRBF, using a sparse matrix benefit from the compactly supported radial basis neighborhood search, can approach O(n 2 ) time complexity in this stage.
The HRBF and multi-CSRBF also requires a visualization algorithm to explicitly represent the isovalue surface, such as the marching cubes algorithm with a worst-case time complexity of O(q), where q depends on the spatial granularity of the visualized surface.
In summary, in terms of time efficiency, although the TIN-OK method costs extra time in constructing the TIN-based search network, this method accelerates the neighbor selection process and leads to an improvement in overall time efficiency. Furthermore, fault displacement and optimization methods are needed after geological surfaces are interpolated from the D-OK, Q-OK, HRBF, and multi-CSRBF methods in a geological surface modeling workflow, which makes the modeling workflow even more complicated and may cause unpredictable time expenditure. As a comparison, our proposed method has no such steps because fault displacements have been considered in the interpolation step.

B. CASE STUDY
The Qianjiaying coal mine is located approximately 15 kilometers southeast of Tangshan, China, with an area of 88 square kilometers (Fig. 14). The strata in this mine have a smaller dip angle (generally 10 • ∼15 • ) and simpler construction. The coal-bearing strata of this field belong to the Upper Carboniferous and the Lower Permian, and the basement strata are the limestone of the Middle Ordovician Majiagou Formation. The total thickness of the coal-bearing strata is about 500 meters, containing more than ten layers of coal with a total thickness of 19.79 meters and a coal content of 3.96%.
Five horizontal coal seams (C05, C07, C08, C09, and C12a) and 26 faults are expected to be modeled in this area. Available data sources are listed in Table 4, including boreholes drilling, underground sketch mapping, and fault observation. VOLUME 8, 2020   Based on the proposed workflow, fault surfaces were modeled (Fig. 15a), and their horizontal cutoff lines were computed for each intersecting fault and coal seam. Then, the sample points extracted from drills and underground surveys and fault horizontal cutoff lines were interpolated using the proposed method. Finally, all faulted coal seams models were reconstructed (Fig. 15b). Fig. 15c illustrated detailed models of two faults intersecting coal seams.

V. CONCLUSION
In this paper, a TIN-constrained OK method is proposed to build the faulted surfaces in geological context. To the best of our knowledge, this is a pioneering attempt to apply TIN to kriging as a neighborhood search constraint. Compared to previous methods, this method makes three important methodological contributions: (1) A search network based on TIN is developed, and neighborhood search is conducted along this network, which ensures that each estimation uses sample points from a stable and continuous area.
(2) The proposed method performs better in both accuracy and efficiency than the OK method using distance-based neighborhood search approaches and the HRBF method in the geological modeling of the faulted surface, particularly regarding near faults.
(3) The fault displacement and the optimization steps after interpolation are simplified in our method, which may also improve the time efficiency of modeling.
We provide a simplified modeling workflow based on the proposed TIN-constrained OK method without introducing the complex kinematics and mechanics of faults. The accuracy of modeling improves as the available data increases, and time efficiency exceeds two OK methods using the distancebased neighborhood method, and two implicit approaches (HRBF and multi-CSRBF).
However, because the fault displacement direction is not always parallel to the slip direction at the fault center, the method of computing cutoff lines does not always apply to all fault types [31]. In the future, other fault modeling methods will be introduced to address this problem [19]. Moreover, as an explicit modeling method, the proposed method present elevation value (z) in terms of x and y coordinates, which means it cannot be applied to present the overfolding or erosion of a surface such as an overturned fold or salt surfaces, etc. A possible solution is to project such structure to a medium plane aligned with the slope direction of such surface [31], or by using a central cylindrical projection.