Facies and Petrophysical Modeling of Triassic Chang 6 Tight Sandstone Reservoir, Heshui Oil Field, Ordos Basin, China

A lot of research has been done on tight sandstone reservoirs by many authors around the globe. Usually, the oil in a tight sandstone reservoir with permeability less than 0.1 mD and porosity less than 10% can be referred to as tight (sandstone) oil [1–9]. Recently, tight (sandstone) oil emerged as an important unconventional oil resource in the petroleum sector. Accumulations of tight oil exist in China in different basins, especially in the Chang 6 reservoir located in the Ordos Basin. Most of the reservoirs in the Ordos Basin are characterized by ultralow permeability with substantial heterogeneity. The current research is focused on reservoir modeling that has become a unique approach to access the reservoir quality, characterization, reserve estimation enhancement, and field development, future production prediction, additional wells placements, and assessment of varying reservoir management situations [10].

The tight reservoir modeling emerged as a new frontier in the geoscience field. A lot of new development is underway, and recently, the new modeling techniques are being used. As the conventional and unconventional reservoirs are different in many other characteristics, there is also difference in related to the modeling aspect. The tight reservoir modeling is different from conventional reservoir modeling mainly because of reservoir characteristics, distribution of facies, scale, heterogeneity, etc. While modeling tight reservoir attention should be paid to honor fine-scale heterogeneities, ultralow petrophysical properties, facies distribution, facies associations (abundance, trend, stacking pattern of facies, interfacies relationships, mutual contacts, extent along with different directions, and geometry), and dimension of geobodies.

Previous studies have reported that the use of 3D computer-based reservoir models plays a vital role in increasing recovery, management of subsurface reservoirs, and optimization of drainage strategy [11–15]. The reservoir modeling mainly incorporates gathering the data, building structural and stratigraphic framework, 3D grids generation, then populating the grid cells with the petrophysical properties. Mostly these petrophysical properties are interlinked with facies attributes and models; those could be further integrated with the conceptual 3D geological model. The facies models also provide a better approach for the representation of petrophysical properties. The models usually incorporate a variety of data sets of the understudied reservoir. These data sets usually include geophysical, geological, production, and reservoir engineering data [16–19]. This data can be further used for facies and petrophysical modeling that is important for concepts of petrophysical and facies properties of the reservoir. Understanding of the structural, stratigraphic, petrophysical, and facies configuration is the application of 3D modeling.

The main aim of establishing the facies models is to understand the characteristics of different environments, including scale heterogeneity, and controlling physical processes [20]. The essential step in constructing a facies model is to fully understand facies, facies architecture (which includes facies-stacking pattern and cyclicity), and diagenesis [21]. Sedimentary facies architecture or spatial distribution [22–25] in clastic depositional systems and petrophysical rock fabric distribution in carbonate reservoirs [26, 27] is the most important parts of facies models. In petroleum reservoirs, all these factors, from a small scale such as facies to a large scale, for example, a depositional system [28], can be analyzed.

The current study incorporates research on the Chang 6 member of the (Upper Triassic) Yanchang Formation, which is one of the significant (oil-bearing) formations in the Ordos Basin. A lot of work has been done by previous authors on Chang 6 related to the sedimentary facies identification [29], diagenesis, and diagenetic facies [30, 31] sand body structural characteristics [32]. The main aim of this study is facies and petrophysical modeling of tight sandstone reservoirs of the target member. Well data was utilized for the interpretation, modeling of facies, and petrophysical parameters (permeability, porosity, and water saturation). At the reservoir scale, these parameters control the fluid flow characteristics and reservoir heterogeneities; therefore, reservoir heterogeneity is scale-dependent as well. Different methods were applied for facies modeling such as indicator kriging, sequential indicator simulation, and multiple-point geostatistics. These powerful methods have made facies modeling a vital technique for exploration in the petroleum industry. The parameters derived from facies models are provided as inputs in subsequent reservoir simulations. The petrophysical behavior of the reservoir can be analyzed by modeling the effect of reservoir properties in the three-dimensional structure [33]. Gaussian random function simulation (GRFS) and sequential Gaussian simulation (SGS) were applied for petrophysical modeling, including permeability, porosity, and water saturation. Precise assessment of petrophysical properties is an important part of reservoir management, and it could be performed by utilizing principles of geostatistics [34]. The hydrocarbon flow depends upon the facies distribution in the reservoir [35]. Choosing the best method for facies modeling requires the implication of the methods and analyzes the differences and similarities in the resulting models.

The Heshui area is not intensely deformed, and geological structures present in this area are gently dipped with dipping angle less than 1° and with the exception of the low-amplitude nose-like bulge that is present in a small region on the west-plunge monocline. Lithological traps are their main exploration targets as the secondary tectonic belt and structural traps lacking oil and gas enrichment. Geostatistical studies of the stratigraphic sequences are requisite tools in making fair predictions of the reservoir properties away from well control points and thereby aiding the exploitation of the huge potential of the reservoirs. A reservoir model needs to be estimated using various combinations of property modeling algorithms for discrete (facies) and continuous (petrophysical) properties. Uncertainties associated with these algorithms, as well as the sensitivity of some input data which could aid the accuracy of the computed volume, will also be estimated. 3D reservoir modeling of the Chang 6 reservoir is required to tackle the complicated problems related to reservoir heterogeneity.

Ordos Basin is in central China, spanning the five provinces of Shaanxi, Gansu, Ningxia, Shanxi, and Mongolia, which covers an area around $25×104km2$ (Figure 1(a)). The Ordos Basin experienced a fundamental transformation from the marine phase to the transitional phase of the continental facies and eventually evolved into it is a large freshwater inland sedimentary lake basin with a mild climate and dense vegetation. It is a multistage cratonic basin with the platform, Taiyuan depression, and Taizhong depression. It has geological features of obvious twisting, depression migration, and structural stability settlement, which can be further divided into six first-order structures: Jinxi Flexural Fold Belt, Tianhuan Depression, West Marginal Thrust Belt, Yishan Slope, Weibei Fault Fold Belt, and Yimeng Uplift Belt (Figure 1(b)). The locality of the study area, Heshui oil field, is in the southwestern part of the Yishan slope of the Ordos Basin which is positioned in the central part of China (Figure 1(b)). The study area is located in Qingcheng and Heshui County of Gansu Province, with an exploration area up to 3,000 km².

The sandstone reservoirs are tight and heterogeneous in this area [8]. According to the vertical distribution law of the reservoir, it is further subdivided into ten members from Chang 1 to Chang 10. The sedimentary environment was mainly dominated by sublacustrine fan deposits at the time of the deposition of Chang 6 member. Chang 6 reservoir of the Triassic Yanchang Formation has a depth of 1428 m ~2000 m and a thickness of 110 m ~130 m. The porosity of the Chang 6 reservoir is 8.56% ~12.3%, and the permeability is less than 1 mD, which belongs to the tight reservoir.

At the time when Chang 6 member was depositing, the rate of sedimentation in the basin was slowed down at a significant level and then started to shrink gradually. Delta and sublacustrine fan deposits of the Yanchang Formation (Ordos Basin) are the two major sedimentary systems developed in the northeast and southwest, that is why this period is considered as an important period of delta deposits construction. The upper section of Chang 6 consists of interbedding of dark-gray fine sandstone and gray-black mudstone, and the lower part is gray-black mudstone with black carbonaceous mudstone (Figure 2). The sedimentary structures in Chang 6 mainly include parallel bedding, massive bedding, trough, plate, and wedge-shaped cross-bedding. The Chang 6 reservoir is marked as important oil-producing layers in Ansai, Xifeng, Ganguyi, and Heshui oil fields due to the sedimentation during the Kuanfeng period of basin evolution.

Based on previous analysis and research [38], sublacustrine fan sedimentary facies is developed in the Chang 6 reservoirs in the study area (Figure 1(b)), and the Chang 6 is further subdivided into three units, namely, Chang 6₁, Chang 6₂, and Chang 6₃. The lithology of these units is dominated by fine sandstone, dark gray mudstones, and siltstones [38] (Figure 2). The sedimentary facies of the channel (main channel and distributary channel), lobe body (lobe main body and lobe margin), interlobe, and interchannel facies have been identified on the basis of well logs identification (Figure 3).

There is a variation in the petrophysical properties of each facies analyzed on the basis of core analysis. In the distributary channel, the porosity ranges from 4.5% ~9%. While in lobe main body facies, the porosity is 3% ~8%. There is 1% ~3.5% porosity in the lobe margin facies. In the distributary channel, the permeability ranges from 0.1 ~ 0.7 mD. However, the permeability range in the lobe main body is from 0.1 mD~0.55 mD. The range of permeability in the lobe margin facies is from 0.01 mD~0.1 mD.

The lithology of the main channel is dominated by fine sandstone and coarse siltstone; the thickness of the main channel is generally more than 3 meters, and finning up grain size; bell-shaped GR values on logging curves, with sudden change at the bottom; the main channel facies mainly develops sedimentary structures such as trough, plate, and wedge-shaped cross-bedding, with the obvious erosional surface at the bottom; the deposition mechanism is the coexistence of traction current and gravity current (Figure 4).

The lithology of the distributary channel is dominated by fine sandstone and siltstone with positive grain size and thickness generally greater than 2 meters. On logging curves, GR curves are toothed bell-shaped with sudden change at the bottom. Block bedding, parallel bedding, and plate cross-bedding can be seen in the core, and the prominent erosional surface can be seen at the bottom. The sedimentary mechanism is the coexistence of traction current and gravity current (Figure 5).

According to the location and thickness of the lobe sand body, lobe sand body can be divided into the lobe main body and lobe margin.

The main lithology of the lobe main body is fine sandstone and coarse siltstone, and the grain size is counter-rhythm or homogeneous rhythm. The thickness of the sand body is generally greater than 2 m. The GR curve is dentate funnel-shaped or box-shaped. The sedimentary structure is dominated by massive bedding with cross-bedding. The sedimentary mechanism is sandy debris flow or upper turbidity flow with lower sandy debris flow (Figure 6).

The lithology of the lateral margin of lobes is mainly argillaceous siltstone with no obvious grain size sorting, and its thickness is about 1-2 meters; the GR curve is finger-like; massive bedding and parallel bedding are mainly developed (Figure 7).

The interchannel/interlobe lithology is dominated by gray-black mudstone and argillaceous siltstone, with occasional thin sand bodies less than 1 m and no obvious rhythm. On logging curves, GR is near the mudstone baseline (Figures 7 and 8).

The data for the present study include core data of 10 wells (900 m); well data of 208 wells with an average distance of about 250 m; analytical test data such as porosity, permeability, and oil saturation data; and abundant logging data such as GR, SP, DEN, and RT logs. This data has been used to obtain structural and stratigraphic framework that provided basis for facies and petrophysical modeling. The workflow adopted to model the Chang 6 tight sandstone reservoir in Heshui area incorporates structural modeling, facies modeling, and petrophysical modeling based on well log interpretation and utilization during modeling process (Figure 9).

Spatial data analysis is important step in the reservoir modeling workflow, it consists of analyze the data, quality control (QC), and provide inputs for facies and petrophysical modeling [10, 40]. Input data transformation, variograms interpretation, and identification of trends are utilized to identify the variation facies (discrete reservoir properties) and porosity, shale volume, permeability, and water saturation (continuous reservoir properties) with respect to spatial continuity. Facies modeling was used for distribution of discrete facies, clear idea of facies stratigraphic framework, and reservoir architecture throughout the modeling grid.

The lithology interpretation of Chang 6 reservoir is based on the gamma ray, resistivity, and AC logs. Petrel (modeling software by Schlumberger) was used in the present study. That involves all the major step taken to complete this research, and the data was loaded into petrel which included well heads, well paths, well tops, and well logs. Then, stratigraphic modeling was initiated consisted of marking well tops and well correlation. In structural modeling, the target of making six horizons (Chang 6_2_1, Chang 6_2_2, Chang 6_3_1, Chang 6_3_2, Chang 6_3_3, and Chang 7_1), five zones (zone 1, zone 2, zone 3, zone 4, and zone 5), and layers (layers thickness varied for each zone) was achieved. Then logs were upscaled in order to use it for initiating facies modeling process.

A stochastic (pixel-based) approach was used to model the facies. This method enabled us to individually set volume fractions and variograms for each facies and distribute channel in a model using fluvial modeling and facies objects using object modeling [41]. It was also utilized to allocate values to the facies that existed in other model cells. Where sparse data is present, this technique is usually utilized. A 3D realistic static model is provided by the accurate distribution of facies with the well data [16]. These models provide an efficient approach to understand reservoir heterogeneity. The SIS method commonly uses a spherical variogram that produces a linear behavior at small separation distances. The variogram is defined by its: type, range, sill, and nugget. The range indicates where the variogram model approaches its plateau. The plateau refers to the point where increasing separation distance no longer changes the correlation degree among data point pairs. Variograms are often measured along the horizontal plane (major and minor directions) and vertical plane. The main direction determines the direction where the highest correlation was found among the data points, whereas the minor direction is defined as perpendicular to the main direction [40]. The sill determines the semivariance where the distance of separation is more than the range. Semivariance refers to the average of the squared deviations of values between data points. The nugget determines the semivariance, where the separation distance is zero [40].

Multiple-point geostatistics is another stochastic simulation method, but instead of using a two-point statistical variogram. MPG derives the local conditional distribution from a training image [42]. Multiple-point geostatistics was designed for the reproduction complex curvilinear features from a training image, but retain the stochastic flexibility to condition hard and soft data [40]. That allows the modeler the capability to capture curvilinear geometries from a user-constructed training image, and then anchor these geometries to hard data. The MPG modeling method is capable of reproducing the observed facies association geometries. Therefore, the MPG method could accurately model the three-dimensional connectivity and continuity of productive and nonproductive reservoir facies. The workflow of the multiple-point geostatistics method consists of the following steps:

• (1)

  Making a training image grid or working on a previous grid (that can be quite large and take a long time to run)

• (2)

  Building a training image (e.g., using object modeling, surface-based modeling, interactive drawing, and user-defined object creation)

• (3)

  Pattern creation by applying the search mask and multigrid concept

• (4)

  Testing the pattern that is optional, we can either test or apply the pattern directly, but it is recommended to test the pattern before applying to the real model, and it can save time and give results more effectively

• (5)

  Facies modeling is the final step that utilizes training image, pattern, rotation/scaling, and well conditioning

Kriging was built on the concept to use the distance-weighted algorithm in order to predict values at unknown points and also using known points. Kriging is usually defined as a stochastic algorithm that calculates the spatial autocorrelation between measured points using a semivariogram, with respect to the distance between them [43]. Then, optimal weights can be determined through that semivariogram. Therefore, it is recognized as one of the best unbiased linear estimators. Kriging is considered as the best technique with large sample populations, and in some particular situations, it can also handle the data that is irregularly spaced data [44]. In the case of discontinuous data sets, kriging does not perform good, e.g., data set discontinuities related to geological modeling even with good data support [45]. Even though kriging remains one of the most widely used geological interpolation methods, and also for interpolating surfaces from borehole data [46–48]. Indicator kriging (IK) is known as a method of spatial interpolation, used to estimate the cumulative conditional distribution function (ccdf) of the variable at an unsampled spot. The result should be a vector that is obtained from the approximation of conditional cumulative distribution function and its corresponding discrete probability density function (cpdf), where the probability of an occurrence of a class is given by each component. After interpolation of these log-odds, the preferred ccdf or cpdf can be conveniently rerepresented.

Model building using a surface-based method may boost up the recent developments in reservoir modeling techniques [49]. The surface-based model building approach shows main stratigraphic surfaces and boundaries of facies as surfaces which can be interpolated between control lines or points deterministically or integrate a stochastic dimension, with sparse control data [50]. This method is based on the anatomical results of reservoir architecture. In this method, the plane distribution of reservoir architecture is used as the plane constraint, and the thickness distribution of reservoir sand body is used as the vertical constraint. And then, the morphological features of each architecture element are combined to perform deterministic modeling of three-dimensional architecture, with emphasis on reproducing the results of reservoir architecture anatomy. The method is guided by multidimensional interaction of single wells, planes, and sections, and the reservoir architecture model is established hierarchically. First, the distribution characteristics of the three-dimensional spatial architecture elements in the study area are reproduced. Then, depending on the interpretation of single well and plane distribution of the argillaceous interlayer, a three-dimensional model of the argillaceous interlayer is established to reproduce the distribution characteristics of each small layer of sand bodies. This method varies from existing modeling approaches used to describe stratigraphic zones in a reservoir using surfaces. A grid is constructed within each zone to model the facies architecture on it.

The use of a sequential Gaussian simulation (SGS) method is known as one of the common methods for determining petrophysical parameters. Due to its simplicity and flexibility in generating real heterogeneities, this approach is used frequently for modeling. The sequential simulation can be performed in 5 steps: (1) main data conversion to a new space, (2) modeling variogram in new space, (3) stochastic route determination to identify the locations lacking the sample, (4) using an alternative way for the estimation of places without a sample, and (5) simulated values conversion in a reverse way [51]. The key concept of SGS is so straightforward. Kriging offers an approximation of the mean and standard deviation of a variable for each grid node, which provides an easy approach to classify the variable by a standard (Gaussian) distribution as a random variable. Rather than taking the mean to estimate every node, SGS chooses the random deviate from the normal distribution chosen by a uniform random number, which reflects the degree of probability. The flow chart for sequential Gaussian simulation is shown in Figure 10. The petrophysical properties can be simulated using SGS. Then, data analysis can be conducted on each property. Before sequential Gaussian simulation porosity goes through a normal score transformation and prior to normal score transformation permeability underwent a logarithmic transformation. SGS is utilized to describe the distribution of the petrophysical properties inside stratigraphic layers and facies in the reservoir.

The algorithm used primarily for simulation of continuous property in petrel is Gaussian random function simulation (GRFS) [53]. This geostatistical method is performed based on kriging [40] in order to determine interpolated values and then add an unconditional simulation to the kriging estimates to simulate a realization conditional to the available data. This method is implemented in Petrel 2015, is computationally efficient, and produces results that are comparable to certain methods of geostatistical simulation including sequential Gaussian simulation (SGS) [40].

All kriging-based methods require the spatial structure of the modeled variable to be quantified. The spatial structure is represented by a variogram, which is modeled by calculating an experimental variogram that quantifies the difference in values at specified distances, called lags. A variogram is calculated in the minor and major directions of spatial correlation and also the vertical direction. Each direction is modeled mostly independently of the others, with the weight or sill values of the different model functions being the only interaction between them.

There are different geostatistical and stochastic methods for facies modeling. These methods include stochastic and geostatistic modeling (Table 1). For facies modeling, the following stochastic modeling algorithms can be used: Variogram-based geostatistics modeling, object-based stochastic modeling, sequential indicator simulation, multiple-point geostatistics modeling, and deterministic modeling algorithms include indicator kriging. The key to using this method is based on the well logs upscaled properties and data analysis in various methods. The variogram parameter settings also play important role in getting the best results, and getting results for best experimental variogram have its own importance. There are three possible ways in which best variogram parameters can be interpreted, e.g., from data analysis in the variogram settings, variograms can be adjusted, the second way involves using variogram maps, and third, it is possible to interpret variogram parameters from previous model built using other techniques/available maps. In variogram analysis variogram type (Spherical, Exponential, and Gaussian), sill, major range, minor range, vertical range, azimuth, and dip calculations are important. In current study, all these methods have been utilized to get the best results for facies models using different methods.

(1) Sequential Indicator Simulation (SIS). Sequential indicator simulation (SIS) is a stochastic, pixel-based, modeling method. Stochastic simulation methods use the same set of data to generate numerous and equiprobable realizations [54].

(a) Interpretation

Varying parameters were used to ensure the full potential and effectiveness of the method also to get the best possible facies model (Table 1). In total, three hundred realizations were generated, and the resulting facies model showed the distribution of all facies but on most of the layers, there was a lack in connectivity between the facies; the dimension of geobodies was not clear enough (Figure 11). There was a representation of each facies in the model, especially the lobe main body and channel facies were well presented. Each type of variable corresponds to an indicator variation function. The cross-sectional view gives more clear idea of the lateral and vertical distribution of facies in each zone also among $I$⁠, $J$⁠, and $K$ directions of the model. The lateral distributions of facies could be easily identified along $K$ direction, especially along $K=19$⁠, there was a representation of channels although not that perfect, also lobe main body, and lobe margin facies. The histograms reproduced from this method were accurate with respect to each facies, and the vertical proportion curves were also preserved during the upscaling process. The facies trends and variations were obvious along $I$ direction especially along $I=153$ and along $J=175$⁠. The resulting model also honored the fine-scale heterogeneities and facies distributions. For three realizations, separate variograms were used for each facie located within each zone. For another three realizations, the variogram settings remained the same for all facies. The facies type, major direction, minor direction, vertical direction, and azimuth for each realization are listed in Table 1.

Compared with the indicator kriging method, the sequential indicator simulation method could not restore the geometry of the channels in a more realistic manner, and the stacking pattern and distribution law of the sand body are not justified. Poor preservation of facies associations and geobody dimensions was encountered, especially connectivity among the facies associations was lacking. More specifically, the two-point-based method was unable to constrain the connectivity of (curvilinear) distributary channels.

(1) Multiple-Point Geostatistics (MPG). MPG is a pixel-based method that makes it more efficient than other Boolean/object-based methods in terms of conditioning to data and connectivity, especially in sparse data conditions. The training image generation and pattern creation are preprocesses for facies modeling through multiple-point geostatistics. The final output can be a pattern or tree explaining the neighbor relationship between facies and training images.

(a) Interpretation

The surface-based model was used to build a training image because of its simplicity, stationarity, perfect facies representation, and minimal processing time. Also, due to lack of good results using other approaches such as object modeling and interactive drawing, Chang 6_1_2 (surface) was used to first built a surface-based facies model to use as a training image, and the reason to choose this surface is that there is certain criteria to build better training images like simple (has to explain facies relationships), stationary (has to deliver good statistics), large enough to cover facies characteristics, and small enough to run within an acceptable time.

A new training image grid was used, i.e., a number of cells (⁠$I$⁠, $J$⁠, $K$⁠) $100∗100∗30$⁠. In most cases, two to six main facies are good enough; four facies were also enough to get better results as in this case. In search mask settings, the ellipsoidal shape was used instead of rectangular because the latter is bigger but slower. Several different parameters were utilized for the same training image to get the best MPG model. The best parameter settings for this model are shown in Table 2. The distributary channels were very clear and had most of the attributes as in the input training image, despite two channels, one originating from the west direction and one in the centered is incomplete possibly due to parameter settings or overtraining of image. Other facies such as lobe main body and lobe margin were identifiable, and this method worked well to model the facies as it is in the training image. Cross-sectional view provides a better idea of facies in all zones including vertical and lateral distributions and trends (Figure 12). The MPG method was capable of reproducing the facies associations; therefore, it had the capability to model three-dimensional facies continuity and connectivity. Similar to other methods like indicator simulation and truncated Gaussian simulation, MPG honors the fine-scale heterogeneities. This method also preserved the dimensions of geobodies especially distributary channels, lobe main body, and lobe margin (Figure 13). The accurate MPG model honors the reservoir heterogeneity that further enhances the interpretation and results of petrophysical parameters. The facies trends and variations were obvious but not accurately represented along $I$ and $J$ directions. The facies were more elongated and stretched; do not entirely follow the geometries of the facies especially along $I$ direction $I=153$ and along $J=175$⁠.

(1) Indicator Kriging (IK). Indicator kriging method offers a suitable interpolation approach for datasets (1) containing various observations below the limit of detection, (2) the strongly skewed histogram, or (3) attribute values specific classes are well connected in space. To apply the indicator kriging with its full potential needs, the multiple indicator semivariogram modeling and tedious inference are also the results postprocessing for attribute estimates retrieval and associated uncertainty measures.

(a) Interpretation

The shape and the development of the river phase were well modeled, but on various points, some facies associations were smoothed out, i.e., not well connected or with distortions. Especially along $J$ direction, modeling results were better as compared to $K$ direction (Figure 14). Along $K$ direction, the model results were better in all zones, and well-developed distributary channel facies were identified but channels were not developed in the lower part, in some layers in the upper and middle part with variation and channel facies supposed to be smoothed out on specific layers on different localities (Figure 15). Lobe main body facie was identifiable but with less abundance and configuration pattern. Lobe margin facie was present but with less abundance as compared to other facies.

In all directions, some layers showed better modeling results like along $K$ direction Index#19, 48, 72, 91, 129; along $J$ direction Index #78, 110, 180, 220, 265; and along $I$ direction Index #223, 268, 313, 324, 340. Along with $I$ and $J$ directions, the distributary channels shape, and contacts with other facies such as lobe main body and lobe margin facies were identifiable especially the key features of sand body configuration with the channels.

Overall, four hundred fifty realizations were generated with different variogram parameters settings (Table 3). The modeling results were suitable enough for interpretations of facies associations like mutual contacts, abundance, geometry, and extent along with different directions. For this technique, different approaches were adopted to get the best results, like calculation of variogram parameters using various methods, discrete data analysis for each facies, and each zone, using the same variogram for all facies and zones.

Better results were achieved using the same variogram for all facies and same zone settings for all zones. Realization results from parameter setting 6 showed much better results as compared to most of the other realizations. There was a more clear representation of distributary channels and lobe main body except for lobe margin facies along $K=72$ (Figure 15). The distribution and geometries of facies (distributary channels and lobe main body) were clear enough along $J=175$ (Figure 14) and along $I=153$⁠.

(2) Surface-Based Method. Surface-based modeling is an innovative geostatistical method that enables improved integration of geological information in clastic reservoir models. Traditionally, the geostatistical tools were limited to the model construction by geometric objects or by pixels. The surface-based method allows the reproduction of trend hierarchies and stacking patterns associated with the sedimentary processes for modeling by sequential fill of stratigraphic layers available for accommodation. Unlike other modeling techniques, surfaces were used to represent the stratigraphic zones and facies boundaries within a reservoir. Five surfaces have been used for building surface-based facies models.

(a) Interpretation

The surfaces used for this method include Chang $4+5_2_1$⁠, Chang $4+5_2_2$⁠, Chang 6_1_1, Chang 6_1_2, and Chang 6_1_3. The representation of each surface-based facies model was distinct for each zone, indicating all the facies dimensions, geometries, and trends. The facies configuration was identifiable along all directions; the spatial connectivity and stacking pattern of facies along all directions were well preserved (Figure 16). The complex geometries were honored including the shape of channels and lobe margins. Morphological features of all architecture elements were combined. Scale-independent stratigraphic surface geometries share a common geometric template. The shape and distribution of channels, lobe main body, and lobe margin were also identifiable along $J$-direction (Figure 17). Overall, the modeling results were very accurate and followed the principles of facies association stacking pattern and interfacies relationships. The surfaces provide a more precise approach as compared to the use of variograms as in some other methods. Each zone of the resulting models has been cross-checked to identify the unusual variation in facies behavior lateral and vertical extent of facies associations. After identifying the infrequent facies presence, interactive facies drawing tool was used to fix the issues.

According to the analysis of petrophysical properties characteristics in different facies, each of the facies has diverse proportions of petrophysical properties. Such as in distributary channel facies, the permeability ranges from 0.1 mD~ 1.5 mD (⁠$mean=0.72mD$⁠, $standarddeviation=0.48mD$⁠), and 90% values were between 0.01 mD and 1.2 mD. In the lobe main body facies, the permeability ranges from 0.01~ 1 mD (⁠$mean=0.38mD$⁠, $standarddeviation=0.35mD$⁠) with 88% values between 0.01 mD and 0.7 mD. In the lobe margin facies, the permeability ranges from 0.01 mD~ 0.4 mD (⁠$mean=0.13mD$⁠, $standarddeviation=0.12mD$⁠) with 95% values were between 0.01 mD and 0.2 mD. The porosity ranges from 2.4 ~ 10% in different facies. In distributary channel facies, the porosity range was 6.98% ~ 10% (⁠$mean=7.9%$⁠, $standarddeviation=1.4%$⁠), and 90% porosity values were between 6.9 and 7.4%. The range of porosity in lobe main body facie was 2.4% ~ 7.3% (⁠$mean=4.6%$⁠, $Std=1.5%$⁠), and 89% porosity values were between 2.4% ~ 5%. In the lobe margin facies, the range of porosity was 1% ~ 3% (⁠$mean=1.96%$⁠, $Std=0.67%$⁠), and 85% porosity values were between 1% ~ 2.5%.

Sequential Gaussian simulation (SGS) is a stochastic algorithm; this is performed in the spatial domain and is considered a very common method used for modeling petrophysical properties. If the reservoir property does not have a Gaussian distribution, then, the modeled property will go through a rank transformation to change it to a normal standardized form.

(1) Interpretation of SGS Porosity Model. The porosity model was built based on the corrected results of porosity logs. The SGS algorithm was used as a statistical method to construct a porosity model. Data analysis was first performed to calculate normal score transforms and experimental variogram ranges. Various variogram parameters were utilized in two hundred different realizations for porosity modeling (Table 4). The porosity modeling results from parameter number one were more accurate; the variogram ranges and azimuth were set according to the variogram calculation of each facies in data analysis. Figure 18 represents the average and frequency of porosity through Chang 6 member of Yanchang Formation.

The porosity in general in this unit was not very high that satisfies the porosity limit in the tight sandstone reservoirs. Porosity was increasing toward the eastern direction of the model. SGS reduced the smoothing effect as compared to the kriging. Small-scale heterogeneities were reproduced. Porosity models showed very low porosity values. The utilization of a variogram was important in porosity modeling. Two hundred eight well log data were analyzed, and porosity analysis was performed. 74% values were <2%, i.e., high frequency of porosity was <2%, 11.7% values were >9%, and 3%-4% average porosity in important wells (Figure 18). The porosity model showed a variation in porosity values but in most of the regions, the porosity was very low (Figure 19); the view from above also showed different ranges of porosity in the model and also cross-sectional provided better understanding of varying porosity. The $J$ direction also showed a better variation in the porosity model, particularly along $J=175$⁠. View from $I$ direction also was captured to show the variability in the porosity model focused on $I=153$⁠. The histogram also enhanced the understanding of the behavior of well logs, upscaled cells, and porosity variation in the model.

(2) Interpretation of SGS Permeability Model. The permeability model was built using a sequential Gaussian simulation (SGS) method. After facies modeling and scaling up well logs, data analysis was performed. While performing the data analysis process, the variogram parameters were calculated and applied on the real model. Different parameters were used for each facies (Table 5). Four realizations were generated to ensure the authenticity of the model and comparing the model to choose the best one. Logarithmic and normal score transformations were used. Got the best results from parameters settings number one (Table 5). The final 3D permeability model showed considerably low permeability values. Less smoothing effect was encountered using this method, and important wells and 3D permeability model are shown in Figure 20. After careful observations along $K$ direction (⁠$K=44$⁠), high permeability profiles were identified showing the high permeability zone in the model. We could get a better idea of the transformation of normal well logs, upscaled well logs, and permeability well logs (generated after modeling) shown in different colors in Figure 20. The lateral variation in permeability could be easily observed along $J$ and $K$ directions. Along $J$ direction (⁠$J=175$⁠), there was obvious representation of permeability indicated by color with legend and histogram. Most the values were around 0.01 to 1 mD, 74% values are 0.01-0.03 mD, ~18% values were >1 mD, and average permeability in important wells was 1-2 mD (Figure 21). Also, along $I$ direction (⁠$I=153$⁠), permeability profiles could be easily observed and identified that showed the heterogeneous behavior of permeability in a tight sandstone reservoir (Figure 22). The histogram provides better approach to understand the behavior of well logs, upscaled cells, and permeability variation in the model.

There were different geostatistical methods available for petrophysical modeling, but Gaussian random function simulation remained the best in terms of performance and random variable generation. The model was generated in Petrel based on effective physical property values and the Gaussian random function simulation algorithm. This simulation belongs to a conditional simulation algorithm that combines kriging and unconditional simulation [55]. Three different parameter sets have been used for modeling porosity and permeability for generating one hundred fifty realizations (Table 6).

(1) Interpretation of GRFS Porosity Model. The porosity model was populated by Gaussian random function simulation. Three realizations were generated using three different parameter settings for each facies. The average porosity values were in the range of 8%-10%. The porosity model with important wells locations was presented in Figure 23. The lowest porosity values were around 1% to 3%. Along $K$ direction (⁠$K=44$⁠) showed varying porosity values range from high medium to low. The $J$ and $I$ directions were also showing clear results to access the porosity trend in the model. Particularly $J$ direction $J=175$ showed high, medium, and low porosity zones. Along $I$ direction $I=153$ displayed various porosity zones (Figure 24). The normal, upscaled, and well logs after modeling were represented in a histogram.

(2) Interpretation of GRFS Permeability Model. Gaussian random function simulation was used for permeability modeling. Three realizations were built to further investigate the comparison among permeability models. The permeability model overall showed low permeability values (0.1 mD) (Figure 25). For a better understanding, the model was analyzed along different directions (⁠$I$⁠, $J$⁠, $K$⁠), and slide player was also used to check the permeability variation. While observations along $K$ direction the heterogeneous behavior of the model was observed (⁠$K=44$⁠). Some low and high permeability zones were also identified along $J$ direction (⁠$J=175$⁠). While along $I$ direction, the permeability was also in a range from high to low (⁠$I=153$⁠) (Figure 26). The permeability histogram also enhanced the understanding of well logs behavior before and after the modeling process. Three realizations showed variability in permeability profiles and behavior with the use of different parameter settings.

The inferences were drawn from the evaluation of facies and petrophysical modeling methods, and statistical results that were performed during this research:

• (1)

  Sequential indicator simulation (SIS) was proved to be more robust and honored the spatial distribution of facies, but the connectivity between facies lacks, and the geobody dimension was not clear. The resulting MPG model exhibits a clear distribution of the distributary channels. Lobe margin and lobe main body facies were identifiable and showed facies variation throughout the model. The indicator kriging (IK) modeling results were suitable for facies interpretation and facies associations. IK model was focused on local accuracy, but some facies were smoothed out

• (2)

  All the three facies modeling algorithms used during current research could be used for facies modeling of a tight sandstone reservoir. But after the evolution of the facies modeling results, the multiple-point geostatistics proved as the best method for facies modeling, despite some drawbacks, it still provides a better approach for facies modeling and fulfills most of the major requirements to build perfect facies model

• (3)

  SGS porosity model has represented reduced smoothing effect as compared to kriging and small-scale heterogeneities were preserved. The permeability model showed low permeability, less smoothing effect, and lateral and vertical variations. GRFS remains the best method for modeling petrophysical properties in terms of random variable generation and performance

• (4)

  SGS and GRFS both provided a better approach for petrophysical features but the GRFS method approach was more realistic with fewer drawbacks. Especially the use of unconditional simulation and cokriging makes it unique in addition to less computation time

The data is confidential and it cannot be provided on request.

The authors declare no conflict of interest.

This study was supported financially by the strategic cooperation technology projects of the China National Petroleum Corporation (CNPC) and the China University of Petroleum, Beijing (ZLZX2020-02) and the National Natural Science Foundation of China (no. 41872107). The authors also greatly appreciate the research institute of exploration and development of the PetroChina Changqing Oilfield Company for their support and permission to publish the results.