© S. Sun et al., Published by EDP Sciences 2020

1 Introduction

After a decade of initial preparation and industrial testing (Zou et al., 2015), China’s shale gas industry reported the production of 7.8 billion m³ in 2016, making it third in world production, with a cumulative production of up to 12 billion m³. The Sichuan Basin in the Middle-Upper Yangtze region is a common location for shale gas exploration in China. In this area, a large amount of shale gas is being recovered from the stratigraphic section extending from the Upper Ordovician Wufeng Formation (O₃w) to the Silurian LIandovery Longmaxi Formation (S₁l) (Zou et al., 2015). Currently, the exploration and development of shale gas in this basin has targeted sub-layers. After several years of efforts, the industry has determined that the WF2-LM5 graptolite zone in the O₃w–S₁l is a favorable target for shale gas development (Chen et al., 2015).

Many scholars have determined the Wufeng Formation (WFF) shale, located in the bottom of the WF2-LM5 targeting zone, has extremely heterogeneous reservoir properties and thicknesses in the lateral and vertical directions, a favorable reservoir for accmulating shale gas (Zou et al., 2015; Wang et al., 2016; Yang et al., 2016; Zou et al., 2016; Chen et al., 2017; Nie et al., 2017). WFF shale presented a high maturity, with Ro value of 2.0 ∼ 3.5%, showing a high total porosity of 4.6 ∼ 8.2%, and the permeability is 0.0002 × 10⁻³ ∼ 0.5000 × 10⁻³ µm², developing rich organic matter pores and the micro-cracks. During the early sedimentary period, with the warm and humid paleo-climate, a large area of undersea oxygen-deficient environment was formed due to sea-level rising. With the values of δ¹³C values of −30.2 ∼ −29.9‰ and P₂O₅/TiO₂ of 0.24 indicating a nutrient-rich surface water, the productivity of plankton (e.g., algae, radiolarian, and graptolite) was high, and complexes of bioclastic particles, organic matter and clay minerals slowly settled in a “marine snow” (Algeo and Lyons, 2006; Hammes et al., 2011) manner, forming organic- and biogenic silica-rich argillaceous, siliceous shale with a TOC of 2.0 to 8.0%. During the middle-late sedimentary periods of the WFF, due to sea level fall (about 50 ∼ 100 m), graptolites feeded on plankton became extinct, δ¹³C began to exhibit a positive drift, reaching −29‰ to −27.6‰, P₂O₅/TiO₂ value reached a peak high of 0.84 where the oxygen-deficient deep water area was narrowed to the depression zone in the southern-eastern-northeastern Sichuan regions, with the TOC of 2.7 ∼ 8.4%.

However, its quantitative distribution and paleo-geographic setting are still unknown, which impedes decision-making in exploration. Traditional paleo-geographic studies of the WFF in southern China have mainly focused on investigating the lithofacies and biofacies, e.g., the lithofacies paleo-geographic map of (Feng et al., 2001) and the lithofacies and biofacies distribution map of Chen et al. (2004). However, these maps are less satisfactory in view of precision and quality due to insufficient original data points, the subjectivity of the thickness contour mapping using traditional qualitative and artificial interpolation methods, and inadequate spatial distribution analysis. Currently, the seismic data available is not sufficient for the demand of the precision exploration of shale gas. At present, there is no reliable method for large-scale high-precision quantitative assessment of shale gas.

With the rapid development of computer technology, the combination of databases and Geographic Information System (GIS) methods provides a new solution for conducting paleo-geographic studies (Fan et al., 2016). The database enables the acquisition, storage, and management of massive amounts of spatial data, and the GIS method provides a spatial interpolation that can predict the attributes of an unknown point in the area of interest from a known point through extrapolation or interpolation. The combination of these tools allows for the effective and objective utilization of massive amounts of data for determining the law and provides a powerful tool for the quantitative restoration of the distribution pattern of the shale.

In this paper, taking the WFF shale in the Upper Yangtze region in southern China as an example, we reconstructed the spatial-temporal lithofacies distribution using ArcGIS and data for 310 sections (including drilling) collected from the Geobiodiversity Database (GBDB). Based on a detailed comparison, we selected the optimal interpolation method from three common methods to quantify the spatial distribution of the WFF shale. In addition, we extracted data from 275 sections to determine the paleo-water-depth and paleo-topography and reconstructed the 2D and 3D paleo-topographic framework using the Ordinary Kriging method of ArcGIS. Our results provide the basis for the precise exploration and development of shale gas in the WFF, and therefore provide realistic theoretical significance.

With the development and deepening of geoscience research, geology research has gradually entered the era of “big data”. Facing with the massive geological data, the methodology of GIS technology is readily exportable to any resource play having spatial distribution pattern. Except for paleo-geography reconstruction, strata thickness, rock mineral composition, TOC, vitrinite reflectance and other geochemical parameters can be applied with this means for quantitative research. The results can be generated from geological data providing a rapid, efficient quantitative means to assess.

2 Geologic setting

The Upper Yangtze region is bordered to the north by the fault on the southern margin of the Qinling, to the south by the Yadu-Ziyun-Luodian fault, to the west by the Longmenshan fault system, and to the east by the Xuefeng Mountains. The region covers a total area of about 3.5 × 10⁵ km² (Deng, 2013). The Sichuan Basin being the main part of the study area, but it also contains the northern Yunnan-Guizhou, eastern Chongqing-eastern Hunan-Hubei, and northern Guizhou-Central Guizhou regions around the basin.

The Yangtze block was located at a tropical to sub-tropical zone on the western margin of Gondwanaland during the Late Ordovician. In addition to the uplift and expansion of the ancient lands and paleo-uplift, the Upper Yangtze Sea closed. The Yangtze platform transitioned from a shallower in the west and deeper in the east structure to a post-uplift, semi-restricted, neritic marine basin, which was surrounded to the west, south, and north by three uplifted areas, i.e., the Central Sichuan uplift, the northern Sichuan uplift, and the Central Guizhou-Xuefeng uplift, and dipped gently to the northwest, forming a depression (Mou et al., 2014). Within the basin, the distribution of the sedimentary facies belt and the provenance was controlled by these uplifts. The closure to the south and the connection to the Qinling Ocean allowed for the formation of a semi-closed neritic marine basin, causing the WFF shale to be distributed over a broad area (Mu, 1954; Rong and Chen, 1987; Feng et al., 2001) .

The Late Ordovician WFF was named by Sun Yunzhu in 1931 based on the graptolite sample recovered in Wangjiawan section, Yuyangguan Town, Hubei Province (Mu, 1954). Later, over nearly a decade of studies of shale gas exploration and development practices produced a massive amount of stratigraphic, paleontology, paleo-geography, and geochemistry data. In the Upper Yangtze region, the WFF was deposited during the Ordovician, from the Late Katian to the Early Hirnantian, i.e., over 3.19 Ma from 444.43 Ma to 447.62 Ma. This process led to the development of four graptolite zones in two stages (Chen et al., 2015) (Fig. 1), i.e., from bottom to top, the Dicellograptus complanatus, Dicellograptus complexus, and Paraorthograptus pacificus in the upper section of the Katian stage, and Normalograptus extraordinariu at the base of the Hirnantian stage. The WFF is underlain by the Linxiang Formation and separated from the overlying Longmaxi Formation by a thin-bedded shelly facies argillaceous limestone, known as the Guanyinqiao Layer.

Some formations which were deposited synchronously but with different lithofacies than the WFF shale include the Daduhe, Tiezufeike, Nanzheng, Tianmashan, and Tianlingkou formations (Fan et al., 2016) (Fig. 1). The WFF shale is distributed across a broad area within the Yangtze platform and ranges from 6 to 20 m thick. It is composed of grey black to black to yellowish-brown weathered, thin-bedded shale, siliceous rocks, and carbonaceous and siliceous shales, intercalated with multiple volcanic tuff layers. This formation contains abundant graptolites, a few phosphatic shell brachiopods, trilobite and brachiopods present at the top, and radiolarians in the siliceous rocks. The Daduhe Formation consists mainly of dolostone interbedded with shale and is present in northeastern Yunnan and southwestern Sichuan (Tang, 2017). The Tiezufeike Formation is distributed in Butuojinyang, western Sichuan (Tang, 2017), with the lithology of grey black dolomitic limestone and dolostone. The Dajing Formation consists mainly of dolomitic limestone and limy dolostone and is distributed in northeastern Yunnan and southern Sichuan (Tang, 2017). The Nanzheng Formation, which was deposited in northern Sichuan, and the Tianmashan Formation, which was deposited in southern Hunan, are mainly composed of clastic rocks. The basal interface of the WFF shale is roughly isochronous, while the top interface is partially diachronous.

Fig. 1 Division and correlation of the Ordovician strata in the Upper Yangtze region (Zhang et al., 2014a). (A) the broad area in the Upper Yangtze region containing the typical black shale of the Wufeng Formation; (B) the area extending from northern Sichuan to southern Shanxi; (C) the area extending from southeastern Sichuan to northeastern Yunnan; (D) the area extending from southwestern Hunan to northeastern Guangxi.

3 Data and programs

The combination of GBDB and GIS allows for the large-scale extraction of applicable data, the effective, high-precision, large-scale, quantitative identification of the spatial distribution pattern and law, and the performance of prediction analysis. This method involves two steps: data preparation and data modeling.

3.1 Data preparation

High-precision modeling of the spatial and temporal distribution of strata across a broad area requires a massive data set. Large-scale stratigraphic and paleontologic databases, such as the GBDB (Fan et al., 2013), have been compiled in recent years and provide a useful tool for the investigation, acquisition, editing, management, downloading, and application of large-scale, standard and quantitative data. The GBDB was established in 2006 by a research team at the Nanjing Institute of Geology and Paleontology, Chinese Academy of Sciences. The GBDB provides a section-focused digital scientific research platform for paleontologic and stratigraphic studies based on internet, databases, and GIS, and has become the world’s largest stratigraphic database and the second largest paleontologic database (Fan et al., 2016; Zhang et al., 2016).

The data used in this paper was mainly collected from scientific articles, regional geological reports, geological books, the GBDB, and core observations. By searching, collecting, preparing, and analyzing data from the GBDB, a total of 257 section points in the WFF shale and synchronous strata in the Upper Yangtze region were selected and exported from the GBDB in the model-required format. Using the shale gas drilling results in southern China, based on core observations and data collection, 18 wells were determined to provide a control on the thickness of the graptolite zone. In addition, due to a lack of data in the GBDB on the Daduhe and Dajing formations on the western margin of the Yangtze region, 35 stratigraphic sections were collected from the literature (Tang, 2017). In total, 310 data points were collected from 7 provinces/municipalities, i.e., Sichuan, Chongqing, Yunnan, Guizhou, Hubei, Hunan, and Shaanxi (Fig. 2). Of these data points, 179 are from the WFF shale, 60 are from the other synchronous strata, and 71 indicate the absence of the formation. Statistics reveal that the thickness of the WFF ranges from 0 to 40 m with an average of 4.1 m.

The sections used in this study can be divided into three types: (1) sections including the WFF black shale; (2) sections including formations deposited synchronously with the WFF shale, but with different lithofacies; and (3) sections that do not include sediments similar to the WFF shale. The data points for the WFF and synchronous formations can be utilized to determine the distribution range of the different lithofacies, while the data points for the absence of these formations provide an opportunity to define the marine-continent boundary at the time and the range of the area that initially received sediments, but was later eroded (Fan et al., 2016). The object data for interpolation includes: (1) the geographic position, i.e., longitude and latitude (must), elevation, and the administrative dividing unit; (2) the attributes, i.e., the name of the litho-stratigraphic unit (or absence), the formation thickness, the analytical and testing parameters, and other discrete spatial data; and (3) remarks, auxiliary instructions and data sources. In general, the larger the volume of the object data point used for interpolation, the more uniform its spatial distribution and the better the effect of the interpolation.

Fig. 2 Geographic distribution of the section points of the Wufeng Formation shale and the synchronous strata in the Upper Yangtze region (The yellow circles represent the Wufeng Formation and the synchronous strata points; the orange crosses represent the sections that do not contain these formations; the brown line outlines the Sichuan Basin; and the red line outlines the Yangtze region).

3.2 Data modeling

GIS provides a computer system for the input, storage, query, management, computation, analysis, displacement, and description of the geographic distribution data in land surface space (Chen et al., 1999; Tsung, 1999), which enables the spatial storage and rapid mapping of the geological data volume as well as spatial data analysis and modeling. GIS has been applied frequently in multiple fields, such as mining, environment, hydrology, geochemistry, mobile location service, and geological modeling.

In essence, GIS-based spatial interpolation is theoretically based on the ‘similarity and proximity’ law, and has proven effective in fitting the function of the spatial position and the attribute value using limited discrete data to predict the attribute features at an unknown point within a geological body and to explore the spatial distribution model and law (Zhang et al., 2016). Compared with traditional practices that determine the contour using artificial interpolation, GIS provides enhanced objectivity, quantification, and analyzability, making it a powerful tool for determining the quantitative and precise spatial and temporal stratigraphic distribution of massive amounts of section data that contains geographic position information and for analyzing this data (Zhang et al., 2016).

Objectivity means that the interpolation result is given by the computation as a mathematical equation. This process is fast and precise and is less affected by human subjectivity. For the same interpolation algorithm and parameters, different operators may yield the same interpolation result, demonstrating the repeatability and replicability of the method. Precise quantification means that the maps prepared based on the massive amount of data are high resolution and high precision in time and space, which allows for the precise identification of the interface between different lithofacies. Analyzability refers to the quantification of the spatial position and attribute, the rational derivation of the unknowns, and the subsequent analysis. For example, it is possible to compute the distribution area and volume, and even the variation ratio of the thickness along a certain direction. In addition, it allows for the prediction of shale gas sweet spots and resources using the total organic carbon (TOC) and porosity of the shale and the increases the precision and depth of any related studies.

In this study, we reconstructed the geographic lithofacies, analyzed and inspected the thickness distribution, and conducted spatial topographical modeling of the WFF using Arc GIS, a software of ESRI company. The quality of the spatial interpolation results depends on the data quality and the interpolation method. If the data quality is fixed, the software provides a variety of interpolation methods. In this study, we selected the optimal interpolation method based on the actual data conditions of the reconstruction of the lithofacies distribution and the related analysis.

4 Reconstruction analysis

4.1 Reconstruction of the lithofacies distribution

A total of 310 sections (containing 8 litho-stratigraphic units) were collected to identify the geographic distribution of the WFF in the Upper Yangtze region, including 179 sections (including 18 drilling sections) from the WFF, 60 sections from synchronous strata with different lithofacies, and 71 sections that did not contain synchronous strata (Tab. 1).

These sections were projected onto the administration map of the Upper Yangtze region, which based on their lithology and paleontological properties, can be divided into four types (Fan et al., 2016) (Fig. 3).

Fig. 3 Geographic distribution of the Wufeng Formation and synchronous formations with different lithofacies in the Middle-Upper Yangtze region.

4.1.1 Type A

This type of section contains typical WFF black shale, which contains Late Ordovician graptolite. In Figure 3, type A1 is the calcareous shale facies zone with a high carbonate content, and type A2 is the siliceous shale facies zone with a high siliceous mineral content.

4.1.2 Type B

This type of section contains synchronous carbonate rocks. In Figure 3, type B1 is the Nanzheng Formation shale interbedded with limestone in northern Sichuan, which consists of mixed sediments containing graptolite, brachiopods, and siphonopods (Fan et al., 2016), and represents a coastal tidal flat facies environment. Type B2 is the tidal flat facies black shale interbedded with limestone and intercalated with the limestone of the Daduhe and Tiezufeike formations in the western Sichuan region (Deng, 2013).

4.1.3 Type C

This type of section shows the presence contains synchronous clastic rocks, which consist of the flysch sediments of the Tianmashan and Tianlingkou formations, which allowed for the development of (silty) sandy shale.

4.1.4 Type D

This type of section is mainly located in the Hannan, central Sichuan, and central Guizhou ancient lands on the periphery of the deposition zones of the WFF and other synchronous formations and in some of the regions of the Upper Yangtze platform, such as the Hunan-western Hubei underwater highland (Chen et al., 2001) on the Hubei-Hunan border and several isolated regions extending from northern Guizhou to southern Chongqing (Rong et al., 2011).

In the Upper Yangtze region, the WFF black graptolite shale is distributed broadly throughout the Yangtze platform, while the other synchronous formations with different lithofacies are distributed along the margin. The shale zone is bordered to the north by the Qinling-Dabieshan orogenic belt and the Tanlu fault and to the west, south, and southeast, by either the ancient land or the absence of the formation. Geographically, it extends across 7 provinces/municipalities, including eastern and western Sichuan, southernmost Shaanxi, eastern Chongqing, northeastern Yunnan, northern Guizhou, the main part of Hubei, and central and northern Hunan. These areas were in a normal neritic epicontinental sea environment, with the exception of the hundreds of meters of flysch sediments in the Tianmashan Formation clastic rock in southern Hunan, which represent the rapid subsidence of the basement (Ge, 2013). In addition, sections lacking the formation can be utilized to precisely define the marine-continent boundary at the time.

4.2 Restoration of the formation thickness distribution

Based on the data set containing 310 sections collected to study geographic distribution, a sub-set comprising 248 sections (including 177 data points for the WFF and synchronous formations, and 71 lacking the formation) was extracted, in which each data point had a precise thickness. To determine the thickness distribution of the WFF shale, it was necessary to set the thickness of the synchronous formations with different facies and the formation absence point section to 0 in the interpolation and to ensure that the thickness of the WFF remained constant. ArcGIS provides a variety of interpolation methods, each of which has its own applicability and limitations. In order to select the best method for the interpolation of the thickness of the WFF in the Upper Yangtze region, three interpolation methods were applied to the geological modeling and interpolation, i.e., the Inverse Distance Weighting, Kriging, and Radial Basis Function methods. The optimal method was selected through general verification and was used for the interpolation and extrapolation.

4.2.1 Interpolation methods

4.2.1.1 Inverse Distance Weighting (IDW)

The IDW method is a precise interpolation method based on the “similarity and proximity” law, which is the first law of geography (Tobler, 1970). This method is simple to understand. The closer the known point is to the interpolating point, the higher the weight. IDW considers the spatial self-correlation, which is directly dependent upon the power of the distance and the hunting zone of the neighborhood. The variation in distance changes the weight of the unknown, and the hunting zone of the neighborhood controls the number and mode of the sample points. This method is relatively applicable to cases in which the sample data is closely spaced and is distributed uniformly. The IDW method yields interpolation results ranging between the maximum and minimum values. However, the weight, uniformity, and abnormal data may result in the occurrence of the “cow-eye” effect, i.e., the appearance of small and closed contour rings.

4.2.1.2 Kriging

The Kriging method is a statistical interpolation method that describes structure and randomness. It utilizes a semi-variogram to conduct an unbiased linear optimal estimation of the regional variation (Sun, 1990). The Kriging method takes into consideration the correlation of the prediction and the known points and the spatial self-correlation and random variation of the variable. Error analysis is also possible with this method. There are various types of Kriging, which can be applied to different objects and statistical data characteristics. The Ordinary Kriging method, which is the most widely used, is adopted in our study. Although the Kriging interpolation does not produce edge effects, it is a local optimal estimation method that pays less attention to the overall spatial correlation. Therefore, this method is applicable to cases in which there is a high spatial self-correlation and the sample points are closely distributed. The statistical probability-based interpolation result may be greater or less than the extreme value of the sample points.

4.2.1.3 Radial Basis Function (RBF)

The RBF method is also known as the radial basis neural network, which is a precise interpolation method for predicting values greater than the maximum measurement and less than the minimum measurement. In essence, it is an interpolation method that uses stochastic approximation technology, which applies the recursion to the back propagation learning algorithm, to create and predict a smooth surface with a large number of known sample data points on the basis of smoothness. The RBF method is considered to have a surface fitting issue in high-dimensional spaces, which is equivalent to the simple neutral network method and has the advantages of a simple computation format, flexible node configuration, small computation amount, and a relatively high precision (Han et al., 2011). In the case where the surface varies gently, this method would deliver a good effect. However, in the case where the surface varies greatly over a short distance and/or the sample values are suspected to have measurement errors or uncertainties, it is not applicable.

4.2.2 Comparison of the interpolation models

In order to determine the influences of the different interpolation methods on the mapping results, we compared the three methods. The parameters required for each method are set in the geostatistical wizard of the geostatistical analysis modulus. In the IDW method, the power was 2, 12 points were involved in the interpolation, the neighborhood was the standard tool, and the size of the output pixel was 4149. In the Kriging method, the Ordinary Kriging method was adopted, the spherical function was used as the semi-variogram function model, the step length was 34 082, and the step number was 12. In the RBF method, the tension spline was selected as the function type, the maximum adjacent element number was 15, the minimum adjacent element number was 10, and the other parameters were set to the defaults.

Figure 4 shows the areal distribution of the thickness points of the WFF. The WFF sections are distributed broadly, but nonuniformly across the Upper Yangtze region. The formation mainly outcrops on the periphery of the Sichuan Basin, such as in northern Guizhou and western Hubei. However, very few cores recovered within the basin reveal its presence. In addition, the thickness of the WFF varies greatly even within the same region. For example, in northeastern Sichuan, the thickness of the WFF shale decreases from 32.3 m in the Chengkou Miaoba section to 7 m in Well WX2 in Wuxi (Liang et al., 2016). Furthermore, the thickness of the WFF shale decreases from 40 m in the Jiuxi Taoyuan section in the Hunan-western Hubei region to 0 m in the nearby Cili Yichongqiao section in Hubei province due to the presence of the Hunan-western Hubei underwater highland.

Based on the thickness data for the WFF (Fig. 4), the thickness contours generated using the IDW, Kriging and RBF methods are shown in Figures 5–7, respectively. As shown in Figure 5, the IDW method provides a precise interpolation with relatively simple principle and parameter settings and delivers interpolation data consistent with the known original data. This method requires that the data points be distributed as uniformly as possible. The interpolation result given by this method is always accurate in areas where abnormal data points vary greatly, such as in the extreme value zone. Therefore, it is applicable in areas where the formation thickness is significantly greater than that of the surrounding areas, such as in the Chengkou Miaoba, Taoyuan Jiuxi, Jishou, and Anhua Yanxi sections. However, the “cow-eye” effect is likely to emerge in the extreme value zone of the data, forming multiple small closed contour circles around the extreme values.

Figure 6 shows two characteristics of the Kriging interpolation. First, since the Ordinary Kriging method is a probability interpolation rather than a precise interpolation, it delivers an interpolation result that is not identical to the thickness of the known original data. Second, since the Kriging algorithm takes into consideration the relationship between the prediction point position and the known position and the spatial correlation of the thickness data, the interpolation result provides a good representation of the overall variation trend in the data, e.g., predicting the merged regional distribution of the thick-bedded formation in Central Hunan.

Figure 7 shows two characteristics of the RBF interpolation. First, since the RBF is a precise interpolation, it allows for the generation of a smooth surface with many data points, instead of predictions between the known maximum and minimum values. Thus, the interpolation yields a smooth and continuous contour. Second, since the thickness varies greatly and the data points are distributed uniformly in central and northern Hunan, there is no merged distribution of the formation in areas that lack the sections, and the method is of limited usefulness in areas where the data varies significantly locally.

A comparison of these interpolation maps reveals that the interpolation results are closely related to the spacing and variation of the data points. The interpolation results are slightly different in areas with closely spaced data points, such as in southwestern Sichuan, northeastern Yunnan, northern Guizhou, and western Hubei, while they are significantly different in areas with widely spaced data points, such as in central Guizhou, central Sichuan, and central Chongqing. In addition, the interpolation results are significantly different in areas where the data points vary greatly and frequently, including in western Hunan-Hubei, while they are slightly different in areas where the data points vary slightly, such as in southern Sichuan. In general, the IDW method is suitable for use with continuously uniformly distributed data, the Kriging method provides a better reflection of the overall variation trend of the data, and the RBF method is less effective in areas where the data varies greatly. All these methods yield good interpolation results for uniformly distributed data and slightly variable values.

To more accurately reflect the lateral variations of these interpolation results, the related five well lateral section (wells W202, N203, JY1, WX2, and J101) (Fig. 8) and a cross-well correlation diagram (Fig. 9) were constructed. The location of the section is shown in Figures 5–7. Since the IDW method is the precise interpolation method, the blue thickness contour always passes through the original data (yellow circles in Fig. 9). The contours generated using the Kriging and RBF methods seldom pass through the original data points. The Kriging method is the probability interpolation method, and hence, it delivers significantly different interpolation results than those of the other methods in areas that lack the control of known points. In general, these 3 methods yield significantly different interpolation results in areas with widely spaced data points, such as the first half of the region (700–830 km) between well WX2 and well J101 in Figure 9, especially for the Kriging method, which delivers significantly different interpolation results. In the second half of the region (830–1000 km), which has uniformly distributed, closely spaced data points, there is no significant difference between the different interpolation results.

Fig. 4 Position of the thickness points of the Wufeng Formation sections in the Middle-Upper Yangtze region.

Fig. 5 2D interpolation map of the Wufeng Formation in the Middle-Upper Yangtze region produced using the IDW method.

Fig. 6 2D interpolation map of the Wufeng Formation in the Middle-Upper Yangtze region produced using the Kriging method.

Fig. 7 2D map interpolation of the Wufeng Formation in the Middle-Upper Yangtze region produced using the RBF method.

Fig. 8 Cross-well correlation of the Wufeng Formation in the Upper Yangtze region through wells W202-N203-JY1-WX2-J101.

Fig. 9 Correlation of the lateral interface of the Wufeng Formation in the Upper Yangtze region for the 3 interpolation methods.

4.2.3 Selection of the optimal model

To determine the optimal model, it was necessary to verify the precision of the interpolation results given by the three methods. Precision inspection is usually accomplished via cross verification or general verification. In this paper, the general verification method is utilized. The training point set and the test point set are selected using the sub-set element tool of the geostatistical tool. In addition, 75% of the section data serves as the training point set for the interpolation, while the remaining 25% of the section data serves as the test point set, which is not be involved in the interpolation computation. Then, the training point set is used for the interpolation. The interpolation result of the training point set is compared with the test point set using the value to point tool of the spatial analysis tool. The difference between the known data and the predicted data is shown statistically, i.e., as the standard deviation (SD) and variance (Var), which provide a measurement of the data fluctuation. A relatively small SD and Var suggest that the data is close to the average and is relatively precise, indicating a good interpolation effect. After the general verification of the three methods, the difference between the predicted values and the original values was determined using the statistical numerical parameter indices. Based on the minimum SD and Var (Tab. 2), the IDW method yields the interpolation result closest to the original values and has the highest precision. Therefore, the optimal model can be derived using the IDW method.

4.3 Spatial law analysis

Based on the selection of optimal interpolation model, a three-dimensional thickness model is quantitatively reconstructed (Fig. 10). This is the first three-dimensional thickness model reconstructed by the quantitative method in the study area, where relevant parameters of sedimentary black shale can be analyzed.

A 3D thickness map of the WFF was derived from spatial interpolation using the IDW method and the Arcscene platform. The WFF shale is distributed broadly and continuously across the Upper Yangtze region, and in general, it thins from southeast to northwest. The shale’s thickness ranges from 0 to 40 m, but is mostly less than 10 m, which reveals a relatively flat topography with topographic changes in localized areas. This map shows that the WFF shale is relatively thick bedded in three areas, including northeastern Chongqing, central Hubei, and northwestern Hunan. These areas surround the western Hunan-Hubei underwater highland on the border of Hunan and western Hubei. Within this highland, the WFF is relatively thin-bedded and is even absent in some areas due to the presence of the Yichang uplift. The formation’s thickness gradually increases with distance from the highland. The thickness of the WFF exceeds 5 m on the periphery of the western Hunan-Hubei highland with a local maximum on the western margin of the Xufengshan uplift in northwestern Hunan, e.g., in Jishou and Taoyuan. This is due to the presence of the southern China depocenter to the south and the transition from the WFF on the platform to the Tianmashan flysch sediments in the depression.

Based on the thickness contours, the WFF black shale was estimated to have a 2D distribution area of 0.67 × 10⁶ km², a total volume of about 3902 km³, and an interpolated average thickness of 6.2 m. This shale was deposited over a 2.83 Ma period. Therefore, the volume of the shale deposited per million years is estimated to be 1379 km³/my and the average thickness of the shale deposited per million years is 2.37 m/my.

Fig. 10 3D map showing the thickness in the Upper Yangtze region derived from the IDW method.

5 Paleo-topographic analysis

In southern China, the Middle Hirnantian strata share a similar paleo-topographic framework with the underlying Late Katian strata (Chen, 2013). Thus, it is possible to infer the paleo-topography of the Late Katian WFF from the paleo-topography of the Hirnantian. Considering that the bentonic organism content, e.g., brachiopods, is an important indicator of the paleo-water depth, we selected the minimum Hirnantian paleo-water depth as the paleo-topographic value for the Wufeng period. In order to study the paleo-topographic distribution characteristics of the WFF in the Upper Yangtze region, a data set containing 275 related sections (including data from 18 wells) was compiled. This dataset included the following attributes: section name, longitude and latitude, lithology, and thickness. Based on the program used to determine the paleo-water depth listed in Table 3 (Zhang et al., 2016), we predicted the paleo-water depth of each section and picked its minimum as the paleo-topographic value to be determined.

After acquiring the paleo-water depth of all of the sections, the Ordinary Kriging method of the geostatistical analysis model was used to create 2D and 3D paleo-water depth maps of the Upper Yangtze region using the following parameters: the Ordinary Kriging as the Kriging type, the semi-variogram as the variable model, the spherical function as the model, 50 maximum adjacent elements, and other default parameters.

Based on the restored 2D paleo-topographic map (Fig. 11) and the Arcscene-based 3D paleo-topographic map (Fig. 12), we determined that the Upper Yangtze region was generally a semi-closed, neritic epicontinental sea environment with relatively gentle topography bordered to the west, south, and north by ancient high lands. The platform was bordered to the west by the central Sichuan ancient land, to the south by the central Guizhou ancient land, and to the north by the Hannan ancient land, exhibiting a paleo-geographic framework characterized by one uplift and four depressions, such as the southeastern Sichuan depression, the central Hunan depression, the southeastern Hubei depression, the northern Hubei depression, and the western Hunan-Hubei underwater highland (Chen et al., 2001).

Due to the presence of the Hunan-western Hubei underwater highland and the adjacent ancient lands, the water body was blocked, retained, and poorly connected to the open sea. This closed, restricted marine environment was one of the major factors controlling the enrichment of the organic-rich shale of the WFF. This almost completely closed reducing water environment is globally unique. The formation of this paleo-geographic framework was the result of the intracontinental orogenic movement in the Yangtze and South China blocks (Zhang et al., 2013). The convergence and compression of these blocks resulted in the formation of local uplifts. For example, the Yichang uplift (Chen et al., 2001) led to the formation of the western Hunan-Hubei underwater highland and directly influenced the development of the graptolite zone and the favorable shale bed.

Well LD1, which was drilled by the China Huadian Corporation on the margin of the western Hunan-Hubei underwater highland, encountered a thin bed of graptolite shale (WF2-LM5 is 20.4 m thick) with unsatisfactory results. However, several wells have been reported to have good amounts of shale gas, including wells N203, W202, and JY1 in the southeastern Sichuan sag and well J101 in the southeastern Hubei depression. Therefore, it is recommended that shale gas wells be placed in areas within the depression and far from the uplift zone.

Fig. 11 Restored 2D paleo-topography of the Middle Hirnantian-Late Katian in the Upper Yangtze region.

Fig. 12 Restored 3D paleo-topography of the Middle Hirnantian-Late Katian in the Upper Yangtze region.

6 Conclusions

Based on the above analysis and discussion, we reached the following conclusions:

1. Compared with traditional artificial interpolation mapping, high-precision thickness mapping based on GBDB and GIS eliminates the disadvantages of the traditional artificial interpolation method through spatial interpolation. With its strong objectivity, large amounts of data, high precision mapping, and high degree of quantification, this method allows for easy subsequent quantitative analysis and the possibility of exploring the spatial distribution mode and law.

2. In the Upper Yangtze region, the WFF and other synchronous sediments contain 8 litho-stratigraphic units and are present in 7 provinces. The black graptolite shale is distributed broadly from east to west across the Yangtze block, while the other synchronous formations with different lithofacies are present on the margins of the platform. These formations were deposited in a normal neritic epicontinental sea environment, with the exception of the flysch sediments in southern Hunan. The sections that do not contain these formations allow for the precise determination of the marine-continent boundary at the time.

3. There is no absolute optimal interpolation method that is generally applicable to the spatial interpolation of geologic data. Therefore, it is necessary to select the ideal interpolation method based on the distribution of the original data and the specific conditions of the study area. We compared the interpolation results of the IDW, Kriging, and RBF methods. Based on general verification, the IDW method is the optimal interpolation method and has the highest precision.

4. Using the IDW method, the WFF in the Upper Yangtze region was estimated to have an average thickness of 6.2 m, an area of 0.67 × 10⁶ km², and a volume of 3902 km³. This shale was deposited over a 2.83 Ma period. Therefore, the volume of shale deposited per million years is estimated to be 1379 km³/my, and the average thickness of the shale deposited per million years is 2.37 m/my.

5. During the Late Katian Wufeng period, the Yangtze platform was bordered by ancient lands to the west, south, and north, exhibiting a paleo-geographic framework characterized by one uplift and four depressions. This setting blocked water circulation, causing the water to be retained and form a closed restricted marine environment, which was one of the major factors controlling the deposition of the organic-rich WFF shale.

6. With the advent of the big data era of geology, the application of GIS-based technology can be applied to almost any resource play with spatial distribution patterns in the world. Based on the geoscience and feature attribute data, the results could be obtained efficiently and rapidly with a quantitative method, presenting broad application prospects.

Funding

This research is funded by China National Science and Technology Major Project (Project No. 2017ZX05035-001).

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

We thank LetPub (www.letpub.com) for its linguistic assistance during the preparation of this manuscript.

References

All Tables

All Figures

Fig. 1 Division and correlation of the Ordovician strata in the Upper Yangtze region (Zhang et al., 2014a). (A) the broad area in the Upper Yangtze region containing the typical black shale of the Wufeng Formation; (B) the area extending from northern Sichuan to southern Shanxi; (C) the area extending from southeastern Sichuan to northeastern Yunnan; (D) the area extending from southwestern Hunan to northeastern Guangxi.

In the text

	Fig. 2 Geographic distribution of the section points of the Wufeng Formation shale and the synchronous strata in the Upper Yangtze region (The yellow circles represent the Wufeng Formation and the synchronous strata points; the orange crosses represent the sections that do not contain these formations; the brown line outlines the Sichuan Basin; and the red line outlines the Yangtze region).
In the text

	Fig. 3 Geographic distribution of the Wufeng Formation and synchronous formations with different lithofacies in the Middle-Upper Yangtze region.
In the text

	Fig. 4 Position of the thickness points of the Wufeng Formation sections in the Middle-Upper Yangtze region.
In the text

	Fig. 5 2D interpolation map of the Wufeng Formation in the Middle-Upper Yangtze region produced using the IDW method.
In the text

	Fig. 6 2D interpolation map of the Wufeng Formation in the Middle-Upper Yangtze region produced using the Kriging method.
In the text

	Fig. 7 2D map interpolation of the Wufeng Formation in the Middle-Upper Yangtze region produced using the RBF method.
In the text

	Fig. 8 Cross-well correlation of the Wufeng Formation in the Upper Yangtze region through wells W202-N203-JY1-WX2-J101.
In the text

	Fig. 9 Correlation of the lateral interface of the Wufeng Formation in the Upper Yangtze region for the 3 interpolation methods.
In the text

	Fig. 10 3D map showing the thickness in the Upper Yangtze region derived from the IDW method.
In the text

	Fig. 11 Restored 2D paleo-topography of the Middle Hirnantian-Late Katian in the Upper Yangtze region.
In the text

	Fig. 12 Restored 3D paleo-topography of the Middle Hirnantian-Late Katian in the Upper Yangtze region.
In the text