A Workflow for Integrated Geological Modeling for Shale Gas Reservoirs: A Case Study of the Fuling Shale Gas Reservoir in the Sichuan Basin, China

Shale gas reservoir evaluation and production optimization both require geological models. However, currently, shale gas modeling remains relatively conventional and does not reflect the unique characteristics of shale gas reservoirs. Based on a case study of the Fuling shale gas reservoir in China, an integrated geological modeling workflow for shale gas reservoirs is proposed to facilitate its popularization and application and well improved quality and comparability. This workflow involves four types of models: a structure-stratigraphic model, reservoir (matrix) parameter model, natural fracture (NF) model, and hydraulic fracture (HF) model. The modeling strategies used for the four types of models vary due to the uniqueness of shale gas reservoirs. A horizontal-well lithofacies sublayer calibration-based method is employed to build the structure-stratigraphic model. The key to building the reservoir parameter model lies in the joint characterization of shale gas “sweet spots.” The NF models are built at various scales using various methods. Based on the NF models, the HF models are built by extended simulation and microseismic inversion. In the entire workflow, various types of models are built in a certain sequence and mutually constrain one another. In addition, the workflow contains and effectively integrates multisource data. Moreover, the workflow involves multiple model integration processes, which is the key to model quality. The selection and optimization of modeling methods, the innovation and development of modeling algorithms, and the evaluation techniques for model uncertainty are areas where breakthroughs may be possible in the geological modeling of shale gas reservoirs. The workflow allows the complex process of geological modeling of shale gas reservoirs to be more systematic. It is of great significance for a dynamic analysis of reservoir development, from individual wells to the entire gas field, and for optimizing both development schemes and production systems.


Introduction
Shale gas reservoirs have low porosity (maximum porosity: 4-5%) and permeability (<1 × 10 −3 μm 2 ). Shale gas accumulation and production depend on two factors. The first is the material base where the shale gas exists, which is mainly controlled by parameters such as organic carbon content and the physical properties of the matrix. The second is the result of hydraulic fracturing of the shale gas reservoir, which is mainly affected by the physical parameters of the rocks and natural fracture (NF) development [1][2][3][4][5]. Determining the underground distribution of the main control parameters for the geological and engineering "sweet spots" of shale gas reservoirs is particularly important and crucial to economical and effective reservoir development [6]. Geological modeling of a shale gas reservoir is a process of visualizing and quantitatively analyzing the constitutive components of the shale gas "sweet spots" from a three-dimensional (3D) spatial perspective [7,8]. A reliable geological model can help us predict the distribution of shale gas reservoir parameters, optimize the placement of wells, effectively increase the reserve utilization rate, and enhance the recovery from the shale gas reservoir. By contrast, determining the spatial distribution of NFs and hydraulic fractures (HFs) in the shale can help guide risk prevention during drilling or fracture utilization and is thus of great significance to shale gas development and production.
Currently, geological modeling techniques are focused mainly on conventional oil and gas reservoirs, such as clastic and carbonate reservoirs. The geological modeling of shale gas reservoirs is still in its initial stage and centers around three types of models: the reservoir matrix, NF, and HF models [9][10][11][12][13][14][15]. The relevant modeling methods still rely on modeling approaches and techniques for conventional oil and gas reservoirs.
However, shale gas formations are both hydrocarbon source rocks and reservoirs and are unique in a number of ways. Geological modeling of shale gas reservoirs is challenging and differs significantly from the modeling of conventional oil and gas reservoirs. For example, shale gas reservoirs are mostly developed with horizontal wells, and their geological modeling relies more on the data acquired from horizontal wells. In contrast to conventional property models, the geological modeling of shale gas reservoirs focuses on depicting the distribution of parameters such as TOC content, brittle mineral content, stress field, and fracability. NFs in shale are of various origins and scales and can affect the results of later hydraulic fracturing. NF models not only reflect geological characteristics but are also required for gas reservoir development. The distribution of HFs is affected by NFs. When simulating HFs, particular consideration should be given to their relationship to NFs. Geological models for shale gas reservoirs contain enormous amounts of information. Their updating and upscaling require considering the mutual constraints between various types of models.
Overall, most studies of the geological modeling of shale gas reservoirs have established geological models for the matrix. Through the classification and description of shale lithofacies, mineral content analysis, and classification of pore space types, distribution models for shale gas reservoir matrix parameters, such as shale lithofacies, total organic carbon (TOC) content, mineral composition, physical parameters, and gas content models, have been built using geostatistical methods [16][17][18]. There are two main approaches for simulating fractures in shale gas reservoirs. One approach detects fractures by seismic attribute analysis and then builds a deterministic model for the fractures. The deterministic model can be decomposed into component models of fractures used varying scale and orientation and mode [19,20]. Another approach statistically analyzes the fracture parameters (e.g., azimuth and dip angle) and then randomly simulates a fracture distribution and builds a discrete fracture network (DFN) model [21,22]. Researchers are more enthusiastic about studying the distribution and prediction of HFs and conducting an extended simulation of the trend of fractures in shale gas reservoirs after hydraulic fracturing using finite element and boundary element methods [23][24][25][26]. However, the existing fracture models rarely reflect the multiorigin and multiscale characteristics of NFs in shale and, more importantly, fail to reflect the relationship between NFs and HFs.
Geological modeling of shale gas reservoirs is complex and unique, and it currently lacks a satisfactory technical workflow. Consequently, researchers build models in their own way, and one fractured reservoir can have many models, whose quality and comparability cannot be ensured. As a result, most research on the numerical simulation of shale gas reservoirs investigates the relevant mechanisms based on conceptual models, and the derived results do not support the efficient exploitation of shale gas. In this study, based on the practical development of the Fuling shale gas reservoir in the Sichuan Basin, China, the main issues associated with the geological modeling of shale gas reservoirs, including data preparation, the scope of work, and relevant solutions, are discussed. In addition, based on the characteristics of shale gas reservoirs, some key issues in geological modeling are identified. The main objective of this study is to provide an integrated workflow for the geological modeling of shale gas reservoirs that is systematic and of reference value to develop more reliable geological models and reduce uncertainty.

Structural and Formation
Characteristics. The Fuling shale gas field is located in the eastern Sichuan Basin, China, and is under the jurisdiction of the Fuling District of Chongqing Municipality (Figure 1(a)). As part of the Sichuan Basin, this area has undergone superimposition and transformation as a result of multistage tectonic movements [27]. The main structural body of this area has been weakly deformed, as reflected by a box-shaped fault anticline morphology that is controlled by a northeast-trending fault and a nearly southnorth-trending fault (Figure 1(b)); that is, the main structural body of this area has a wide and gentle top, where the formations have shallow dip angles and fewer faults, and two steep sides, where faults have developed (Figure 1(c)).
The Fuling shale gas field developed in the Upper Ordovician Wufeng Formation and the Lower Silurian Longmaxi Formation, which are buried at depths of 2,240-4,200 m. In the main area, the burial depth is less than 3,000 m in the north, and it increases gradually southward. The Wufeng Formation is composed mainly of siliceous shale intercalated with carbonaceous shale (generally 5-7 m thick). The Longmaxi Formation exhibits distinct trichotomous characteristics: the lower segment (generally 89 m thick) is composed of dark carbonaceous and siliceous mud shale; the middle segment is composed of turbidite sandstone; and the upper segment is composed of silty mudstone. The Wufeng Formation and the lower part of the Longmaxi Formation were formed from deep-water shelf deposition [28][29][30], maintained an anaerobic environment for a relatively long time, and are the main production formations in the Fuling shale gas field (Figure 2).

Reservoir Characteristics.
Drilling data show that the Fuling shale formations contain a TOC content of 0.17-6.79%. Longitudinally, the lower shale section is the most favorable, and the organic matter (OM) abundance generally decreases from bottom to top ( Figure 2). The brittle mineral components of the shale in this area include quartz, potassium feldspar, and carbonatite. Core measurements show that the 2 Geofluids shale in this area has an average brittle mineral content of 58.1%, is highly brittle, and is prone to fracturing. Organic, clay, and intra-and interfragment pores as well as fractures or structural fractures (including shale bedding fractures (SBFs)) constitute the main shale gas reservoir space in this area [31,32]. The TOC content, brittle mineral content, and porosity of the shale are all spatially heterogeneous, which affect the distribution of the initial gas content. There is no consensus on the classification of shale lithofacies [33,34]. In this area, based on three key indices, i.e., lithological composition, carbonaceous content, and siliceous content, shale formations are classified into nine lithofacies sections of various types [35], which are indicated by sublayers ①-⑨ for ease of modeling ( Figure 2). Of these intervals, due to their carbon-rich, high-siliceous content and high-gas content characteristics, lithofacies sublayers ①-⑤ are the main production intervals in the Fuling shale gas field.
The Jiaoshiba block in the Fuling shale gas field is selected as the study area for geological modeling. The 3D seismic area selected for the study has an area of approximately 85 km 2 . The study area contains dozens of well platforms and a total of 26 wells, including three pilot wells (JY2, JY3, and JY4) and 23 horizontal wells that are spaced at an average distance of approximately 600 m (Figure 1(b)). Three seismic horizons have been interpreted for this area: the bottom and top surfaces of the shale gas reservoir and a surface in the middle of the shale gas reservoir (which, according to our analysis based on time-depth conversion, is the interface between sublayers ⑦ and ⑧). Drilling data include well logs and well test rate data. Well log interpretation can provide various types of information (e.g., TOC and mineral con-tents) for the shale gas reservoir matrix models. Core samples have already been collected from the three pilot wells and subsequently analyzed and tested. The results provide information regarding the microscopic characteristics of the shale gas reservoir. Imaging logs have not been obtained from the study area. Nevertheless, the relevant fracture parameters can be determined based on the imaging logs obtained from the adjacent wells. In addition, microseismic and hydraulic fracturing data obtained from the three wells on well platform A (Figure 1(b)) are available and can be used to simulate and validate HFs.

Content and Workflow of Shale Gas Geological Modeling
The planar distribution of the Fuling shale gas content is mainly controlled by deep and large faults. The shale in the main region contains a higher gas content than that in the peripheral fault area (e.g., east region). However, the gas content varies insignificantly within each region [36][37][38]. In comparison, the shale gas production capacity varies significantly between different locations within each region (particularly the main region), which is mainly related to the burial depth of the structure, the results of fracturing, and the proportion of wells that pass through the high-quality shale sublayers. Therefore, there is a need to build an accurate, finescale structural model for the shale gas reservoir and to simulate the spatial distribution of the shale property parameters within its framework to facilitate the better prediction of shale gas production capacity and well trajectory design.    [36]. The Fuling shale gas field was divided into three regions: eastern, western, and main, according to the structural and development characteristics [37]. (c) Seismic profile across the Fuling shale gas field.

Geofluids
In addition, large-scale multistage fracturing of the long well sections of the horizontal wells is a precondition for the commercial exploitation of a shale gas reservoir, which is, in essence, a process in which a large number of HFs connect to various types of preexisting NFs in the shale gas reservoir. NFs in shale gas reservoirs are of various origins and vary significantly in scale; in addition, they significantly control the generation, propagation direction, and scale of HFs [39]. Therefore, the distribution of fractures in a shale gas reservoir is an extremely important factor that affects shale gas production and is the main object of geological modeling.
Shale gas reservoirs are self-generating, self-storage reservoirs and have extremely low permeability, which makes the shale gas flow between different well shafts difficult. Therefore, researchers have concluded that a geological model for only a single well (or a well group) is needed for the geological modeling of a shale gas reservoir. However, there is still a need to build regional geological models for the following reasons. (1) Having more drilling data can facilitate the more effective use of geostatistics in an analysis (e.g., variograms).
(2) Regional geological models can facilitate a better analysis of the spatial distribution pattern of shale gas reservoir parameters. (3) Regional geological models can also facilitate the prediction of shale gas reserves and yields.
The general approach and workflow used in this study are to integrate various types of data and build four types of geological models: a shale structure-stratigraphic model, shale gas reservoir parameter distribution models, NF models, and HF models (the former two are regional geological models, and the latter two are mainly geological models for single wells (or well groups)) ( Figure 3). In this large workflow, a specific workflow is designed for each type of geological model. Different types of models are connected to one another, and the former models constrain the latter models. Ultimately, these component models are combined to form a complete geological model for the shale gas reservoir. A structurestratigraphic model reflects the basic spatial framework of a shale gas reservoir and is essential for the subsequent modeling of shale gas reservoir parameters. The spatial distribution characteristics of shale gas reservoir properties that are simulated based on a fine, reliable structure-stratigraphic model are more in line with the actual situation.

Geological Modeling for Shale Gas Reservoirs
The development of the Fuling shale gas field focuses mainly on the most favorable shale intervals (lithofacies

Geofluids
sublayers ①-⑤), and most of the horizontal wells pass through sublayers ①-③. Particular attention should be paid to the relationships between the shale intervals and the well trajectories when building the structure-stratigraphic model so that the model can be as consistent with existing shale gas wells as possible. Only when the model is built in this way can it be used to make predictions. In this study, the structure-stratigraphic model for the shale gas reservoir is built using a horizontal-well lithofacies sublayer calibrationbased method ( Figure 4). The top horizon data of a shale gas reservoir can be obtained from all the wells (the pilot wells and the vertical

Geofluids
sections of the horizontal wells). By making full use of the individual interval data combined with the seismically interpreted horizons, models for the key horizons of the shale gas reservoir can be built, and the overall framework of the shale gas reservoir can be determined. In this large framework, a fine simulation of the distribution of each interval is the key in shale structure-stratigraphic modeling. The Fuling shale gas reservoir sublayers are related to the lithofacies. The geological and well log characteristics of each lithofacies are different. Therefore, based on their cores, cutting logs, and well log curves, the lithofacies in the pilot wells can be classified, and comparable lithofacies sublayers and their characteristic features can be determined. It is difficult to observe the log characteristics of lithofacies in the horizontal wells and to compare the classifications of the lithofacies in the horizontal and pilot wells. However, it is possible to first vertically project the trajectories of the horizontal sections of the wells and the well log curves onto a vertical plane and then classify the lithofacies sublayers [35]. Ultimately, the lithofacies sublayers which the trajectories of the horizontal sections of the wells pass through and their variation can be determined ( Figure 5).
By calibrating the lithofacies in each horizontal well, the corresponding interval data can be obtained. One horizontal well can occasionally provide data for the interfaces between multiple sublayers and for the locations of multiple interfacial points of a certain sublayer, based on which the respective data points of multiple sublayers inside the shale can be located ( Figure 5). In addition, based on the calibration results for the lithofacies in the pilot wells and the vertical sections of the horizontal wells, the thickness proportion of each interval in the shale formation is determined. On this basis, a horizon model is built for each sublayer in the shale gas reservoir under the constraint of seismically interpreted horizons.
This structure-stratigraphic model can basically reflect how the wells pass through the shale intervals. However, in a high-precision geological model, an interval is generally longitudinally divided into multiple fine grid layers. If the spatial location of a certain section of a horizontal well inside a certain interval is indeterminate, it will affect the accuracy in the discretization of its well log data into the grid. When this occurs, the distance between the trajectory of the horizontal well and the structural plane of the sublayer encountered during drilling (mainly the top surface of the sublayer) can be determined based on the induction log data on the interval interface to further correct the model ( Figure 6).
The stratigraphic model established in this study contains 9 sublayers, of which the average thickness is about 8 m in sublayer ①, 1 m in sublayer ②, 15 m in sublayer ③, 10 m in sublayers ④-⑧, and 15 m in sublayer ⑨. The established stratigraphic model is divided into vertical grids. Among them, sublayers ①-⑤ are high-quality shale formations, and the vertical grid accuracy of sublayers is 0.5 m, and the grid accuracy of sublayers ⑥-⑨ is 1 m. The total number of model grids in the study area is 1,502,451. The established formation model is cut along the well trajectory to observe the travel relationship between drilling and shale sublayers ( Figure 7). The results show that the formation in the model conforms to the geological understanding and division of sublayers in drilling, and the established stratigraphic framework model is more accurate. On the basis of this model, further modeling of reservoir matrix properties and fractures can be carried out.

Shale Property Parameter Models.
The property parameters of a shale gas reservoir are modeled to analyze and simulate the relevant parameters of the shale matrix, to determine the spatial heterogeneity of each property parameter of the shale gas reservoir, and to predict the distribution of shale gas "sweet spots." The main constitutive parameters for a shale gas "sweet spot" include lithofacies, brittle mineral content, TOC content, physical parameters, and rock mechanical parameters. Most of these parameters can be determined based on well log curves or through calculations using relevant equations.
It is the view of previous researchers that modeling of shale gas reservoir parameters is similar to modeling of conventional oil and gas reservoirs. However, there is a certain correlation between the parameters of a shale gas reservoir (   6 Geofluids Therefore, the lithofacies model should be established on the basis of the TOC model and mineral content models. (2) Both porosity and TOC are positively correlated with gas content (Figure 9(a)). Therefore, the abundance of organic matter is positively correlated with the content of brittle minerals. (3) Organic matter can form pores in the evolution process, and the porosity in shale reservoirs has a positive correlation with TOC to some extent because it includes    (Figure 9(b)). (5) Siliceous minerals contain quartz, and the content of brittle minerals is related to petrophysical parameters. Based on the correlation between property parameters, the hierarchical constraint sequence of all kinds of continuous properties under the control of discrete lithofacies can be clarified, and the source-reservoir coupling constraint framework of shale gas can be established ( Figure 8).
Statistical data analysis is essential prior to model building. Relatively good data analysis can facilitate an understanding of the heterogeneity and distribution trend of the reservoir parameters. The property parameters of the Fuling shale gas reservoir exhibit significant longitudinal variability.
For example, statistical analysis shows the following longitu-dinal distribution characteristics of the TOC content: sublayers ①-⑤ range from 2.6 to 4.5%; intervals ⑥ and ⑦ range from 1.35 to 2.5%; and intervals ⑧ and ⑨ range from 1.2 to 2.0%. Clearly, the TOC content significantly decreases from bottom to top (Figure 10(a)). The average porosity is significantly lower in the eastern region (1.85%) than in the main region (3.25%), where the average porosity is higher in the north (4.5%) than in the south (2.45%) of the plane. Based on the spatial distribution pattern of these parameters, a corresponding trend model can be built that is constrained in the subsequent modeling process (Figure 10(b)).
The reliability of a parameter model is, to a large extent, determined by the modeling method (and/or algorithm). Different modeling methods are applicable to different property parameters. If a property parameter is highly correlated to a certain seismic attribute, then a deterministic method   can be directly used to model this property parameter (e.g., fracture facies). Currently, there are two main types of random geological modeling methods: object-and pixel-based modeling. The sequential Gaussian simulation (SGS) method and the sequential indicator simulation (SIS) method can be used to model lithofacies. The SGS method and the truncated Gaussian simulation method can be used to simulate the brittle mineral content, TOC content, and porosity ( Figure 11). The variation function and trend model of each property are used to control the distribution proportion of the parameter in the three-dimensional space. Meanwhile, the nextlevel property model in the constraint framework is used as the second data for the collaborative constraint, and the property models are established successively (Figures 12(a)-12(c)). TOC is established according to trend model constraints. The brittle mineral is established according to the trend model constraint, and the TOC model is the second data. According to lithofacies classification, the discrete model of lithofacies is established according to the TOC model and brittle mineral model. The porosity is controlled by the trend model and controlled by lithofacies. The compressible index is constrained by the in situ stress trend model, and the brittle mineral content model is the second data.
Finally, by comprehensively analyzing the spatial distribution of the property parameters of the shale gas reservoir in question, the geological pattern of shale gas "sweet spots" is summarized. In addition, the distribution of shale gas "sweet spots" is predicted based on the contribution of various property parameters to the shale gas content and production capacity. The TOC and high-quality lithofacies with high gas content in the study area are mainly carbon-rich calcium-bearing siliceous shale and carbon-rich mixed shale. The statistical correlation between TOC, porosity, and brittle mineral content and gas content in high-quality lithofacies is 8.31, 7.56, and 6.89, respectively. In other words, the contribution of organic matter abundance, physical property, and lithology to shale gas "sweet spots" is about 36.5%, 33.2%, and 30.3%, respectively. The established TOC model, porosity model, and brittle mineral content model are, respectively, multiplied by the corresponding sweet spot contribution ratio, and the distribution model of shale gas sweet spots is obtained by fusion (Figure 12(d)). Through the analysis of the model, it is found that the shale gas "sweet spots" are distributed continuously in the northeast area, while the scale of the northwest area is smaller, and the shale gas "sweet spots" are scattered in the southern tectonic slope, which is consistent with the actual production capacity of shale gas in the study area. The natural fracture (NF) system plays a decisive role in the dissipation and effective development of shale gas and significantly affects the effectiveness of subsequent fracture-induced transformation. Therefore, NF models play an important role in the geological modeling of shale gas reservoirs. The effective integration of the fracture data obtained from multidisciplinary (e.g., geological, geophysical, and geomechanical) analysis is thus required to accurately simulate the spatial distribution of NFs.
NFs in shale have various origins and scales. Most largescale fractures are of structural genesis, whereas most relatively small-scale fractures are horizontal shale bedding fractures (SBFs) [10]. Therefore, when building NF models, it is necessary to analyze the origins of fractures and classify their scales. Currently, the common fracture modeling methods first decompose into fracture models of various scales and orientations [20], use them to mutually constrain one another, and ultimately integrate them.
The scale, quantity, and distribution of fractures vary with different scales. According to the special properties of different reservoirs, the grading and standard of fracture scales are different. The observation and description methods of fractures of different scales are not consistent, which largely controls the corresponding modeling methods and modeling constraints. The different scales of fractures observed in shales include microfractures, joints (or bedding and foliation), and faults, with scales ranging from microns to kilometers. According to the fracture extension length, the fractures developed in the shale can be divided into three scales in this study, e.g., large-scale fractures (hundred meters to kilometers), medium-scale fractures (meters to hundred meters), and small-scale fractures (micrometers to centimeters).  (Figure 13(a)). The in situ stress model should be analyzed firstly based on the strain characteristics indicated by structural deformation and strain variation in 3D space. Then, this in situ stress model is used to predict the properties (e.g., fracture density, azimuth, and dip angle) of NFs based on the geomechanical properties (e.g., Young's modulus). On this basis, a large-scale DFN model is generated (Figure 13(b)).

Medium-Scale Fracture Modeling.
Fractures that are observed on imaging logs (which cannot be easily identified from seismic data) are mostly sub-first-grade fractures. Medium-scale fracture modeling usually takes the seismic prediction fracture attribute as the spatial constraint of fracture development density; the spatial distribution of discrete fractures is simulated by the stochastic method through logging fracture interpretation parameters.
The establishment of the DFN model of these fractures requires fracture density as the simulation generation driver of fracture spatial distribution. Core and outcrop observations reveal that mesoscale fractures are usually associated with large fractures or faults. Therefore, the occurrence of fractures of this scale is similar to that of large-scale fractures, and their density is, to a certain extent, related to the scale of and distance from large-scale fractures (the shorter the distance from large-scale fractures is, the more accompanying fractures there are).
In the Fuling shale gas field, the extent of development of NFs is controlled by lithofacies, strata thickness, and 11 Geofluids structural location [39]. A clear understanding of the scale and orientation of fracture development is of particular importance prior to fracture modeling [19,20]. Aiming at the modeling of medium-scale natural fractures in shale gas, multiple fracture facies can be divided by considering stress factors and combining with fracture distribution, tectonic stress, and shale gas production analysis (Figure 14(a)). Fracture facies can be regarded as a comprehensive representation of fracture properties. Different fracture facies have great differences in the fracture azimuth, dip angle, strength, and scale ( Figure 14(b)). The maximum principal stress direction of the shale in the study area is the NW-SE compression  Figure 13: Large-scale fracture distribution of shale gas reservoirs in the study area. Based on seismic data interpretation and in situ stress analysis, six large-scale faults, namely, F1-F6, were characterized in the study area. Geofluids direction, and it has the characteristics of a synclineanticline alternate landform on the box-type anticline tectonic background (Figure 14(c)). According to the correlation between the development law of faults in the study area and the tectonic stress, combined with the distribution of faults, the study area is divided into five fracture facies ( Figure 14, Table 1). There are differences in shale gas productivity of different fracture facies. For example, the average cumulative production of natural gas in fracture facies 2 is above 8,000 × 10 3 m 3 , while that in fracture facies 3 is below5,000 × 10 3 m 3 . Based on the theory of fracture genesis and geomechanics, the fracture development characteristics of each fracture facies can be defined as follows:

Geofluids
(1) Compressional High-Strength Fracture Facies. The fracture facies is mainly located in the southeast of the study area. Affected by the fault (F1), the fractures are mainly compressional fractures, which are parallel to the F1 in occurrence and distributed in the NE strike, and the fracture intensity is high.
(2) Extensional Medium-Strength Fracture Facies. The fracture facies is mainly located in the southwest of the study area. Due to the strong east-west structural compressional deformation, a low uplift is formed. Some NE-trending extensional fractures are developed in the curved part, and the occurrence of fractures is mostly NE trending. Fractures in this region have large extension length and good continuity, and the fracture intensity is not as good as that of fracture facies 1, but it is also relatively developed.
(3) Compressional Low-Strength Fracture Facies. The fracture facies is developed in the northwest of the middle part of the study area and is generally located in the subsyncline of the main area. This area is a syncline of low amplitude, with no obvious compressional deformation and weak in situ stress. Two clusters of fractures are developed in the EW and NW strike, and the overall fracture intensity is low.
(4) Medium-High-Strength Fracture Facies. The fracture facies is located in the intersection area of F1 and F4 faults in the eastern part of the study district. The fractures in this fracture facies are controlled by both the compressive and shear properties, and the overall occurrence is chaotic. The fractures are mainly developed in two clusters of the NE and NWW strike, and the fracture intensity is relatively strong. Based on the obtained fracture occurrence and fracture strength parameter data of each fracture facies, the parameter input for fracture modeling of each fracture facies was obtained according to the corresponding fracture development law and occurrence (Table 1). Discrete fracture fragments were generated randomly by the conditional simulation method for each fracture facies. The occurrence of fracture fragments was controlled by the input parameters of different fracture facies. The number and distribution of fracture fragments were uniformly controlled by the fracture strength model. A unified and complete discrete fracture network (DFN) model in the study area was constructed by fusing the discrete fracture segment models simulated by the five fracture facies (Figure 15(a)). According to the fracture model of the study area, the number of fractures encountered by the drilling wells was calculated to be 83 at most and 2 at least. The fracture model reflected the heterogeneity of fracture development (Figure 16). Each fracture in the discrete fracture network (DFN) model has relevant fracture properties, including fracture aperture, fracture permeability, and fracture conductivity, in which the number, direction, length, and area of fractures are all known data, and the fracture aperture data can be obtained according to FMI information.
Based on the properties of fractures in different fracture facies, different parameters of fracture aperture are used. For example, the input values of fracture aperture in fracture facies 1 and facies 5, which are dominated by compressional fractures, are relatively small. However, in fracture facies 2, which mainly develops extensional fractures, the fracture aperture needs to be given a larger value so as to reflect that the connectivity coefficient (which affects the fracture permeability) is different due to the property of fractures (Figure 15(b)).

Small-Scale Fracture
Modeling. The shale bedding planes are the most prominent characteristic of a shale gas reservoir. Under proper conditions, shale bedding planes are prone to fracturing, after which they become a type of shale gas-transporting fracture-horizontal SBFs. There are two key areas in the simulation of this type of fracture: the acquisition of SBF parameters and the factors that affect the development of SBFs (or constraints for modeling). Because SBFs are relatively small in scale and are difficult to identify from seismic and imaging log data, core observation and analytical testing are currently the primary means for  14 Geofluids determining their parameters. Through microcomputed tomography (digital core construction) and immersion tests of core samples, the development characteristics and aperture of SBFs in a shale gas reservoir can be determined visually [40,41]. In addition, high-resolution scanning electron microscopy Maps analysis can be employed to determine the frequency, aperture, and fillings of SBFs. By statistically analyzing the relationships between their parameters (e.g., fracture length, aperture, and density) and various geological variables (e.g., the mineral content of the rocks), the development pattern of and the factors that affect horizontal SBFs can be determined. For example, the comprehensive index

Geofluids
(fracture density × fracture aperture) for SBFs in a shale gas reservoir was positively correlated with the siliceous mineral content and was negatively correlated with the clay mineral content [39].
The discrete fracture modeling method can be used to model small-scale fractures, and the DFN model can be established by random simulation by giving lamellation fracture parameter data obtained from tests. The other method is to establish the lamellation fracture model by means of equivalent matrix reinforcement by clarifying the lithofacies types and their quantitative relationship with the matrix reservoir parameters (organic, etc.).

Modeling of the Fracture Network after Hydraulic
Fracturing. The horizontal drilling-induced hydraulic fracture is the key to the efficient development of shale gas. In other words, it is necessary to generate a large-scale, complex fracture network in a shale gas reservoir to increase the contact area between the fractures and the matrix. After volumetric fracturing of the shale gas reservoir, the initial production is not only high but also more conducive to long-term stable production. Fracturing and stimulation will be the most important technical means to develop shale gas. The numerical simulation shows that the larger the volume of shale gas reservoirs is, the better the stimulation effect will be [42]. However, due to the multimedia texture of shale gas reservoirs after fracturing, especially the complex flow relationship between natural fractures (NFs) and hydraulic fractures (HFs), there are multiple flow patterns in the space of shale gas reservoirs [43]. Current conventional gas reservoir models are not necessarily suitable for shale gas reservoirs. Numerical simulation of shale gas reservoirs requires further research on the coupled flow between shale gas reservoirs and fractured horizontal wells so as to clarify the near-wellbore seepage law of fractured horizontal wells. Both HFs and NFs exist in developed shale reservoirs. However, most geological models used in gas reservoir engineering are prefractured models without HFs or simple postfractured conceptual models. Therefore, the simulation of HFs is essential to achieve a more accurate numerical simulation of a shale gas reservoir and a more effective evaluation of the effectiveness of hydraulic fracturing.
The conventional simulation of HFs requires a number of assumptions, e.g., that the shale formation is homogenous or the existence of NFs is not considered and that the HFs are two-winged planar fractures that are distributed symmetrically on the two sides of the well [44,45]. However, both hydraulic fracturing tests and microseismic monitoring data show that the distribution of HFs is neither symmetric nor regular [46]. The distribution of HFs in shale is significantly affected by the intrinsic properties of the shale, including the physical properties of the shale rocks and the distribution of NFs [47][48][49]. Therefore, the key to the accurate simulation of HFs lies in the determination of the relationship between NFs and HFs ( Figure 17).
In this study, two technical methods of forward and inversion modeling are discussed to construct HF modeling.
The forward simulation approach is employed to simulate the propagation and trend of HFs using hydraulic fracturing data based on the NF models. This approach includes three main steps: (1) statistically analyzing the hydraulic fracturing curves and fracture monitoring results, evaluating the parameters (e.g., azimuth and altered volume) of the HFs, and describing an HF distribution model for each area based on the heterogeneous physical properties of the rocks in the Fuling shale gas field; (2) inputting relevant hydraulic fracturing parameters into the NF model for each area to simulate the propagation of the HFs; and (3) analyzing the settling of the proppant in the HFs, evaluating the flow conductivity of the HFs, and building property models for the HFs. The advantages of this approach are that it is fully combined with the NF models and reflects the factors that control the distribution of HFs. However, this approach, to a large extent, passes the uncertainty in the simulation of NFs onto the simulation results, which are therefore overly random.
Another approach for simulating HFs is to employ a certain algorithm to directly simulate HFs based on microseismic data (inverse simulation). Here, well platform A, from which the microseismic data were obtained, is used as an example. The inverse simulation approach includes the following steps: (1) loading the drilling-induced hydraulic fracturing-related microseismic event data, calibrating the locations of the HFs, and building a 3D distribution map for each hydraulic fracturing interval (a "nephogram" of the microseismic event); (2) determining the constraints for the geometric parameters of the main fracture based primarily on the time attribute of the occurrence of the microseismic event points and simulating the possible development process of the fractures (fracture paths); (3) describing the density distribution attribute of the microseismic event point set and the extent of development of HFs within the space where the microseismic event is distributed; (4) calculating the locations where HFs are generated and their distribution directions during the hydraulic fracturing process based on the azimuth attribute and energy level (seismic magnitude or amplitude and signal-to-noise ratio) of the event points; and (5) building a DFN model, calculating the distribution area and swept volume of the HFs, and building HF property models. This approach can effectively integrate microseismic data. The simulation results obtained using this approach are, to a certain extent, consistent with the actual distribution of the HFs. However, microseismic signals can be easily affected by ambient noise, resulting in the generation of some false signals. In addition, the simulation results obtained using this approach fail to reflect the effects of NFs and the properties of the shale gas reservoir on HFs.
This study results from two simulation approaches for HFs. In general, the study of hydraulic fracture modeling is still relatively independent of the shale reservoir matrix model and natural fracture model, and the simulation and morphology analysis of hydraulic fracture propagation do not effectively connect the model research foundation of the shale reservoir matrix property parameter heterogeneity field and natural fracture distribution. The existing hydraulic fracture simulation methods can only simulate the fracture path based on mathematical forward modeling or 16 Geofluids microseismic inversion, and the simulation results cannot guarantee the textural rationality of the spatial distribution of the fracture and the agreement with the actual monitoring position ( Table 2). Because of the complexity of hydraulic fractures and the strong spatial heterogeneity, it is difficult to obtain a comprehensive and accurate understanding of fractures only by using a single method. Therefore, in order to objectively explain the distribution characteristics of the fracture system, it is necessary to fuse multiple types of data to obtain the distribution and characteristic parameters of hydraulic fractures in the reservoir and to predict the spatial distribution of hydraulic fractures by integrating multidisciplinary information.

Relationships between Various Types of Models.
A complete geological model for a shale gas reservoir includes four types of models: a structure-stratigraphic model, reservoir parameter models, NF models, and HF models. These four types of geological models differ in their objectives, functions, and contributions to shale gas development and production ( Figure 18). The structure-stratigraphic model is mainly used to finely describe the morphology and internal sequential structure of the shale gas reservoir. The property parameter models focus on reflecting the spatial heterogeneity of the properties (e.g., lithofacies, mineral content, TOC content, and rock mechanical properties) of the shale gas reservoir. The NF models show the distribution characteristics of NFs   The simulation results cannot guarantee that the hydraulic fracture distribution is consistent with the actual fracture location The simulation results do not reflect the influence of natural fractures and shale reservoir properties on hydraulic fractures The simulation results cannot guarantee that the distribution of hydraulic fracture conforms to the theory of elongation 17 Geofluids of various types and scales in 3D space. The HF models focus more on analyzing the factors that affect the distribution of HFs and evaluating the results of hydraulic fracturing.
Various types of geological models for a shale gas reservoir are closely related to one another. Attention should be paid to their sequences when building these models ( Figure 19). The structure-stratigraphic model is the first model that must be built and provides the basis for building other types of models. The quality of the structure-stratigraphic model is determined by the accuracy of hard data (e.g., well log and core data) entry into the grid, which affects the reliability of the subsequent shale parameter and fracture models. Based on the structurestratigraphic model, property parameter distribution models are the second models that must be built for the shale gas reservoir. Of these models, the brittle mineral content and rock physical parameter models are important foundations for analyzing the geomechanical characteristics of the shale and can provide constraints for NF development density and fracabil-ity for HFs. The first two types of models are mainly used to build the framework of the shale gas reservoir matrix, in which NFs and HFs exist as another type of medium. The distribution, scale, occurrence, and spatial arrangement of NFs will affect the density, trend, and aperture of HFs. Therefore, it is necessary to simulate HFs based on the NF models. If the hydraulic fracturing data are still difficult to match after performing multiple simulations, then the NF models must be corrected. Each type of model represents a stage of geological modeling of the shale gas reservoir and reflects the data and geological understanding at the time of model establishment.
In the later part of the modeling process, more geological and engineering data may become available; consequently, the models must be updated. When adjusting and updating the models, attention should be paid to the relationships between the models, and efforts should be made to avoid correcting the models independently.    Figure 19: The key links in the establishment of the geological model of shale gas reservoirs. Each link involves the fusion of models, and the method and quality of fusion control the reliability of the final geological model. 18 Geofluids Through the geological modeling of a shale gas reservoir, we can provide oil reservoir engineers with a complete data volume that contains various types of information. In the entire workflow, there are multiple model integration processes, which are the key nodes of model building. Of the parameter models used for the shale gas reservoir, the lithofacies, TOC content, mineral content, and rock physics models are integrated to form a "sweet spot" distribution model for the shale gas reservoir; the large-, medium-, and small-scale NF models are integrated to form a multiscale NF distribution model for the shale; the NF and HF models are integrated to construct a fracture model for the shale; and the matrix parameter models and the fracture models are ultimately integrated to form a complete, comprehensive geological model for the shale gas reservoir (Figure 18). The integration methods used for different types of models vary. The matrix parameter models are integrated through certain calculations. The fracture models are integrated using relevant algorithms. The integration of the fracture and matrix models requires grid transformation.
Model integration determines whether the information in the models can be accurately passed on, and the information quality not only determines the accuracy of the final geological model for the shale gas reservoir but also affects the reliability of the numerical simulation.

Problems and Future
Work. This study presents a systematic geological modeling workflow that is specifically applicable to shale gas reservoirs. In the practical implementation process, this workflow faces a number of problems and has many areas that require improvement.
Geological modeling of a shale gas reservoir is a complex process that involves a number of modeling methods. These modeling methods each have their respective limitations. For example, although NF simulation reflects the multiscale characteristics of NFs, it fails to display the actual spatial arrangement of NFs of various scales. In future work, more attention should be paid to the selection and innovation of modeling or simulation algorithms.
Robust and reliable geological models, particularly NF and HF models, are not possible without the continuous optimization of modeling software. Various software platforms have their respective modeling advantages and/or capacities. The geological modeling of a shale gas reservoir requires simulating both the reservoir matrix and the fractures. Therefore, building a complete geological model for a shale gas reservoir requires the coordination of (or transfer between) multiple software platforms. Developing compatible modeling algorithms and/or integrated software platforms is an area that warrants future investigation.
As mentioned previously, different types of geological models for a shale gas reservoir mutually constrain one another, resulting in the transfer of model uncertainty. In fact, it is essential to attach importance to the evaluation of model uncertainty in a geological modeling workflow for shale gas reservoirs due to the following benefit. By analyzing the uncertainty factors for each modeling step and determining the weight of influence of each factor, the model quality can be more specifically improved during the model updating, integration, and upscaling processes.

Conclusions
This study proposed an integrated geological modeling workflow for shale gas reservoirs to facilitate its popularization and application and well improved quality and comparability. The main issues associated with the geological modeling of shale gas reservoirs are discussed: (1) For structure-stratigraphic modeling of shale gas reservoirs with a large number of horizontal wells, this study provided a method to build a stratigraphic model that makes full use of the lithofacies sublayers of horizontal wells in shale gas reservoirs. On the basis of the lithofacies sublayer calibration of horizontal wells, the changes of lithofacies encountered by horizontal well drilling can be determined. More lithofacies sublayer interface information of horizontal wells can be obtained and used to correct the sublayer structural plane. The workflow ensures that the path of drilling in shale gas reservoirs is consistent with the actual work area and improves the reliability of the stratigraphic model (2) For matrix modeling of shale gas reservoirs with dual roles of source rocks and reservoirs, this study provided a method to model shale gas properties including control parameters of shale gas "sweet spots." Through geostatistical analysis of reservoir parameters, analysis of the spatial heterogeneity distribution trend of reservoir parameters, and superimposition of the source-reservoir coupling relationship among reservoir parameters, a hierarchical constraint framework of reservoir parameters was formed to realize the establishment of the matrix model that can describe the spatial distribution of shale "sweet spots" (3) For natural fracture modeling of shale gas reservoirs with multiscale fracture systems, this study provided a method to establish the natural fracture model of shale gas reservoirs based on fracture facies characterization. Based on the analysis of geological structure and the prediction of geophysical in situ stress, the fracture facies was divided. Combined with the relationship between geomechanics and fracture genesis, the fracture development characteristics of each fracture facies were determined and the fracture distribution parameters can be obtained as the condition data of fracture modeling. The DFN model of fractures was established by using the idea of faciescontrolled modeling (4) There are two main methods for the modeling of hydraulic fractures in shale gas at present. One is forward modeling of the initiation and propagation trend of hydraulic fractures. The other is to construct a discrete network model of fractures by using microseismic data. However, it is still difficult for the existing hydraulic fracture modeling methods to simulate both the distribution morphology of the hydraulic fractures and the actual monitoring location. It is one of the important development trends in the field of shale gas exploration and development to establish the postfracturing geological model including natural and hydraulic fracture networks (5) Building a robust geological model for a shale gas reservoir involves the building of four types of models: a structure-stratigraphic model, reservoir (matrix) parameter distribution models, NF models, and HF models. The methods used to build different types of geological models vary. Different types of geological models mutually constrain one another. The entire geological modeling workflow for shale gas reservoirs involves the processing and integration of multisource and multiscale data as well as multiple model integration processes

Data Availability
All references to the article are available. Some new ideas for geological modeling of shale gas reservoirs introduced in this paper can be shared. The 3D model data in the study area (i.e., Figures 5,7,11,12,and 15) are not available because of the commercial secret rules by SINOPEC corporation. Some of the raw data for shale gas modeling (i.e., Figures 2,9,10,and 17) are available by contacting the corresponding author (X. Shang, email: shangxf17@163.com).

Conflicts of Interest
The authors declare that they have no conflicts of interest.