Geostatistical Three-Dimensional Modeling of a Tight Gas Reservoir : A Case Study of Block S 6 of the Sulige Gas Field , Ordos Basin , China

In this study, three-dimensional (3-D) geostatistical models were constructed to quantify distributions of sandstone and mudstone. We propose a new method that employs weight coefficients to balance the sandstone and mudstone data from irregular well patterns during stochastic modeling. This new method begins with classifying well groups according to well distribution patterns; areas with similar well distribution patterns are classified within the same zone. Then, the distributions of sandstone and mudstone for each zone are simulated separately using the sequential indicator simulation (SIS) method, and the relevant variogram parameters for each zone are computed. In this paper, we used block S6 of the Sulige Gas Field in Ordos Basin in China as a case study. We evaluated the quality of each set of parameters through the vacuation checking method; certain wells were removed to generate equiprobable realizations using different seed numbers. Subsequently, the variogram parameters for the entire S6 area were obtained by assigning different weight coefficients to the parameters of each zone. Finally, a quality assessment of the sandstone and mudstone models of the S6 area was conducted using the horizontal wells, which were not involved in the stochastic modeling process. The results show that these variogram parameters, which were calculated using weight coefficients, are reliable.


Introduction
It's an obvious worldwide trend that unconventional hydrocarbon accumulations are gaining more and more popularity in the oil and gas industries.Unconventional energy sources such as tight gas or oil have become essential factors that affect the energy consumption of the world [1].It's an associative issue to quantitatively assess unconventional gas or gas systems.As Schmoker [2] pointed out, it's a basic approach for resource-assessment based on wells and reservoir-simulation models.Therefore, it's necessary to establish reliable geologic models for reservoir simulation.As basic types of geologic models, sandstone and mudstone models are important, and investigating sandstone distribution provides significant information for studies related to fluvial sedimentary system in terms of flow units, physical properties, and net reservoirs.Sand body connectivity directly affects well pattern design, field development schemes, and other engineering considerations.Therefore, determining the distributions of sandstone and mudstone is key for petroleum engineering.Three-dimensional (3-D) geological modeling of sandstone and mudstone plays an important role in reservoir prediction.Geologists worldwide have performed reservoir prediction research studies based on well-logging data, dynamic testing data, seismic data, and the modern sedimentary record; modeling methods have included stochastic simulation methods such as modular neural network simulation, SIS (sequential indicator simulation), multiple-point geostatistical methods, object-based facies modeling methods, and other modeling methods.Holden et al. proposed an improved object-based method for the modeling of fluvial reservoirs conditioning on volume ratios, well observations, and seismic data [3].Lynds and Hajek proposed a process-based conceptual model for understanding and predicting the distribution and geometries of mudstone in sandy braided-river deposits that was based on modern braided rivers and ancient braided-river deposits [4].Deustch designed a sequential indicator simulation program named BlockSIS to gain geologically and statistically homogeneous subdivisions for soft secondary data such as porosity and permeability [5].Pyrcz et al. designed an event-based and conditioning-flexibility program, ALLUVSIM, especially for fluvial depositional settings [6].So far, more and more stochastic modeling methods have been created, and they have been applied to many case studies [7][8][9][10][11][12].For the Sulige gas field, the largest onshore tight gas field, numerous reservoir geological modeling studies using different methods have been published.Jia et al. established a geological model of the sandy braided-river system for an infilling zone in Sulige gas field based on the study of outcrops in the city of Datong of Shanxi Province, and reservoir heterogeneity was analyzed [13].Chen et al. [14,15], Zhao et al. [16], and Yang et al. [17] built geological models for Sulige gas field according to analyses of electrofacies and cores; and predictions on potential areas were made, which were summarized from numerous stochastic realizations.Han et al. [18], Wu, and Fu [19] applied the multiple-point method to build geological models for the tight gas reservoir in Sulige gas field, and horizontal wells were used to check the models' quality.Guo et al. simulated the distributions of rock facies, sedimentary microfacies, and physical properties of the braided-river system of Sulige gas field using a particular workflow based on logging data and seismic data.Sandstone distribution models were validated by vacuation method and dynamic data [20].However, among these studies, research specifically on sandstone and mudstone distribution is rare.
As the case studies mentioned above, the well density of the studied area was not taken into account when simulating distributions of sandstone and mudstone.When SIS is used for geological unit modeling, the distances between sampling data in the major direction (generally refer to the source direction) do affect the fitting process of the experimental variogram function.Meanwhile, within a fluvial-dominated area, there may be several sub-directions due to frequent flow changes, and it's much more reasonable to take these sub-directions into consideration during geological modeling.
In this paper, we propose a weight-based method based on the SIS method, to describe and evaluate the sandstone and mudstone distributions of block S6 of the Sulige gas field, which will take the existing information on the irregular well patterns of this field into account and integrate directional wells with horizontal wells.Although this method is not that novel in terms of algorithm, the proposed workflow is valuable for post-processing modeling parameters.This is especially true for large-sized fields that are of good structural integrity, or for areas with greater fault development, under the condition that the fault systems formed after the reservoir was shaped.
The Sulige gas field is located northwest of the Yishan Slope in Ordos Basin (Figure 1a), and its main gas-bearing formations are the Shanxi Formation and the Lower Shihezi Formation of Permian age (Figure 1b).The Sulige gas field is a tight gas field, typical to those in China, characterized by low porosity, permeability, and pressure, and high heterogeneity [1].The study area, block S6, lies in the center of Sulige gas field, and exhibits a gentle and west-inclining monoclinal structure (Figure 1a).The target reservoirs are the Shihezi and Shanxi Formations (Figure 1b).The former is divided into two parts; the upper part includes the He1, He2, He3, and H4 members, and the lower part consists of the H5, H6, H7, and H8 members.The H8 member is subdivided into two parts: the upper member (H8S), and the lower member (H8X) (Figure 1b).The Shanxi Formation comprises the Shan1 (S1) member, and the overlying Shan2 (S2) member.

Materials
In the S6 area, a total of 390 directional wells (including the vertical wells) and 48 horizontal wells have been drilled to date (Figure 2a).All of the sandstone and mudstone interpretation data from these directional wells were used as basic data for stochastic modeling (Figure 2b).The corresponding data from the horizontal wells were not used in stochastic modeling; they were used only to verify the quality of the sandstone and mudstone model.The mudstone is subdivided into two types: interlayer A, which refers to the mudstone between two individual adjacent layers, and interlayer B, which refers to the mudstone within an individual layer (Figure 3a).This classification is based primarily on the difference in thickness between these two types of mudstone layers (Figure 3c,d).Within the study area, a total of six zones were subdivided, and their boundaries roughly determined (Figure 2a).The statistical information of these zones is shown in Table 1.According to the well density and weight coefficient data, these zones reflect the considerable heterogeneity of well distribution in the study area.

Materials
In the S6 area, a total of 390 directional wells (including the vertical wells) and 48 horizontal wells have been drilled to date (Figure 2a).All of the sandstone and mudstone interpretation data from these directional wells were used as basic data for stochastic modeling (Figure 2b).The corresponding data from the horizontal wells were not used in stochastic modeling; they were used only to verify the quality of the sandstone and mudstone model.The mudstone is subdivided into two types: interlayer A, which refers to the mudstone between two individual adjacent layers, and interlayer B, which refers to the mudstone within an individual layer (Figure 3a).This classification is based primarily on the difference in thickness between these two types of mudstone layers (Figure 3c,d).Within the study area, a total of six zones were subdivided, and their boundaries roughly determined (Figure 2a).The statistical information of these zones is shown in Table 1.According to the well density and weight coefficient data, these zones reflect the considerable heterogeneity of well distribution in the study area.

Method
A new method based on SIS has been developed to quantitatively process the variogram parameters (major range and minor range) during stochastic modeling, and the horizontal wells were integrated to evaluate the quality of the sandstone and mudstone model (Figure 4).The SIS method

Method
A new method based on SIS has been developed to quantitatively process the variogram parameters (major range and minor range) during stochastic modeling, and the horizontal wells were integrated to evaluate the quality of the sandstone and mudstone model (Figure 4).The SIS method

Method
A new method based on SIS has been developed to quantitatively process the variogram parameters (major range and minor range) during stochastic modeling, and the horizontal wells were integrated to evaluate the quality of the sandstone and mudstone model (Figure 4).The SIS method was elaborately chosen for the simulation of sandstone and mudstone.The distribution or boundary of sandstone or mudstone is not simply equal to the distribution or boundary of sandy or muddy depositional elements, respectively.Moreover, the sandy depositional elements in the study area are multiple, and the boundaries between them are hard to recognize.Therefore, it's reasonable to select SIS when there are no clear genetic shapes that could be put into object-based modeling, and the SIS models are reasonable in settings where there are no large-scale curvilinear features [5].To implement this new approach, the target reservoir was first subdivided into several zones according to the well patterns, where areas with different well densities were classified as different zones (Figure 2a; Table 1).Then, all lithofacies (sandstone and mudstone) interpretation data from the directional wells were used to stochastically predict the sandstone and mudstone distribution for each zone to acquire the corresponding variogram parameter set.Subsequently, based on well density, each variogram parameter set was assigned different weight coefficients to obtain a variogram parameter set for the target area.The weight coefficient is defined as follows: and the well density is calculated as follows: W i -Weight coefficient of each subdivided zone; D-Mean well density of the target area, wells/km 2 ; D i -Mean well density of the subdivided zone, wells/km 2 ; N-Number of wells; S-Area, km 2 .
was elaborately chosen for the simulation of sandstone and mudstone.The distribution or boundary of sandstone or mudstone is not simply equal to the distribution or boundary of sandy or muddy depositional elements, respectively.Moreover, the sandy depositional elements in the study area are multiple, and the boundaries between them are hard to recognize.Therefore, it's reasonable to select SIS when there are no clear genetic shapes that could be put into object-based modeling, and the SIS models are reasonable in settings where there are no large-scale curvilinear features [5].To implement this new approach, the target reservoir was first subdivided into several zones according to the well patterns, where areas with different well densities were classified as different zones (Figure 2a; Table 1).Then, all lithofacies (sandstone and mudstone) interpretation data from the directional wells were used to stochastically predict the sandstone and mudstone distribution for each zone to acquire the corresponding variogram parameter set.Subsequently, based on well density, each variogram parameter set was assigned different weight coefficients to obtain a variogram parameter set for the target area.The weight coefficient is defined as follows: and the well density is calculated as follows: Wi-Weight coefficient of each subdivided zone; D-Mean well density of the target area, wells/km 2 ; Di-Mean well density of the subdivided zone, wells/km 2 ; N-Number of wells; S-Area, km 2 .The final variogram parameter set (PS) for the target area was obtained via the following equation: The final variogram parameter set (PS) for the target area was obtained via the following equation: Thus, the variogram parameter set for the entire study area was computed directly rather than via fitting with a spherical, Gaussian, or exponential variogram model.Finally, to evaluate the quality of the sandstone and mudstone distribution models, named "model (PS)", numerous vertical profiles of the predicted sandstone and mudstone models through the horizontal sections were made to compare with the 1-D interpretation models of horizontal sections.Meanwhile, models of sandstone and mudstone of the whole area simulated directly from the SIS method were also conducted, following the same modeling workflow of each subdomain.These kinds of models, named "model (PSs)" in the workflow, underwent the same quality checking process and were compared with model (PS).The workflow of this new method is shown in Figure 4.
Variogram parameters are important factors when a two-point statistical simulation method is applied for geologic units modeling [5,21,22].In fluvial reservoir modeling, distribution trends of sediments are significantly affected by the flow direction.When predicting sandstone or mudstone dimensions, it's important and necessary to take the flow direction or source direction into account, as it has important influences on the variogram fitting, such as variogram azimuth and ranges [4,20].As shown in Figure 5, the sandstone distributions of the five subdivided zones are characterized by different trends, which are indicated by six black arrows.If the whole area is not subdivided and the standard method is used, the variogram fitting will be affected by a united north-south tendency that is not consistent with the geological knowledge (Figure 5).However, in a geologic sense, it is more reasonable to respectively capture the ranges of variables with different directions.Therefore the proposed weight coefficient-based method subtly solves the "direction problem".
Energies 2017, 10, 1439 6 of 16 Thus, the variogram parameter set for the entire study area was computed directly rather than via fitting with a spherical, Gaussian, or exponential variogram model.Finally, to evaluate the quality of the sandstone and mudstone distribution models, named "model (PS)", numerous vertical profiles of the predicted sandstone and mudstone models through the horizontal sections were made to compare with the 1-D interpretation models of horizontal sections.Meanwhile, models of sandstone and mudstone of the whole area simulated directly from the SIS method were also conducted, following the same modeling workflow of each subdomain.These kinds of models, named "model (PSs)" in the workflow, underwent the same quality checking process and were compared with model (PS).The workflow of this new method is shown in Figure 4.
Variogram parameters are important factors when a two-point statistical simulation method is applied for geologic units modeling [5,21,22].In fluvial reservoir modeling, distribution trends of sediments are significantly affected by the flow direction.When predicting sandstone or mudstone dimensions, it's important and necessary to take the flow direction or source direction into account, as it has important influences on the variogram fitting, such as variogram azimuth and ranges [4,20].As shown in Figure 5, the sandstone distributions of the five subdivided zones are characterized by different trends, which are indicated by six black arrows.If the whole area is not subdivided and the standard method is used, the variogram fitting will be affected by a united north-south tendency that is not consistent with the geological knowledge (Figure 5).However, in a geologic sense, it is more reasonable to respectively capture the ranges of variables with different directions.Therefore the proposed weight coefficient-based method subtly solves the "direction problem".

Sandstone and Mudstone Modeling for Subdivided Zones
The overall structural model was upscaled into six smaller models for the subdivided zones, and subsequently, 1-D sandstone and mudstone data were discretized into each stratigraphic framework for data analysis.The SIS method was applied to sandstone and mudstone modeling based on variograms generated from lithofacies curves that represent both the horizontal and vertical distributions of proportions for each layer [11,16].The spherical variogram model was chosen to

Sandstone and Mudstone Modeling for Subdivided Zones
The overall structural model was upscaled into six smaller models for the subdivided zones, and subsequently, 1-D sandstone and mudstone data were discretized into each stratigraphic framework for data analysis.The SIS method was applied to sandstone and mudstone modeling based on variograms generated from lithofacies curves that represent both the horizontal and vertical distributions of proportions for each layer [11,16].The spherical variogram model was chosen to simulate the range and 3-D spatial frequency of the sandstone and mudstone.The experimental variogram computation was completed based mainly on the average well spacing, and the spherical variogram model was then fitted to determine the major and minor ranges of sandstone and mudstone (Figure 6; Table 2).No other geological data, aside from the vertical sandstone and mudstone proportion curves, were used as constraints to make the simulated results reflect the original distribution patterns of sandstone and mudstone.Taking layer H8S1 in zone B as an example, Table 3 shows the variogram ranges of different directions that result from the variogram model fitting process (Figure 6).Ranges of sandstone and mudstone for layers in other zones were captured the same way.In this step, stochastic sandstone and mudstone distribution models for the six subdivided zones were established (Figure 7).simulate the range and 3-D spatial frequency of the sandstone and mudstone.The experimental variogram computation was completed based mainly on the average well spacing, and the spherical variogram model was then fitted to determine the major and minor ranges of sandstone and mudstone (Figure 6; Table 2).No other geological data, aside from the vertical sandstone and mudstone proportion curves, were used as constraints to make the simulated results reflect the original distribution patterns of sandstone and mudstone.Taking layer H8S1 in zone B as an example, Table 3 shows the variogram ranges of different directions that result from the variogram model fitting process (Figure 6).Ranges of sandstone and mudstone for layers in other zones were captured the same way.In this step, stochastic sandstone and mudstone distribution models for the six subdivided zones were established (Figure 7).2).No other geological data, aside from the vertical sandstone and mudstone proportion curves, were used as constraints to make the simulated results reflect the original distribution patterns of sandstone and mudstone.Taking layer H8S1 in zone B as an example, Table 3 shows the variogram ranges of different directions that result from the variogram model fitting process (Figure 6).Ranges of sandstone and mudstone for layers in other zones were captured the same way.In this step, stochastic sandstone and mudstone distribution models for the six subdivided zones were established (Figure 7).It is essential to verify the quality of static reservoir models to ensure that each set of variogram parameters is reliable [23,24].The vacuation checking method [25,26], in which the data of certain wells are removed during the stochastic modeling process, was applied to evaluate the quality of the sandstone and mudstone distribution models.Only the seed number was repeatedly defined to generate equiprobable realizations to compare with the 1-D interpreted models for statistical analysis.In this paper, the quality checking results for zone B are presented.The sandstone and mudstone data of the wells shown in Figure 8(a2) were used to generate prediction models (A1, A2, . . ., A5 and B1, B2, . . ., B5), which were then used for comparison with the 1-D interpreted models, shown in the leftmost panels of Figure 8c,d.All of the prediction models are the results of using same variogram parameter sets, which are based on the original well pattern shown in Figure 8(a1).According to this comparison, the prediction distributions of interlayer A show more favorable results than those of interlayer B, which show poor connectivity between wells.The modeled distributions of sandstone are largely in accordance with the 1-D interpreted model, whereas predictions of the thin sets of sandstone are not well matched (Figure 8c,d).It is essential to verify the quality of static reservoir models to ensure that each set of variogram parameters is reliable [23,24].The vacuation checking method [25,26], in which the data of certain wells are removed during the stochastic modeling process, was applied to evaluate the quality of the sandstone and mudstone distribution models.Only the seed number was repeatedly defined to generate equiprobable realizations to compare with the 1-D interpreted models for statistical analysis.In this paper, the quality checking results for zone B are presented.The sandstone and mudstone data of the wells shown in Figure 8(a2) were used to generate prediction models (A1, A2, …, A5 and B1, B2, …, B5), which were then used for comparison with the 1-D interpreted models, shown in the leftmost panels of Figure 8c,d.All of the prediction models are the results of using same variogram parameter sets, which are based on the original well pattern shown in Figure 8(a1).According to this comparison, the prediction distributions of interlayer A show more favorable results than those of interlayer B, which show poor connectivity between wells.The modeled distributions of sandstone are largely in accordance with the 1-D interpreted model, whereas predictions of the thin sets of sandstone are not well matched (Figure 8c,d).In general, realizations from vacuated well data largely coincide with the 1-D interpreted data, which implies that the variogram parameter sets are reliable.The variogram parameter sets of other zones were analyzed in the same way.In this step, reasonable variogram parameter sets for each zone were obtained (Table 3).As shown in the Figure 9a,b, the ranges of the same layers of different subzones vary greatly.The variation trend represented by the dotted arrow indicates that the spatial correlation of sand bodies developed in a sandy braided river system is relatively higher when compared with other meandering intervals (H8S1, H8S2, S1-1, S1-2, and S1-3); however, their geological properties are not discussed in this paper.More details about their depositional system could refer to Zhao et al. [16] and Wang et al. [27].This kind of correlation is consistent with the high sandstone: mudstone ratio expressed by the proposed geologic model.Figure 9c directly shows the strong heterogeneity of well numbers and well densities in six subzones.The subzone B has the highest well density among all of the subzones (Figure 9c), and both the major-direction ranges of sandstone and interlayer A are smaller than those of other subzones.Therefore, it could be inferred that variogram ranges are affected by the well patterns, and the proposed geostatistical method is meaningful.
In general, realizations from vacuated well data largely coincide with the 1-D interpreted data, which implies that the variogram parameter sets are reliable.The variogram parameter sets of other zones were analyzed in the same way.In this step, reasonable variogram parameter sets for each zone were obtained (Table 3).As shown in the Figure 9a,b, the ranges of the same layers of different subzones vary greatly.The variation trend represented by the dotted arrow indicates that the spatial correlation of sand bodies developed in a sandy braided river system is relatively higher when compared with other meandering intervals (H8S1, H8S2, S1-1, S1-2, and S1-3); however, their geological properties are not discussed in this paper.More details about their depositional system could refer to Zhao et al. [16] and Wang et al. [27].This kind of correlation is consistent with the high sandstone: mudstone ratio expressed by the proposed geologic model.Figure 9c directly shows the strong heterogeneity of well numbers and well densities in six subzones.The subzone B has the highest well density among all of the subzones (Figure 9c), and both the major-direction ranges of sandstone and interlayer A are smaller than those of other subzones.Therefore, it could be inferred that variogram ranges are affected by the well patterns, and the proposed geostatistical method is meaningful.

Variogram Parameter Computations and Applications
The variogram ranges (major direction, minor direction, and vertical direction) for the entire S6 area were computed using Equation (3), and the results are shown in Table 4.All of the parameters were directly applied to generate the sandstone and mudstone distribution models for the entire S6 area (Figure 10).In this step, the spherical variogram model is used, and the seed numbers are stochastic.

Variogram Parameter Computations and Applications
The variogram ranges (major direction, minor direction, and vertical direction) for the entire S6 area were computed using Equation (3), and the results are shown in Table 4.All of the parameters were directly applied to generate the sandstone and mudstone distribution models for the entire S6 area (Figure 10).In this step, the spherical variogram model is used, and the seed numbers are stochastic.

Model Evaluation
The quality of models (PS) was evaluated with the assistance of horizontal wells, for which the 1-D sandstone and mudstone data were not upscaled during the stochastic modeling process.

Model Evaluation
The quality of models (PS) was evaluated with the assistance of horizontal wells, for which the 1-D sandstone and mudstone data were not upscaled during the stochastic modeling process.Horizontal wells in the S6 area are mainly designed to produce gas from the H8X member, and the production data show that all of these horizontal wells are capable of producing industrial gas flow.The 1-D sandstone and mudstone models of these horizontal wells were interpreted using the only available electrical logging curve, gamma ray (GR), which was captured by the logging while drilling technology.A total of 50 cross-sections of predicted sandstone and mudstone models for each horizontal well were made to compare with the corresponding 1-D interpreted sandstone and mudstone models, and three classes (I, II and III) of the comparison results were defined to evaluate the prediction models (Table 5; Figure 11).For ease of analysis, if a horizontal interval is dominated by sandstone (mudstone), the mudstone (sandstone) part is ignored when it comes to the evaluation of class II.
Finally, the numbers of each class of each well were collected to analyze model quality (Figure 12).Comparison between the 1-D interpreted model and models (PS) of well SH1 show that the models of type I and II account for a dominant proportion of the prediction models (Figure 12), which reflect the relatively high reliability of the variogram parameters obtained from this proposed method.Meanwhile, the same comparison for models (PSs) of well SH1 was conducted (Figure 13).As a result, types I and II are also dominant.However, the average frequencies of type I, II, and III for models (PS) and models (PSs) are different.As shown in Table 6, the average frequencies of type I and II of models (Ps) are higher than those of models (PSs), while the average frequency of type III of models (Ps) are lower than that of models (PSs).To some extent, the variogram parameters obtained from the proposed method are proved to be more reliable than the standard method.
Horizontal wells in the S6 area are mainly designed to produce gas from the H8X member, and the production data show that all of these horizontal wells are capable of producing industrial gas flow.The 1-D sandstone and mudstone models of these horizontal wells were interpreted using the only available electrical logging curve, gamma ray (GR), which was captured by the logging while drilling technology.A total of 50 cross-sections of predicted sandstone and mudstone models for each horizontal well were made to compare with the corresponding 1-D interpreted sandstone and mudstone models, and three classes (I, II and III) of the comparison results were defined to evaluate the prediction models (Table 5; Figure 11).For ease of analysis, if a horizontal interval is dominated by sandstone (mudstone), the mudstone (sandstone) part is ignored when it comes to the evaluation of class II.
Finally, the numbers of each class of each well were collected to analyze model quality (Figure 12).Comparison between the 1-D interpreted model and models (PS) of well SH1 show that the models of type I and II account for a dominant proportion of the prediction models (Figure 12), which reflect the relatively high reliability of the variogram parameters obtained from this proposed method.Meanwhile, the same comparison for models (PSs) of well SH1 was conducted (Figure 13).As a result, types I and II are also dominant.However, the average frequencies of type I, II, and III for models (PS) and models (PSs) are different.As shown in Table 6, the average frequencies of type I and II of models (Ps) are higher than those of models (PSs), while the average frequency of type III of models (Ps) are lower than that of models (PSs).To some extent, the variogram parameters obtained from the proposed method are proved to be more reliable than the standard method.

Discussion
Weight coefficient-based methods of stochastic modeling for reservoirs have rarely been reported.The gentle regional structural characteristics, the sandy braided river depositional system, and the strong heterogeneity of the well distribution are all important factors that make this method feasible for this study.In terms of local tectonics, the undeveloped fault system in this region results in good structural integrity, which greatly reduces the structural uncertainty among the wells.The sandy depositional sedimentary system with a stable provenance direction is also favorable for the indicator simulation method of lithofacies.The strong heterogeneity of well distribution greatly affects the variogram model-fitting process, and to some extent, applying weight coefficients can weaken or counterbalance the influence of this type of heterogeneity on predictions.The randomly drilled horizontal wells act as an opportune evaluation tool to enable the quantitative judgment of variogram parameters.For further development in the immense Sulige gas field, this novel method can provide meaningful guidance for reservoir prediction based on 3-D geological modeling technology.Moreover, this new method is also lithofacies modeling in areas with greater fault development, under the condition that the fault systems formed after the reservoir was shaped.
However, this method still has several disadvantages.(1) Workloads are multiplied.Modeling for the same layers but for multiple zones makes the simulation procedure considerably more time-consuming.(2) Quality checking of sandstone and mudstone models for the subdivisions of zones is strictly qualitative, and the results therefore lack scientific rigor.Furthermore, the vacuation checking method only allows a small number of wells to be removed, to keep the generated models meaningful, which results in limited 1-D interpreted data available for comparison with the model predictions.Especially for zones like zone F, few wells are present, and therefore such quality checking methods are not ideal for this zone.
Because of the low porosity and permeability, and strong heterogeneity in the study area, several challenges remain for predicting sandstone and mudstone distributions in the S6 area.(1) The 3-D seismic data of the S6 area have limited coverage and, and because of the low quality of the data, they cannot be used to recognize lithofacies and distributions.(2) The quality of the static sandstone and mudstone models could hardly be effectively evaluated with the regular dynamic analysis because of the tightness of the reservoir.For example, in terms of reservoir connectivity, the effective discharge area of each well could be quantitatively evaluated using the numerical simulation method; however, the simulated results are based on the practical production data, which are controlled by the limitations of existing fracturing technology [28].Therefore, the effective discharge area calculated with the numerical simulation does not in fact represent the inherent reservoir connectivity.In other words, the effective discharge area of a well may become larger if the fracturing technology becomes more advanced.Publications on hydraulic fractures or cracks could be referred to Rabczuk and Belytschko [29], Zhuang et al. [30], Ma et al. [31], Wang et al. [32], and Nguyen et al. [33].Other dynamic analysis methods, such as interference tests implemented in the infill well area, face similar challenges.
(3) Comparison between prediction models and 1-D interpreted models is conducted manually and semi-quantitative; the development of automatic and quantitative approaches should be prioritized in the future.

Conclusions
The 3-D distribution of sandstone and mudstone was modeled using the proposed weight coefficient method based on the SIS method.This new method is specifically proposed for the S6 area, where there are no faults and where wells are irregularly distributed, to balance the heterogeneity of well density for more reasonable variogram parameters.The primary variogram parameters of the six subdivided zones were determined based on the variogram model fitting, and quality controlled using the vacuation checking method.Then, 1-D interpreted sandstone and mudstone models of horizontal wells were used as standards for comparison between these models and the equiprobable realizations generated from the weight coefficient variogram parameters.The results show that types I and II are predominant among the prediction models, which indicates that the variogram parameters computed using the weight coefficient method are reliable.

Figure 1 .
Figure 1.(a) Location map of S6 area in Sulige gas field; (b) Stratigraphy of upper Paleozoic; only the H8 and S1 layers are shown.

Figure 1 .
Figure 1.(a) Location map of S6 area in Sulige gas field; (b) Stratigraphy of upper Paleozoic; only the H8 and S1 layers are shown.

Figure 2 .
Figure 2. (a) Classification of the subdivided zones, showing the locations of directional wells and horizontal wells.While subdividing, the boundaries of the six zones are roughly drawn; (b) Discretized sandstone/mudstone data used as input data for stochastic modeling.

Figure 3 .
Figure 3. (a) 1-D interpreted sandstone and mudstone data of four directional wells; (b-d) Statistical histograms of sandstone, interlayer A, and interlayer B, respectively.Note that the mean thicknesses of interlayer A and interlayer B are quite different.

Figure 2 .
Figure 2. (a) Classification of the subdivided zones, showing the locations of directional wells and horizontal wells.While subdividing, the boundaries of the six zones are roughly drawn; (b) Discretized sandstone/mudstone data used as input data for stochastic modeling.

Figure 2 .
Figure 2. (a) Classification of the subdivided zones, showing the locations of directional wells and horizontal wells.While subdividing, the boundaries of the six zones are roughly drawn; (b) Discretized sandstone/mudstone data used as input data for stochastic modeling.

Figure 3 .
Figure 3. (a) 1-D interpreted sandstone and mudstone data of four directional wells; (b-d) Statistical histograms of sandstone, interlayer A, and interlayer B, respectively.Note that the mean thicknesses of interlayer A and interlayer B are quite different.

Figure 3 .
Figure 3. (a) 1-D interpreted sandstone and mudstone data of four directional wells; (b-d) Statistical histograms of sandstone, interlayer A, and interlayer B, respectively.Note that the mean thicknesses of interlayer A and interlayer B are quite different.

Figure 4 .
Figure 4. Workflow of the weight coefficient method of modeling sandstone and mudstone distribution.

Figure 4 .
Figure 4. Workflow of the weight coefficient method of modeling sandstone and mudstone distribution.

Figure 5 .
Figure 5. Different types of sand bodies and mudstone distribution trends of the subdivided zones.The different trends are indicated by six black arrows and the blue arrow refers to the general source direction of the study area.F1 and F2 refer to the sandy facies, F3 and F4 refer to the muddy facies.

Figure 5 .
Figure 5. Different types of sand bodies and mudstone distribution trends of the subdivided zones.The different trends are indicated by six black arrows and the blue arrow refers to the general source direction of the study area.F1 and F2 refer to the sandy facies, F3 and F4 refer to the muddy facies.

Figure 6 .
Figure 6.Variogram model fitting of sandstone and mudstone for layer H8S1 in zone B.

Figure 7 .
Figure 7. 3-D sandstone and mudstone distribution models for the six subdivided zones.It's noted that each distribution model is just one of the equiprobable realizations.

Figure 6 .
Figure 6.Variogram model fitting of sandstone and mudstone for layer H8S1 in zone B.

Figure 6 .
Figure 6.Variogram model fitting of sandstone and mudstone for layer H8S1 in zone B.

Figure 7 .
Figure 7. 3-D sandstone and mudstone distribution models for the six subdivided zones.It's noted that each distribution model is just one of the equiprobable realizations.Figure 7. 3-D sandstone and mudstone distribution models for the six subdivided zones.It's noted that each distribution model is just one of the equiprobable realizations.

Figure 7 .
Figure 7. 3-D sandstone and mudstone distribution models for the six subdivided zones.It's noted that each distribution model is just one of the equiprobable realizations.Figure 7. 3-D sandstone and mudstone distribution models for the six subdivided zones.It's noted that each distribution model is just one of the equiprobable realizations.

Figure 8 .
Figure 8.(a) The original (the upper one) and the vacuated well pattern (the lower one), showing the locations of eight removed wells; (b) Cross-section of sandstone and mudstone model, within which well A and well B did not participate in the modeling process; section location is drawn in Figure 8(a2); (c,d) Comparison between five stochastic realizations or prediction models (A1, A2, …, A5 and B1, B2, …, B5) and 1-D interpreted sandstone and mudstone model of well A and well B, respectively.

Figure 8 .
Figure 8.(a) The original (the upper one) and the vacuated well pattern (the lower one), showing the locations of eight removed wells; (b) Cross-section of sandstone and mudstone model, within which well A and well B did not participate in the modeling process; section location is drawn in Figure 8(a2); (c,d) Comparison between five stochastic realizations or prediction models (A1, A2, . . ., A5 and B1, B2, . . ., B5) and 1-D interpreted sandstone and mudstone model of well A and well B, respectively.

Figure 9 .
Figure 9. Relationships between layers and major direction ranges of different subzones.As shown in the figure, both the major direction ranges of sandstone and interlayer A for the six subzones are characterized by obvious horizontal and vertical fluctuations.(a) The ranges of sandstone for the braided intervals H8X1 and H8X2 are relatively larger than those for the meandering intervals (H8S1, H8S2, S1-1, S1-2 and S1-3); (b) The ranges of interlayer A for the braided intervals are relatively smaller than those for the meandering intervals; (c) Statistics of well numbers and well densities for the six subzones.

Figure 9 .
Figure 9. Relationships between layers and major direction ranges of different subzones.As shown in the figure, both the major direction ranges of sandstone and interlayer A for the six subzones are characterized by obvious horizontal and vertical fluctuations.(a) The ranges of sandstone for the braided intervals H8X1 and H8X2 are relatively larger than those for the meandering intervals (H8S1, H8S2, S1-1, S1-2 and S1-3); (b) The ranges of interlayer A for the braided intervals are relatively smaller than those for the meandering intervals; (c) Statistics of well numbers and well densities for the six subzones.

Figure 11 .
Figure 11.Classified cross-sections of prediction models for horizontal well SH1.The uppermost figure is a cross-section of a simulated sandstone and mudstone model, and the lower figure next to it is the 1-D interpreted sandstone and mudstone model.The lowest cross-sections are eight stochastic realizations (R1, R2, ..., R8) generated via the same variogram parameters, but using different seed numbers.The location of horizontal well SH1 is shown in Figure 2a.Lsm (Lmm) refers to the length of sandstone (mudstone) obtained from stochastic modeling; Lsi (Lmi) refers to the length of sandstone (mudstone) obtained from well logging interpretation.

Figure 11 .
Figure 11.Classified cross-sections of prediction models for horizontal well SH1.The uppermost figure is a cross-section of a simulated sandstone and mudstone model, and the lower figure next to it is the 1-D interpreted sandstone and mudstone model.The lowest cross-sections are eight stochastic realizations (R1, R2, ..., R8) generated via the same variogram parameters, but using different seed numbers.The location of horizontal well SH1 is shown in Figure 2a.Lsm (Lmm) refers to the length of sandstone (mudstone) obtained from stochastic modeling; Lsi (Lmi) refers to the length of sandstone (mudstone) obtained from well logging interpretation.

Figure 12 .
Figure 12.Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the proposed method.As is shown in the Figure, type I and II are dominant.

Figure 13 .
Figure 13.Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the standard SIS method.Frequencies of type I and II are similar to those of models (PS).

Figure 12 .
Figure 12.Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the proposed method.As is shown in the Figure, type I and II are dominant.

Figure 13 .
Figure 13.Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the standard SIS method.Frequencies of type I and II are similar to those of models (PS).

Table 1 .
Well densities of six zones and the whole study area.

Table 1 .
Well densities of six zones and the whole study area.

Table 2 .
Data used to model distribution of sandstone and mudstone for layer H8S1 in zone B.

Table 2 .
Data used to model distribution of sandstone and mudstone for layer H8S1 in zone B.

Table 4 .
Computed variogram ranges for the entire zone.

Table 4 .
Computed variogram ranges for the entire zone.

Table 5 .
Evaluation criterion for prediction model.

Table 6 .
Average frequencies of type I, II, and III for models (PS) and models (PSs).Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the proposed method.As is shown in the Figure, type I and II are dominant.

Table 5 .
Evaluation criterion for prediction model.

Table 6 .
Average frequencies of type I, II, and III for models (PS) and models (PSs).Statistics of comparison results of 48 horizontal wells in S6 area.These comparison results are based on the standard SIS method.Frequencies of type I and II are similar to those of models (PS).

Table 5 .
Evaluation criterion for prediction model.

Table 6 .
Average frequencies of type I, II, and III for models (PS) and models (PSs).