Loose layers are the locus of human activities. The high-quality 3D modeling of loose layers has essential research significance and applicability in engineering geology, hydraulic and hydroelectric engineering, and urban underground space design. To address the shortcomings of traditional 3D loose-layer modeling based on borehole data, such as the lack of bedrock surface constraints, simple strata pinch-out processing, and the higher fitting error of the strata surface, a 3D loose-layer modeling method based on the stratum development law is proposed. The method mainly uses three different virtual boreholes, bedrock-boundary virtual boreholes, pinch-out virtual boreholes, and densified virtual boreholes, to control the stratigraphic distribution. Case studies demonstrate the effectiveness of this 3D loose-layer modeling method in the Qinhuai District of Nanjing and Hangkonggang District of Zhengzhou. Compared to the previous methods that interpolated stratigraphic surfaces with elevation information, the method proposed in this article interpolates the stratum thickness based on stacking, which could improve the interpolation accuracy. In the area where the loose layers and exposed bedrock are alternately distributed, stratigraphic thickness errors’ mean and standard deviation decreased by 2.11 and 2.13 m. In the pure loose-layer area, they dropped by 0.96 and 0.33 m. In addition, the proposed approach allows us to infer the different stratigraphic distribution patterns accurately and complete 3D loose-layer model construction with higher accuracy and a good visualization effect.
Three-dimensional models of geological bodies can truly and intuitively reflect the 3D shape and the spatial distribution law of geological bodies, accurately express the boundary conditions of complex geological bodies, and show the spatial distribution law of different geological structures. These models have essential research significance and application value in the fields of 3D geological analysis [1,2,3,4,5,6,7,8,9], underground engineering-aided design [9,10,11,12,13], and urban 3D evaluation [14,15,16,17,18]. According to the classification method of strata in China, the layers composed of related sediments or accumulations in Quaternary and Neogene, such as alluvium and diluvium, are defined as loose layers [19,20]. As one of the ordinary and widely distributed geological bodies on the Earth’s surface, the loose layer is the locus of human activities . Its resources and environmental characteristics affect and enhance human life, production, and social development. The rational use of this region is necessary to ensure human survival and development. Therefore, the high-precision and high-quality 3D modeling of loose layers has become a trending research topic because of its essential reference value in engineering geology [22,23], hydraulic and hydroelectric engineering , and urban underground space design .
Engineering drilling is the primary method used to obtain underground 3D spatial information and is an important data source for estimating the deformation modulus of a stratified rock mass [26,27], underground planning and analysis [8,28,29], and 3D geological modeling [2,12,30,31,32]. Therefore, the modeling method using borehole data has always been the mainstream 3D model construction technology for loose layers. Currently, common modeling methods based on borehole data mainly include the horizons-to-solids algorithm , the strata-framework algorithm [33,34,35], the boreholes-interfaces algorithm [36,37], the boreholes-voxel-solids algorithm [38,39], the virtual boreholes algorithm [35,37,40], and the machine learning algorithm [2,41,42]. The horizons-to-solids algorithm  is a method used for uniformly numbering the stratigraphic contact, interpolating each stratigraphic surface based on the control points of each stratigraphic contact to construct triangulated irregular network (TIN), cutting the stratigraphic surfaces, and finally generating a 3D model of the loose layer. Based on constructing a 3D geological entity using Boolean operations, this method further integrates the user-defined cross-sections to optimize the strata pinch-out position. However, in the fitting process of stratigraphic surfaces, when there are intersections or gaps between the interpolation stratigraphic surfaces, the elimination method of intersections or gaps by adjusting the elevation may make the modeling results inconsistent with the actual distribution of the strata. The strata-framework algorithm [33,34,35] and the horizons-to-solids algorithm are consistent in the stratigraphic numbering rules and the intersection processing method of interpolation stratigraphic surfaces. The only difference is that the stratigraphic surfaces of the latter are expressed in the grid. This method is effective for geological bodies with clear sedimentary sequences and horizontal or nearly horizontal distributions. However, the modeling effect is poor for strata with complex geological tectonics and serious deformation. The borehole-interface algorithm [36,37] is a 3D modeling method using measured boreholes and virtual boreholes. This method can introduce expert knowledge through the use of virtual boreholes and support the modeling of faults. However, the pinch-out positions of the missing strata in this method are all located on the borehole. When the strata are continuously missing, the missing strata are all pinched out to the same point, which does not conform to the law of gradual pinch-out of strata. The borehole-voxel-solid algorithm [38,39] is a 3D modeling method using a generalized tri-prism (GTP), which is generated by expanding down the three boreholes in a mesh. In this method, the stratum model takes GTP as the basic voxel, which can realize the 3D automatic modeling of the stratum with faults and can be easily transformed into a tetrahedral structure. However, this method cannot effectively solve the expression of the intersection position of stratigraphic surfaces and the treatment of abnormal fragments. The virtual borehole algorithm [35,37,40] is a method that adds multiple hypothetical virtual boreholes based on expert experience or other survey data. This method can generate a more refined and reasonable model but simply determines the pinch-out position according to the stratum thicknesses. Continuous pinch-out strata whose thickness values are in the same range pinch out to the same position on the same virtual borehole, which does not follow the law of gradual pinch-out. The machine learning algorithm [2,41,42] transformed the strata 3D modeling problem into a process of geological attribute classification of the underground spatial grid units. On the basis of the classification algorithms of support vector machine and deep neural network, they can use different spatial interpolate methods to fit geological surfaces to realize automatic 3D geological modeling from borehole data. This method can use limited and sparse borehole data to build a three-dimensional stratigraphic model and reduce manual intervention efficiently. However, they rely too heavily on the quality of samples, and the quality of modeling in different regions is prone to polarization. Therefore, the relevant modeling process needs to be further verified. In conclusion, the related modeling methods given earlier provide a positive exploration of loose layer modeling and basically meet the needs of scientific research and application.
However, when addressing actual engineering application requirements, the existing loose layer modeling methods based on borehole data have some gaps and deficiencies in bedrock surface boundary constraints, strata pinch-out treatment, and stratigraphic surface interpolation. (1) Currently, the restriction of bedrock on the loose layer is basically not considered in the modeling process. This means the strata spatial distribution law of gradual pinch-out at the boundary of the loose layer is not fully reflected (Figure 1a). In this case, there are usually intersections and gaps between the bedrock model and the loose layer model, and they cannot be seamlessly integrated. (2) In terms of stratum pinch-out, the processing results of many studies [37,40,43,44] do not conform to the spatial distribution law of missing strata that continuous missing strata pinch out gradually and discontinuous missing strata pinch out alternately (Figure 1b). (3) Currently, the top and bottom surfaces of each stratum are fitted by interpolating the measured borehole inner strata elevation. On the one hand, this treatment easily causes the intersection of the top and bottom surfaces (Figure 1c) due to high elevation error [32,35]. On the other hand, due to the influence of later tectonic deformation, the stratum elevation usually changes greatly. Therefore, the stratum surface generated based on elevation interpolation is quite different from the actual situation.
The cause of the aforementioned problems is that the relevant modeling methods do not utilize the development law of the loose layer. First, the loose layer is usually developed in the low-lying part of the bedrock and superimposed on the bedrock. Due to the constraints and limitations of the bedrock surface shape, from the center of the loose layer to the surrounding boundary, the relevant loose layers gradually become thinner and pinch out one by one until they are all pinched out at the contact boundary between the loose layers and the bedrock. Second, the development of loose layers follows the principle of superposition , and the thickness of each stratum is continuous. Even under the influence of later tectonic deformation (folding, faulting, and so on), the surface elevation of each stratum has changed greatly, but its thickness is still continuous.
Therefore, based on summarizing the development law of loose layers, this article proposes a 3D modeling method. First, bedrock-boundary virtual boreholes (BBVBs) are used to restrict the distribution of the loose layer. Second, pinch-out virtual boreholes (POVBs) are used to express the five stratigraphic spatial distribution modes between boreholes. Finally, assuming the boreholes in the modeling are distributed unevenly and the number is too small, a certain number of densified virtual boreholes (DVBs) can be used to keep the stratum extending continuously with a reasonable thickness. Using the measured boreholes and the aforementioned three types of virtual boreholes, the development laws of the loose layer can be applied in 3D modeling to establish a more accurate 3D model of the loose layer.
Compared with the previous modeling methods, the modeling method in this study would focus on the following three aspects for in-depth exploration: (1) It would be studied for geological modeling of a large-scale complex modeling area with alternating distributions of bedrock and loose layers. (2) It would support multiple stratigraphic spatial distribution patterns to avoid the situation in multiple strata pinch-out to the same position, which are more consistent with the actual distribution of underground strata. (3) It would interpolate the stratigraphic thickness with the stacking, avoid the intersection issue between a stratum’s top and bottom surfaces, and improve modeling accuracy.
2.1 Development law of the loose layer
First, the loose layer is mainly developed in the low-lying parts of the surface and superimposed on the bedrock. Due to the constraints and limitations of the bedrock surface morphology, the relevant loose layers gradually become thinner from the center of the loose layer to its edge and pinch out one by one until all pinch out at the contact boundary between the loose layer and the bedrock. Figure 2 shows the development characteristic of the loose layers.
In addition, under the different geological processes, the sedimentary strata show different spatial distribution forms. There are five stratigraphic spatial distribution patterns according to the genesis and spatial distribution characteristics of strata (Figure 3) . In loose layer modeling based on borehole data, it is necessary to analyze the development characteristics of borehole strata, reasonably judge the stratigraphic distribution pattern (SDP), and carry out stratigraphic correlation and modeling constraints based on the pattern law.
Finally, according to the principle of superposition, the several loose layers developed in the area are stacked together regularly from bottom to top according to the sedimentary sequence. Therefore, the 3D model of the loose layer is often called the “layer cake” model . During the development of the loose layer, due to the influence of later geological processes, surfaces of the stratum may fluctuate to varying degrees, resulting in significant changes in the surface elevation. However, the thickness of each stratum is basically the same before and after deformation (Table 1). In addition, during the actual measurement, the thickness measurement error will be smaller than the elevation measurement error, and using the thickness information to fit the surfaces of the stratum will be more accurate. Therefore, the surface of each stratum should be interpolated and fitted by the thickness information in the measured boreholes rather than the elevation information.
|Geological process type||Surface elevation||Stratum thickness||Example|
|Folding Intrusion||Continuous but inconsistent||Continuous|
|Faulting||Discontinuous and inconsistent||Continuous on the cutaway faulted block|
Therefore, in constructing the 3D model of the loose layer, we should start from the development law of the loose layer and fully consider the constraints of the bedrock surface, different SDPs, and the characteristics of stratum superposition development to construct a reasonable and high-precision 3D model.
2.2 Input data
In this article, the geological data needed to build a 3D model of loose layers mainly include vector data, measured borehole data, and digital elevation model (DEM) (ESRI (Environmental Systems Research Institute) grid format). The vector data (ESRI shapefile format) include the modeling area (user-defined polygon features) and bedrock exposed area (bedrock polygon features in the modeling area or connected to the boundary). The measured borehole data mainly include the basic information about boreholes and their strata (Figure 4).
In addition, the strata in the modeling area are uniformly numbered according to the stratigraphic sedimentary sequence . The numbering rule is “increase from the top to bottom or from young to old.” Therefore, the youngest (or the topmost) stratigraphic unit is denoted as with a stratum identifier = 1. If there are strata in the modeling area, the oldest (or the bottommost) is denoted as .
Meanwhile, to meet the conditions that all strata between boreholes can match each other, the borehole model shown in Figure 5 is used to organize the borehole strata data, which has the following two characteristics:
If a stratum is missing in the borehole, add the corresponding virtual stratum to the borehole strata (the thickness of each virtual stratum is zero, and the elevation values of the top and bottom surfaces are the same).
The number of strata in each borehole is consistent with that in the modeling area.
2.3 Modeling boundary constraints based on BBVBs
The previous 3D model construction methods of loose layers seldom consider bedrock exposure. When there is exposed bedrock in the modeling area or the modeling boundary is the bedrock burial boundary, the methods will be difficult to apply [32,33,34,35,36,37]. According to the development law of loose layers, the edge of loose layers will be constrained and limited by the surface morphology of bedrock and exhibit obvious pinch-out characteristics.
To control loose layer’s stratigraphic distribution form and pinch-out position at the bedrock boundary, virtual boreholes can be set on the contact boundary between the loose layer and the bedrock. These boreholes have the same geological information expression ability as the measured boreholes and are called BBVBs. These BBVBs can participate in determining the stratum pinch-out position and can control the stratum thickness. The location and the number of BBVBs are determined by the nodes on the bedrock contact boundary, and each node generates a BBVB.
In addition, the strata information of each BBVB refer to the first stratum with a thickness greater than 0 in the nearest measured borehole. There should be stratum whose thickness tends to be close to 0 in the BBVB, while the thickness of other strata is 0. In this article, the thickness of stratum S i is set as 0.1 m, and the ground elevation of the BBVB is obtained from the DEM. Figure 6 shows the location (Figure 6a) and strata information (Figure 6b) of the BBVBs when there is exposed bedrock in the modeling area.
2.4 Modeling support of various SDPs based on POVBs
In this article, POVBs are mainly used to control the emergence and disappearance of strata and to realize the reasonable expression of various SDPs. The processing method mainly includes the following three steps.
2.4.1 Recognition of SDPs
By analyzing the distribution of strata in boreholes, there are five SDPs (Figure 3) . The effective recognition of SDPs based on the stratigraphic development information of adjacent boreholes is a vital prerequisite for loose layer modeling.
The strategy recognizing the SDPs of any pair of adjacent measured boreholes is as follows: First, the distribution of each stratum is determined by comparing the thickness value of the same stratum in the two boreholes and classifying the strata type (missing in the right borehole, missing in the left borehole, missing in both boreholes, or existing in both boreholes). Then, the adjacent strata of the same type are merged to generate different types of stratigraphic intervals. Finally, recognize the SDP according to the type and the number of strata in each stratigraphic interval.
The specific steps of recognition of SDPS are as follows:
Defining stratigraphic distribution types. Obtain a pair of adjacent measured boreholes, obtain the thickness of each stratum in each borehole, and quantitatively judge the stratigraphic distribution types between boreholes according to the following formula (Figure 7) and save them into SL = .
where is the thickness of in borehole a on the left, is the thickness of in borehole b on the right, and is the stratigraphic distribution type between the two boreholes.
Generating stratigraphic interval . First, the stratigraphic distribution type is read from SL in turn to form a string; then, the substrings with consistent adjacent characters in the string are extracted to form a substring set; and finally, each substring set is seen as a stratigraphic interval and added to list DS = , , where means the index of , means the number of s, is the stratum corresponding to the first character of substring, is the stratum corresponding to the first character of substring, is the number of characters in the substring, and is any character in the substring.
where represents that the strata in exist continuously between two adjacent boreholes, represents that the strata in miss continuously in the left one of the two adjacent boreholes, represents that the strata in miss discontinuously in the left one of the two adjacent boreholes, represents that the strata in miss continuously in the right one of the two adjacent boreholes, represents that the strata in miss discontinuously in the right one of the two adjacent boreholes, and represents that the strata in miss continuously between two adjacent boreholes.
2.4.2 Determining the location of POVBs
Between adjacent measured boreholes, the number of POVBs and strata pinch-out positions are different according to different SDPs. The rules for determining the number of POVBs and strata pinch-out positions under different SDPs are shown in Table 2.
|SDP||Number of POVBs||Location of POVBs||Examples||Description|
|SO or SU||0 (in the DS, all s whose SDPs are SO or SU)||—||—||—|
|SR or SP||(in the DS, only one whose SDP is SR or SP)||①Equidistant insertion between adjacent boreholes||(a)||A exists whose SDP is SP with three missing strata|
|ST or SQ||1 (in the DS, only one whose SDP is ST or SQ)||(b)||A exists whose SDP is ST with one missing stratum|
|②According to the development sequence, the stratum determines the pinch-out point on each POVB from far to near|
|Multiple SDPs||Max( ) (in the DS, except that SDP is all s of SO and SU, when any two s do not exist continuously)||(c)||One whose SDP is SP and two strata are miss, one whose SDP is SQ and one stratum is missing. Due to discontinuity, two s do not need to be merged, so the number of POVBs is the maximum missing strata number|
|Max( ) (in the DS, except that SDP is all s of SO and SU, when two or more s exist continuously, the continuous s are merged to form a new DS)||(d)||Two s whose SDPs are SQ, one whose SDP is ST and one whose SDP is SR. Merge continuous s and take the maximum number of strata missing as the number of POVBs|
2.4.3 Calculating the stratum thickness in POVBs
In this article, the process of calculating the thickness of each stratum in POVBs mainly involves the following steps: (1) generating stratigraphic thickness matrix PT = , where is the index of the stratum, is the index of the POVB, and is the thickness of stratum in POVB (Figure 9); (2) calculating the thickness of each stratum in each POVB and saving in PT based on the thickness and of that in the measured boreholes on both sides; and (3) completing the ground elevation of each POVB and the top or bottom surface elevation of each stratum in this POVB according to the data structure of measured boreholes.
The methods of calculating the thickness of each stratum in the POVBs are explained in detail:
Get a POVB between boreholes and , then pick up a stratum in it.
According to Section 2.4.2, determine which the is located in and then calculate its thickness. There are two situations:
If is located in the whose SDP is SO, the thickness of is calculated using the following formula and saved into .(3)
If is located in the whose SDP is not SC, the thickness of is calculated using the following formula and saved into .(4)
where is the mark value of the in the in the borehole and is the mark value of the in the in the borehole . The method used to determine the mark value is to obtain , in which ( ) in the borehole (b), and incrementally marking from 1, which will be marked as ( ). For example, we assume that there is a between boreholes and ; if , and , then , and , and if , then .
Increase the identifier and cycle step (2) to complete the calculation of all stratigraphic thickness information in POVB .
Increase the identifier and cycle steps (1)–(3) to complete the calculation of all POVBs.
2.5 Adding DVBs based on stacking
The 3D solid model constructed directly by the measured boreholes, POVBs, and BBVBs has a rough surface, many sharp corners, and a poor visualization effect, which cannot meet practical application needs. Therefore, DVBs can be used to optimize the model and make the modeling results more refined and reasonable . The key to solving this problem is to reasonably determine the location of DVBs and the thickness of each layer.
Many associated surfaces resemble one another in nature . When a material is deposited on an old surface, the resulting new surface tends to exhibit the same or similar features as the old surface, although the new features tend to be attenuated. Therefore, the actual vertical thickness between adjacent surfaces will remain relatively stable. The process by which we take advantage of these facts to obtain improved interpretation is called stacking or conformable mapping (Figure 10). The loose layer follows the geological law of overlapping development. Controlling the stratigraphic information in DVBs based on the stacking rule is an excellent solution to overcome the aforementioned disadvantages. The determination of DVBs mainly includes two steps: location determination and stratigraphic thickness calculation.
The process of determining the DVBs is mainly as follows: (1) Take the nodes on the modeling boundary as densified points to construct DVBs; (2) discretize the interior of the modeling area by grid according to the user-defined discrete spacing ; (3) take the grid discrete points as the densified points to construct DVBs and remove the DVBs that are too close (less than the spacing threshold specified by the user) to the existing boreholes. Figure 11 shows the effect of determining the position of DVBs.
The thickness of each stratum in the DVBs can be calculated using the following steps:
Superimpose all DVBs and the TIN, and then traverse each subtriangle in the TIN in turn. If the thickness of of three boreholes corresponding to the three vertices of the subtriangle is 0, the thickness of of the DVBs inside the subtriangle will also be set to 0, and other conditions will not be handled.
Increase the identifier and cycle steps (2) and (3) to complete the calculation of all stratigraphic thickness information in each DVB.
Complete the ground elevation (obtained from DEM) of each DVB and its top and bottom surface elevations of each stratum according to the data structure of the measured boreholes.
2.6 Testing model quality
The quality inspection of the 3D loose-layer model is not only an important link to verify the advantages and disadvantages of modeling methods but also the primary method to judge the applicability of modeling methods. Aiming at the modeling method proposed in this article, the quality of the 3D loose-layer model will be tested using the following two methods:
Method Ⅰ: This test method is based on the stratigraphic thickness errors in the measured boreholes. First, a certain number of measured boreholes are picked up as the validation boreholes from the modeling area, and the remaining measured boreholes are used to construct the 3D model by the horizons-to-solids method (named M1)  based on elevation interpolates and the modeling method of this article (named M2). Then, the stratigraphic thickness value in the validation borehole (named MT, measured thickness) and the thickness value at the location of that validation borehole in the model (named ST, simulated thickness) are obtained to calculate the difference between them (named TE, thickness error). Finally, the mean and standard deviation of multiple thickness errors are calculated.
Method Ⅱ: This test method is based on the cross sections: the stratigraphic development law inside the model, an effective way to verify, can observe the profile visually. The cross sections between multiple measured boreholes can verify the accuracy of the SDPs. The cross sections on both sides of the bedrock exposed area can verify the stratigraphic gradual pinch-out distribution law.
3 Case studies
3.1 Case 1
3.1.1 Study area and dataset
The first study area is located in the Qinhuai District, Nanjing, China. The modeling area is approximately 57.11 km2, and scattered exposed bedrock areas are distributed (Figure 12a). There are 341 measured boreholes in this case, and the 3D models of some measured boreholes are shown in Figure 12b. In this case, there are seven stratigraphic units, and the strata are denoted as , , and from top to bottom.
3.1.2 Modeling performance
The modeling parameters are presented in Table 3, and the modeling results are shown in Figure 13. The model quality inspection methods proposed in Section 2.6 are used to test the modeling performance. The boreholes W120, W201, W227, and W39 (see Figure 12a) were selected as the verification boreholes, and the stratigraphic thickness error table of the two modeling methods is shown in Table 4. The mean and the standard deviation of the method of this paper are 1.23 and 1.10, respectively and all are less than the value of the comparison method, which are 3.24 and 3.23, respectively. Therefore, compared with the horizons-to-solids method using elevation interpolation, the method of this article, which interpolates stratigraphic surfaces by stacking, can fit the strata thicknesses more accurately and improve the accuracy of the model.
|Modeling parameters||DEM resolution||Discrete spacing ( d )||Spacing threshold (td)||Interpolation method||Exposed bedrock||3D data model|
|Value||5 × 5(m)||500||550||IDW||Yes||B-Rep|
|M1 (mean = 3.24, standard deviation = 3.23)||1.2||3.25||2.05||0.4||3.35||3.05||2.6||4.7||2.1||2.1||5.7||3.6|
|M2 (mean = 1.23, standard deviation = 1.10)||1.2||2.13||1.07||0.4||2.55||2.15||2.6||1.1||1.5||2.1||4.5||2.4|
MT = measured thickness, ST = simulated thickness, TE = thickness error, TE = |MT-ST|, M2 is the method of this paper, and M1 is the comparison method.
Furthermore, two stratigraphic cross sections were generated (Figure 13e and f). The pinch-out mode of the strata on both sides of the exposed bedrock area is consistent with the development characteristics of the loose layers (Figure 2), and the locations of strata pinch-out follow the SDPs. Therefore, the modeling method of this article is applicable to the area where bedrock and loose layers are alternately distributed.
3.2 Case 2
3.2.1 Study area and dataset
This case is located in the Hangkonggang District, Zhengzhou, China. The modeling area is approximately 426.5 km2 and includes 69 measured boreholes (Figure 14a). Figure 14b shows the 3D models of some measured boreholes. There is no exposed bedrock area, and the terrain is flat. There are 35 stratigraphic units in this case, and the strata are denoted as ,…, , and from top to bottom.
3.2.2 Modeling performance
The modeling parameters are presented in Table 5, and the modeling results are shown in Figure 15 (Z-axis exaggeration 20 times for clear visualization). This case selects boreholes ZG36, G128, G92, and G25–44 (Figure 13a) as the verification boreholes and then chooses seven key strata to make the stratigraphic thickness error table shown in Table 6. The mean and the standard deviation of the method of this paper are 1.30 and 0.70, respectively, and all are less than the value of the comparison method, which are 2.06 and 1.03, respectively.
|Modeling parameters||DEM resolution||Discrete spacing ( d )||Spacing threshold (td)||Interpolation method||Exposed bedrock||3D data model|
|Value||5 × 5( m )||400||400||IDW||No||B-Rep|
|M1 (mean = 2.06, standard deviation = 1.03)||5.7||4.53||1.17||2.46||4.12||1.66||0||0||0||3.45||1.64||1.81|
|M2 (mean = 1.30, standard deviation = 0.70)||5.7||4.22||1.48||2.46||4.12||1.45||0||0||0||3.45||2.33||1.12|
MT = measured thickness, ST = simulated thickness, TE = thickness error, TE = |MT-ST|, M2 is the method of this paper, and M1 is the comparison method.
In addition, because the modeling area does not contain bedrock, Method II is not carried out completely. However, Figure 15d shows that the pinch-out patterns are still valid when there are many stratigraphic units in the boreholes.
4.1 Factors affecting model quality
4.1.1 The number and distribution of the measured boreholes
The measured boreholes are the data basis for constructing the 3D loose-layer model, and their quantity and distribution are the primary factors for determining the quality of the model. In the modeling area, a large number of measured boreholes with a uniform distribution can reveal the loose layer’s distribution characteristics reasonably and ensure the accuracy of the model. In contrast, too few or unevenly distributed measured boreholes will make it difficult to accurately define the pinch-out positions of missing strata and ensure the model accuracy in the sparse area of boreholes. In China, the general principles of the measured borehole layout can refer to the national standard code for engineering geological drilling (DZ/T 0017-91).
In addition, at least one or more measured boreholes are required in the area between exposed bedrocks to control the spatial distribution of loose layers. If the area between the exposed bedrock lacks measured boreholes, the loose layer in this area will not be modeled correctly.
4.1.2 Screening of the measured boreholes
If boreholes that are too shallow are involved in constructing the 3D loose-layer model, the illusion of a missing local stratum will appear. Therefore, it is necessary to screen the shallow boreholes to ensure that the borehole depth is basically the same and to meet the requirements of 3D modeling. According to the bottom stratum ID and the depth of each borehole, the shallow borehole can be automatically screened.
4.1.3 Discrete spacing
The number of DVBs directly determines the smoothness of the model. If the discrete spacing is set too large, the number of DVBs is small. For the model established in this case, the strata surfaces are rough, the boundary edges and corners are sharp, and it has a poor visualization effect. If the discrete spacing is set to small, the number of DVBs will be large. In this case, the strata surfaces will be smooth, and the boundary transition will be gentle.
However, it is meaningless to pursue the model’s accuracy further when the fluctuation of the local level is indistinguishable to the human eye, which will reduce the efficiency of model construction, increase the amount of model data, and impact the execution efficiency of model display and 3D spatial analysis applications. Therefore, the appropriate discrete spacing can achieve a balance among the visualization effect, modeling efficiency, and the amount of model data.
4.1.4 DEM accuracy
During the construction of a 3D loose-layer model, the ground elevation of virtual boreholes needs to be obtained from the DEM to ensure that the surface of the topmost stratum is consistent with the DEM. Therefore, the DEM accuracy has an important influence on the modeling accuracy of the 3D loose-layer model.
First, the DEM used in this article can be generated by interpolating the ground elevation of the measured boreholes. In this case, the constructed model surface is relatively smooth, and the surface fluctuation is continuous. The quality of the DEM entirely depends on the number and density of measured boreholes, and there are two deficiencies: (1) When the number of measured boreholes is small or unevenly distributed, the constructed DEM is rough, and the surface fluctuation is not apparent. Especially in the area where a large number of measured boreholes are missing, there is a significant error between the constructed DEM and the actual ground. (2) If the drilling time of the measured boreholes is too long, their ground elevation may be inconsistent with the actual surface elevation, which makes the current situation of the constructed surface DEM poor, and it is difficult to realize the seamless integration of the constructed 3D model and the latest DEM.
Second, the latest high-precision DEM is used for modeling. In this case, the advantage is that it can easily realize the seamless integration of the 3D loose-layer model and the latest DEM. The disadvantage is that the ground elevation of the measured borehole is inconsistent with that of the DEM and needs further optimization.
The actual drilling time is a reference standard for choosing which DEM to use. If the time is too early, the first DEM is suitable. In contrast, the time is very close, and the second DEM is acceptable. In addition, if the 3D loose-layer model does not need to consider the actual surface fluctuation and is mainly used for excavation analysis and bearing capacity calculation, using the first DEM is more appropriate. In contrast, if the 3D model is mainly used for the above-ground and underground 3D scene visualization, the second DEM is recommended to ensure the effective integration between underground, ground, and above-ground 3D models.
4.2 Applicability of the method
4.2.1 Applicability of modeling objects
First, the modeling method of this article is designed for a large-scale modeling area with alternating distributions of bedrock and loose layers. Starting from the development law of the loose layer, the method fully considers the five stratum distribution modes of the loose layer.
Second, a pure development area of a loose layer is also applicable. When the modeling area does not contain bedrock, there is no need to provide bedrock boundary data in the modeling process.
Finally, the method of this article is also applicable for 3D modeling in bedrock areas with relatively simple geological structures. For example:
Both horizontal and monoclinic structures have prominent overlapping distribution characteristics, similar to the distribution mode of loose strata. Figure 16 shows a 3D model of a bedrock area that contains six kinds of rocks. The strata of this area are stacked horizontally and can be divided into 14 stratigraphic units according to the stratigraphic sequence.
Vertical fold, inclined fold, and other folds without stratigraphic inversion. For those geological structures, their rock mass is in a continuous distribution state, which is very suitable for using superposition technology to calculate stratigraphic thickness.
4.2.2 Applicability of modeling area
When a modeling area is too large, the modeling efficiency will decline and generate a model file with an excessive amount of data. In this case, the modeling area can be divided into multiple modeling units using the administrative units or geomorphic units and then modeling unit by unit. In addition, if there are islands or gaps in the modeling area (Figure 17), the modeling area can also be dispersed into multiple modeling units for modeling unit by unit.
4.2.3 Applicability of modeling data
In this article, the measured boreholes in the modeling area need to meet certain constraints to participate in the modeling method.
In a modeling area, the numbering rule of the strata identifiers of the measured boreholes should be consistent.
The direction of the measured boreholes should be vertical. The stratum thickness in the inclined borehole is inconsistent with the actual stratum thickness. The inclined boreholes need to be processed as vertical boreholes before participating in modeling.
When the measured boreholes are insufficient, the virtual boreholes can be generated based on the geological profile and participate in the 3D modeling together with the measured boreholes.
4.3 Relationships between the proposed method and existing methods
There are six basic categories of methods currently existed for 3D geological modeling as discussed in Section 1. Both the proposed method in this study and the existing methods use borehole data, but from different perspectives to construct 3D models.
First, the existing methods are mainly based on the elevation information to interpolate the stratigraphic surfaces, while the proposed method is based on stratigraphic thickness to fit surfaces. They just use different interpolation strategies and do not conflict. In other words, we provided a new stratigraphic interpolation strategy for users to choose. For the superficial loose sedimentary layer, the accuracy difference of 3D models constructed by various methods is slight. However, for the loose-layer areas affected by tectonic action, the proposed method is more accurate in fitting the stratigraphic thickness.
Then, the existing methods are mainly applied to the geological modeling of cities in plain or coastal areas, where there is usually no bedrock exposed. This study’s method is suitable in large-scale urban areas where loose layers and bedrock are distributed alternately.
Finally, for the 3D model constructed by this study, its top surface is controlled by ground DEM, and the bedrock surface constrains the boundary. Compared with the 3D models constructed by the existing methods, it could be implemented more easily to realize seamless integration with bedrock geological bodies and ground-building models.
This article proposes a 3D solid model construction method for loose layers based on the stratum development law. The proposed method adds three different virtual boreholes, BBVBs, POVBs, and DVBs, to realize the gradual pinch-out of the loose layers near the exposed bedrock area, support different SDPs, and control the thickness change of each stratum, respectively. This method could address the shortcomings, such as the lack of bedrock surface constraints, superficial strata pinch-out processing, and the higher fitting error of the strata surface. Especially for a modeling area with significant topographic fluctuation, the method can accurately describe the morphological characteristics of the strata.
Case studies demonstrated the effectiveness of this 3D loose-layer modeling method in the large-scale exposed bedrock area and the pure loose-layer area. Compared to the previous method with elevation interpolation, the stratigraphic thickness fitted by the proposed method based on stacking is more accurate. The stratigraphic thickness errors’ mean and standard deviation of case 1 decreased by 2.11 and 2.13 m, and they in case 2 declined by 0.96 and 0.33 m. Moreover, in the discussion, we used the proposed method to construct the 3D model of bedrock, which also achieved a good effect.
In conclusion, the proposed method is suitable for geological modeling in large-scale urban areas where loose layers and bedrock are distributed alternately. The constructed 3D models follow the stratum development law, ensuring the effective integration between underground, ground, and above-ground 3D models. This study’s results help evaluate the megacity’s underground space carrying capacity and the development and utilization of urban underground space.
Funding information: This study was supported by the National Key R&D Program of China (Project No. 2021YFE0112300) and the National Natural Science Foundation of China (Project Nos. 41971068 and 41771431).
Author contributions: Yangen Shen performed the study, analyzed the data, developed the proposed algorithm, and drafted the original manuscript. Anbo Li managed the entire research project and offered supervision. Jianchu Huang developed of the partial modules of the prototype system. Guonian Lv and Kaili Li supported in collecting the experimental data and in the way of writing.
Conflict of interest: The authors declare no conflict of interest.
 Guo T, Wang L, Wang M, Dai W, Wu X, Li L, et al. Three-dimensional geological modeling and spatial analysis from geotechnical borehole data using an implicit surface and marching tetrahedra algorithm. Eng Geol. 2021;284:106047. 10.1016/j.enggeo.2021.106047.Search in Google Scholar
 Collon-Drouailler P, Steckiewicz-Laurent W, Pellerin J, Laurent G, Caumon G, Reichart G, et al. 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France. Comput Geosci. 2015;77:29–43. 10.1016/j.cageo.2015.01.009.Search in Google Scholar
 Zanchi A, Francesca S, Stefano Z, Simone S, Graziano G. 3D reconstruction of complex geological bodies: Examples from the alps. Comput Geosci. 2009;35(1):49–69. 10.1016/j.cageo.2007.09.003.Search in Google Scholar
 Caumon G, Collon-Drouaillet PL, Le Carlier de Veslud C, Viseur S, Sausse J. Surface based 3D modeling of geological structures. Math Geol. 2009;41:927–45. 10.1007/s11004-009-9244-2.Search in Google Scholar
 Zhao B, Hua H, Chen X, Liang D, Liu P, Liu G. New method for estimating strike and dip based on structural expansion orientation for 3D geological modeling. Front Earth Sci. 2021;15:676–91. 10.1007/s11707-021-0903-z.Search in Google Scholar
 Wu X, Liu G, Weng Z, Tian Y, Zhang Z, Li Y, et al. Constructing 3D geological models based on large-scale geological maps. Open Geosci. 2021;13(1):851–66. 10.1515/geo-2020-0270.Search in Google Scholar
 Aygar EB, Gokceoglu C. Analytical solutions and 3D numerical analyses of a shallow tunnel excavated in weak ground: a case from Turkey. Int J Geo-Eng. 2021;12(1):1–26. 10.1186/s40703-021-00142-7.Search in Google Scholar
 Liu XY, Li AB, Chen H, Men YQ, Huang YL. 3D Modeling Method for Dome Structure Using Digital Geological Map and DEM. ISPRS Int J Geo-Inf. 2022;11(6):339. 10.3390/ijgi11060339.Search in Google Scholar
 Hao M, Wang H, Deng C, He W, Zhang L, Xue J, et al. 3D geological modeling and visualization of above-ground and undergro-und integration — taking the Unicorn Island in Tianfu new area as an example. Earth Sci Inform. 2019;12:465–74. 10.1007/s12145-019-00394-z.Search in Google Scholar
 Wycisk P, Hubert T, Gossel W, Neumann C. High-resolution 3D spatial modelling of complex geological structures for an environmental risk assessment of abundant mining and industrial megasites. Comput Geosci. 2009;35(1):165–82. 10.1016/j.cageo.2007.09.001.Search in Google Scholar
 Olivier K, Thierry M. 3D geological modeling from boreholes cross-sections and geo-logical maps application over former natural gas storages in coal mines. Comput Geosci. 2008;34(3):278–90. 10.1016/j.cageo. 2007.09.005.Search in Google Scholar
 Song R, Qin X, Tao Y, Wang X, Yin B, Wang Y, et al. A semi-automatic method for 3D modeling and visualizing complex geological bodies. Bull Eng Geol Environ. 2019;78(3):1371–83. 10.1007/s10064-018-1244-3.Search in Google Scholar
 Hou S, Yang L, Deng C, Ye J, Clarke K, Yang J, et al. Assessing quality of urban underground spaces by coupling 3D geological models: The case study of Foshan city. South China. Comput Geosci. 2016;89:1–11. 10.1016/j.cageo.2015.07.016.Search in Google Scholar
 He H, He J, Xiao Z, Zhou X, Liu Y, Li C. 3D geological modeling and engineering properties of shallow superficial deposits: A case study in Beijing, China. Tunn Undergr Sp Technol. 2020;100:103390. 10.1016/j.tust. 2020.103390.Search in Google Scholar
 Tonini A, Guastaldi E, Massa G, Conti P. 3D geo-mapping based on surface data for preliminary study of underground works: a case study in Val Topina (Central Italy). Eng Geol. 2008;99(1):61–9. 10.1016/j.enggeo.2008.02.010.Search in Google Scholar
 De Rienzo F, Oreste P, Pelizza S. Subsurface geological-geotechnical modelling to sustain underground civil planning. Eng Geol. 2008;96(3–4):187–204. 10.1016/j.enggeo.2007.11.002.Search in Google Scholar
 Chen J, Zhou F, Guo W. Study on the relationship between subsidence coefficient and mining degree under a thick alluvium stratum. J Min Saf Eng. 2012;29:250–4 (in Chinese). 10.3969/j.issn.1673-3363.2012.02.018.Search in Google Scholar
 Chang Z, Wang J, Chen M, Ao Z, Yao Q. A novel ground surface subsidence prediction model for sub-critical mining in the geological condition of a thick alluvium layer. Front Earth Sci. 2018;9:330–41. 10.1007/s11707-014-0467-2.Search in Google Scholar
 Nichols G. Sedimentology and stratigraphy. 2nd edn. Hoboken, New Jersey, United States: John Wiley & Sons; 2009.Search in Google Scholar
 Dong M, Neukum C, Hu H, Azzam R. Real 3D geotechnical modeling in engineering geology: a case study from the inner city of Aachen, Germany. Bull Eng Geol Environ. 2015;74(2):281–300. 10.1007/s10064-014-0640-6.Search in Google Scholar
 Bozzano F, Andreucci A, Gaeta M, Salucci R. A geological model of the buried Tiber River valley beneath the historical centre of Rome. Bull Eng Geol Environ. 2000;59(1):1–21. 10.1007/s100640000051.Search in Google Scholar
 Alemdag S, Gurocak Z, Cevik A, Cabalar AF, Gokceoglu C. Modeling deformation modulus of a stratified sedimentary rock mass using neural network, fuzzy inference and genetic programming. Eng Geol. 2016;203:70–82. 10.1016/j.enggeo.2015.12.002.Search in Google Scholar
 Bahaaddini M, Hosseinpour Moghadam E. Evaluation of empirical approaches in estimating the deformation modulus of rock masses. Bull Eng Geol Environ. 2019;78(5):3493–507. 10.1007/s10064-018-1347-x.Search in Google Scholar
 Aygar EB, Karahan S, Gullu S, Gokceoglu C. Analytical and numerical analyses of the support system for a large-span tunnel in challenging and seismically active ground conditions. Transp Infrastruct Geotech. 2022;1–44. 10.1007/s40515-022-00251-5.Search in Google Scholar
 Saceanu MC, Paluszny A, Zimmerman RW, Ivars DM. Fracture growth leading to mechanical spalling around deposition boreholes of an underground nuclear waste repository. Int J Rock Mech Min Sci. 2022;152:105038. 10.1016/j.ijrmms.2022.105038.Search in Google Scholar
 Nonogaki S, Masumoto S, Nemoto T, Nakazawa T. Voxel modeling of geotechnical characteristics in an urban area by natural neighbor interpolation using a large number of borehole logs. Earth Sci Inform. 2021;14(2):871–82. 10.1007/s12145-021-00600-x.Search in Google Scholar
 Zhu F, Zhang J, Li J, Pan X, Sun Z. Building 3D solid models of sedimentary stratigraphic systems from borehole data: an automatic method and case studies. Eng Geol. 2012;127(1):1–13. 10.1016/j.enggeo.2011.12.001.Search in Google Scholar
 Zhu F, Wu C, Liu G, Shang G. Reconstruction of 3D strata model based on borehole data. Geogr Geo Inf Sci. 2004;20(3):26–30 (in Chinese).Search in Google Scholar
 Zhang Y, Bai W. An approach of 3D stratum modeling based on tri-Prism volume elements. J Image Graph. 2001;6:285–90 (in Chinese).Search in Google Scholar
 Guo J, Wang Z, Li C, Li F, Jessell MW, Wu L, et al. Multiple-point geostatistics-based three-dimensional automatic geological modeling and uncertainty analysis for borehole data. Nat Resour Res. 2022;31:2347–67. 10.1007/s11053-022-10071-6.Search in Google Scholar
 Bi Z, Wu X, Li Z, Chang D, Yong X. DeepISMNet: three-dimensional implicit structural modeling with convolutional neural network. Geosci Model Dev. 2022;15(17):6841–61. 10.5194/gmd-15-6841-2022.Search in Google Scholar
 Guo J, Zhang R, Wu L, Yang Y. An automatic method for generating curvilinear geological section with stratal pinch-out. 2011 IEEE International Geoscience and Remote Sensing Symposium. IEEE; 2011. p. 2965–8. 10.1109/IGARSS.2011.6049838.Search in Google Scholar
 Steno N. The prodromus of Nicolaus Steno’s dissertation concerning a solid body enclosed by process of nature within a solid. London: University of Michigan Press; 1916.10.5962/bhl.title.54340Search in Google Scholar
 Walton D. Ground water resource evaluation. New York: McGraw-Hill; 1970.Search in Google Scholar
 Zhang F, Zhu H, Ning X. Modeling method of 3D strata suitable for massive data. J of Geo-Inf Sci. 2006;25(sup.1):3305–10 (in Chinese).Search in Google Scholar
 Tearpock J, Bischke E. Applied subsurface geological mapping with structural methods. New Jersey: Upper Saddle River; 2002.Search in Google Scholar
 Wang Z, Ma J. Layer-constrained triangulated irregular network algorithm based on ground penetrating radar data and its application. J Beijing Inst Technol. 2018;27:150–8. 10.15918/j.jbit1004-0579.201827.0118.Search in Google Scholar
 Shepard D. A two-dimensional interpolation function for irregularly-spaced data. Proceedings of the 1968 23rd ACM National Conference; 1968. p. 517–24 10.1145/800186.810616.Search in Google Scholar
 Deutsch V, Journel G GSLIB geostatistical software library and user’s guide. New York: Oxford University Press; 1998.Search in Google Scholar
© 2022 the author(s), published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.