Introduction

The reservoir is formed under the combined influence of various factors such as structure, sedimentation and diagenesis. The uncertainty and stochasticity of the reservoir modeling are inevitable when it is constructed based on the limited data such as well logging and seismic data1,2,3. An important step in reservoir modeling is to build reservoir models that conform to the current cognition, minimize the uncertainty in modeling and evaluate the uncertainty of the model.

The uncertainty in the model can be categorized into response uncertainty, spatial uncertainty and local uncertainty4. The local uncertainty can be characterized by plotting the cumulative probability curve5. The spatial uncertainty can be described by using Euclidean distance or by using the correct point-to-percentage method in the scatter plot6,7. By constructing target models with different modeling methods and comparing the characteristics of standard deviation of simulation results under different conditions, the response uncertainty can be characterized8. Additionally, response uncertainty can also be quantified by variance analysis9,10.

The application of the processed seismic data to the model construction can effectively reduce the uncertainty of the model11,12. The adoption of new technical means to collect data can also minimize uncertainty13. Additionally, using higher-order modeling methods, such as multi-point geostatistics, can also effectively reduce the uncertainty14,15. Furthermore, using reservoir numerical simulation to verify and adjust the model, the uncertainty of the model can also be effectively reduced14,16.

Modeling by well pattern sparsing is a common evaluation method for model evaluating17, and it is also a desirable method to verify the model by means of three-dimensional probability distribution18. The introduction of experimental design methods, commonly used in medical research, has significantly improved the rationality of model verification and uncertainty evaluation19,20,21. Once the uncertain sources and evaluation factors in the modeling process are clearly defined, the uncertainty factors affecting the modeling results can be quantitatively analyzed through reasonable experiment design, and a suitable reservoir model can be selected22,23,24,25,26.

After identifying the main sources of uncertainty, this study quantitatively evaluates the uncertainty factors in the modeling process through the method of experimental design, ultimately optimizes the appropriate target model.

Geological setting

The study area is located in Xizhang Town, Kenli County, Dongying City, Shandong Province. It is structurally located in the south of Chenan fault and the north of Lijin-Binnan fault system. It is located in the first and second steps of the structural slope break belt, which is bordered to the west by the Binxian Uplift, to the north by the Chenjiazhuang Uplift, to the south by the Lijin oil-generating depression, and to the east by the Qingtuozi Uplift (Fig. 1).

Fig. 1
figure 1

Regional geological map of study area.

Full size image

The structure of the study area is relatively simple. To the west, the terrain is characterized by higher elevations in the north and lower elevations in the south, while to the east, it is high in the northeast and low in the southwest. The oil reservoir presents a trend of oil-bearing linkage from west to east, and the main oil-bearing sand body is the conglomerate complex fan in the near-shore underwater fan, located within the steep slope zone of the faulted lacustrine basin in the fourth member of Shahejie Formation. The internal lithology of the sand-gravel rock mass is complex, reflecting multi-stage rapid accumulation or re-collapse field deposition27,28. The electrical properties of the reservoir are greatly affected by the tight conglomerate. Core analysis showed that the porosity was concentrated in the range of 5.3-13.5% with an average of 7.8% (Fig. 2), and the permeability was 0.72 × 10 − 3µm2 to 39.07 × 10 − 3µm2 with an average of 4.2 × 10 − 3µm2 (Fig. 3). Due to the deep burial, the effective and non-effective reservoir are mixed, and the heterogeneity is serious. As a result, it is classified as an ultra-low porosity and ultra-low permeability reservoir.

According to the previous analysis of the source of ancient gully sediments in the sedimentary period of the upper submember of the fourth member of Shahejie Formation, there are three primary directions of source in the study area, which are the northern source, the lateral source in the west and the eastern source, which are typical multi-source deposition29,30,31.

Fig. 2
figure 2

The histogram of porosity distribution for the reservoir in study area.

Full size image
Fig. 3
figure 3

The histogram of permeability distribution for for the reservoir in study area.

Full size image

Materials and methods

Partition blocks

The seismic response characteristics across different zones of the study area are distinctly different. In the root area of the fan delta, the reservoir is thick with almost no shale, exhibiting broad and weak reflections, without distinct reflection axes or with reflection axes of large dip angles. In the middle area of the fan delta, the reservoir is relatively thick, with mudstone layers showing strong reflections and poor continuity. In the distal area of the fan delta, the reservoir is even thinner, with developed shale leading to local strong reflections and poor continuity(Fig. 4). In the pure shale areas, reflections are clear and continuous. Based on these observations, seismic facies patterns are established to provide a comprehensive basis for the subsequent boundary delineation of the fan body.

Fig. 4
figure 4

North-south seismic section of Y928C well.

Full size image

Based on the seismic division of blocks, more detailed boundary delineation can be achieved by analyzing the sand body thickness. In the study area, the effective reservoir is located in the fan middle facies zone, while the poor reservoir is in the fan root facies zone. The reservoir quality in the fan middle facies zone is significantly better than that in the fan root facies zone. Therefore, the primary focus is on the stratigraphic thickness of the fan middle zone. Analysis reveals that there are two thickness centers in the study area. In the western Yong935 area, influenced by a sediment source direction approximately from the north-northeast, the effective thickness is relatively thin in the 7th stage. The center of the fan middle zone’s stratigraphic thickness is located to the north of the Yong928C well, with a maximum thickness of 60 m. In the 8th stage, the center of the stratigraphic thickness shifts, roughly toward the northeast of the Yong928C well, with a maximum thickness of 240 m. In the eastern Yong936 area, influenced by a northward sediment source, the stratigraphic thickness is greater in the 7th stage, with the thickest center located near the Yong936 well, reaching a maximum thickness of 260 m. In the 8th stage, however, the stratigraphic thickness in the Yong936 area decreases significantly, shrinking to 140 m, and the center of the thickness shifts toward the northwest (Fig. 5).

Fig. 5
figure 5

Contour map of sand body thickness.

Full size image

Finally, integrating seismic insights and stratigraphic thickness, the boundaries of the fan bodies are delineated, dividing the study area into the eastern Yong 936 region, the central region, and the western Yong 935 region. Considering that defining a single boundary involves certain subjectivity and uncertainty, which cannot objectively represent the fan body boundary characteristics, the original boundaries were appropriately expanded and contracted. Three different boundary levels were established: pessimistic, optimistic, and most likely (Fig. 6). These boundaries serve as the foundation for subsequent construction of sand-mud facies models and attribute models.

Fig. 6
figure 6

Schematic diagram of partition boundary division.

Full size image

Collect sand-mud ratio

Within the study area, the lower limit of porosity is determined to be 4.7% and the lower limit of permeability is 0.7 mD using the forward and inverse cumulative method. The calculated lower limit of porosity is 5.3% and the permeability is 0.7 mD using the resistivity lower limit inversion method. Verification through core sampling and oil testing result in lower limit of porosity at 5.3% and a permeability of 1.6 mD. Selecting a single set of physical property limits as the criterion for evaluate sand-mud facies cannot effectively assess the impact of the effective property limits on the construction of sand-mud facies models and reserve calculations. Therefore, the physical property limits derived from the three methods are set as three different levels: Physical Property Lower Limit One, Physical Property Lower Limit Two, and Physical Property Lower Limit Three. Under the previously delineated three boundary conditions, different effective property limits are used to statistically determine the percentage of sand and mud within each sub-layer under different property limit conditions (Table 1). During the statistical process, when the percentage of mudstone facies reaches 100% under a given condition, the sand-mud ratio for that condition is set to 1.

Table 1 Statistics of sand-silt ratio under different conditions.
Full size table

Constructing conditional data

Constraint surface of azimuth

Due to the presence of multiple source directions in the study area, it is necessary to impose constraints on the azimuth of the variogram during the modeling process. Based on the defined boundaries of the western Yong935 area, the central area, and the eastern Yong936 area, and referencing the source analysis results, constraint surfaces for the source directions of the three regions were assigned. The sediment source direction in the Yong 935 area is 20° west of north, while in the Yong 936 area, it is oriented to the northeast. Therefore, the sediment source direction constraint for the Yong 935 area is set to 340°, and for the Yong 936 area, it is set to 45°. For the central region, the constraint is set to 0° (Fig. 7). For the three different levels of fan body boundary conditions, the same procedure is required to ensure that the azimuth data in the azimuth constraint files used for each boundary condition level are consistent.

Fig. 7
figure 7

Schematic diagram of constrained surface of azimuth.

Full size image

Constraint surface of lithology proportion

By processing the statistically obtained sand-mud rock percentage data under different conditions, the sand-mud rock percentage data for the Yong 935 area, central region, and Yong 936 area in the study area can be effectively utilized. Taking the 7 − 1 sub-layer as an example (Fig. 8), under the pessimistic boundary and effective property limit one conditions, the mudstone% in the western Yong 935 area is 76.3% and the sandstone% is 13.7%; in the central region, the mudstone% is 78.4% and the sandstone% is 21.6%; in the eastern Yong 936 area, the mudstone% is 64% and the sandstone% is 36%. Surface files are generated using the boundaries of the western Yong 935 area, central region, and eastern Yong 936 area, and the statistically obtained sandstone proportion values are assigned to the surface files of the three regions. New surface files are then created for the three regions, and the statistically obtained mudstone proportion values are assigned to the surface files of the three regions, ensuring that on the vertical direction, the sum of the sandstone constraint surface file values and the mudstone constraint surface file values for a given area is 1. Figures 3, 4, 5, 6 and 7 illustrates the sandstone facies constraint surface file for the 7 − 1 sub-layer under the pessimistic boundary and property limit one conditions. Since the sand-shale percentages vary under different fan boundary conditions using various effective physical property lower limits, it is necessary to construct sand-shale facies constraint files for each small layer under different constraint conditions. These surface constraint files serve as conditional data for building lithofacies models.

Fig. 8
figure 8

Constrained surface of lithology ratio (pessimistic boundary and materiality limit).

Full size image

Full factor experiment

In the preparation of basic condition data, the setting of boundaries affects the extent of the fan body, which in turn has a certain impact on the calculation of reserves. The physical property lower limits influence the statistical proportion of sand-mud rock within the boundaries, which in turn affects the construction of the lithofacies model and subsequent reserve calculations. Different variograms models affect the spatial relationships in the specific modeling, thereby impacting reserve calculations. Therefore, there are three main influencing factors during the experimental modeling process: boundaries, physical property lower limits, and variogram range. Based on the above analysis, a three-factor, three-level full-factorial experimental design can be developed as follows(Table 2):

Table 2 Header of full factor experiment.
Full size table

Results

Models with constrained data

By combining boundary conditions, physical property lower limits, and variogram range conditions, 27 distinct sand-mud facies models are constructed by Petrel201832. Consistent with the overall characteristics reflected by the sand-mud ratio statistical parameters, the sand-mud rock percentages vary under different boundary conditions, and the proportion of sand-mud facies in the models differs accordingly. Under the same boundary conditions and effective reservoir property lower limits, the characteristics of the sand-mud facies models show similarities. However, when any of these conditions change, the characteristics of the constructed sand-mud facies models will change significantly.

Using the sand-mud ratio file obtained under optimistic boundary conditions and effective property lower limit one as constraints, the sand-mud facies models built with different variogram range conditions are analyzed (Fig. 9). Under these conditions, it was observed that the sandstone facies characteristics were not prominent in the 7th sub-layer of the western Yong 935 area, but were distinct and continuous in the 7th sub-layer of the eastern Yong 936 area. The sandstone facies characteristics under the three different range conditions were similar, but with minor differences in local areas. A comparison using histograms shows that the percentage of sand-mud facies fit under the three different variogram conditions differs only slightly. Under the pessimistic variogram condition, the sandstone facies accounted for 15.0% of the model, under the possible variogram condition, the sandstone facies accounted for 15.4%, and under the optimistic variogram condition, the sandstone facies accounted for 15.5%. The change in the percentage of sand-mud facies may be that as the range of the variogram increases, the variogram’s representation of the spatial structural characteristics of the sand-mud facies changes, leading to subtle differences in the final model.

Fig. 9
figure 9

Lithofacies model under conditions of optimistic boundary and lower limit32.

Full size image

Models without constrained data

In contrast, when local variably constrained variograms are not applied, the characteristics of the sand-mud facies model in the study area change (Fig. 10). Spatially, the Yong 935 block in the west and the Yong 936 block in the east remain the primary distribution areas for sandstone facies, but the distribution of sandstone facies becomes more discrete. Cross-sectional observations reveal a significant abrupt change in the thickness of sandstone facies at the boundaries of the zones. Comparisons indicate that the lithofacies models constructed under the constraints of variograms with local variations are more consistent with geological understanding.

Fig. 10
figure 10

Lithofacies model without constraints of local variogram32.

Full size image

Property model

After completing the data variation and variogram analysis for the physical property data, a facise-controlled modeling approach was adopted, using the facies model as the control constraint. The sequential Gaussian simulation method was used to construct the porosity model under different constraint conditions (Fig. 11a), while the permeability models(Fig. 11b) and water saturation(Fig. 11c) models were constructed under the constraint of the porosity model. All kinds of property models are constructed using Petrel201832.

Since the study area has two source directions, the direction was constrained by generating constraint surfaces. The preparation of the relevant constraint files has been discussed in the section on condition data processing. To enhance the reliability of subsequent reserve calculations, five corresponding physical property models were generated under the same facies model constraint.

The NTG model (Fig. 12) is obtained by calculating using the lower limit of Physical property lower limit. The calculation formula under the condition of lower limit 1 is formula1, under the condition of lower limit 2 is formula 2, under the condition of lower limit 3 is formula 3:

$$\:NTG1=IF(\text{P}\text{O}\text{R}>=4.7\:\text{A}\text{N}\text{D}\:\text{P}\text{E}\text{R}\text{M}\:>=\text{0.7,1},0))$$
(1)
$$\:NTG2=IF(\text{P}\text{O}\text{R}\ge\:5.3\:\text{A}\text{N}\text{D}\:\text{P}\text{E}\text{R}\text{M}\:\ge\:\text{0.7,1},0))$$
(2)
$$\:NTG3=IF(\text{P}\text{O}\text{R}>=5.3\:\text{A}\text{N}\text{D}\:\text{P}\text{E}\text{R}\text{M}\:>=\text{1.6,1},0))$$
(3)
Fig. 11
figure 11

Propety model under conditions of optimistic boundary and lower limit32:(a) porosity model and (b) permeability model and (c) water saturation model.

Full size image
Fig. 12
figure 12

NTG model under conditions of optimistic boundary and lower limit32.

Full size image

Reserve calculation

The fundamental idea of applying a reservoir model to calculate reserves is to use the three-dimensional grid, defined during the modeling process, as the basic unit. Reserves are then calculated using Eq. 4, based on specific information within each grid, such as porosity, permeability, and water saturation.

$$\:STOIIP=\varSigma\:\left[Vbulk*NTG*POR*\right(1-SW)/Bo]\:$$
(4)

Vbulk: Bulk Volume; NTG: Net to Gross Ratio; POR: Porosity;

SW: Water Saturation; Bo: Oil Formation Volume Factor。.

Under each facies model constraint, five corresponding physical property models were implemented, resulting in 135 reserve calculation values(Table 3). The final evaluation metric is the average reserve value for each facies model condition.

Table 3 Processing table of experimental result.
Full size table

Uncertainty analysis

Full-factor experiment analysis

In the main effect diagram of the full-factor experiment (Fig. 13), under the pessimistic boundary conditions, the reserve value is the largest; under the possible boundary conditions, the reserve value is the smallest; and under the optimistic boundary conditions, the reserve value is in the middle. The reserves under the three boundary conditions do not exhibit a monotonically increasing or decreasing trend with the expansion of the boundaries, but rather show a pattern of decreasing first and then increasing as the boundaries expand. The reason for this characteristic is that as the boundaries increase, the area of the fan body also increases, and the corresponding statistical values of sand-mud percentages change accordingly. Under the possible boundary conditions, the percentage of sandstone reaches its minimum, leading to the smallest proportion of sandstone facies in the model, and the characteristics of the analysis results are consistent with the results of the lithofacies modeling.

Fig. 13
figure 13

Main effect plot of full factor experiment.

Full size image

Similarly, under the conditions of effective reservoir properties, the reserves are the smallest under pessimistic physical property lower limits and the largest under optimistic physical property lower limits. As the effective reservoir physical property lower limits change from optimistic to pessimistic, the reserve values also change from the largest to the smallest, showing a positive correlation between the trend of reserve values and the trend of effective reservoir physical property lower limits. The reason for this characteristic is that as the effective reservoir properties shift from an optimistic level to a pessimistic level, the lower limit values of effective reservoir properties increase, leading to a decrease in the percentage of sandstone, ultimately reducing the proportion of sandstone facies in the model, and thus decreasing the reserve calculation values.

In the range of variogram, when the range is at its minimum, the reserve value is also at its minimum, and when the range is at its maximum, the reserve value is at its maximum. The trend of reserve values shows a positive correlation with the increase in range. In the specific process of lithofacies modeling, as the range of the variogram increases, the spatial distribution range and connectivity of sandstone show a certain trend of improvement, which also well explains the characteristic in the analysis results that reserves increase with the increase in range.

Table 4 Variance analysis of reserve factor.
Full size table

The R-Sq statistics for the model based on analysis of variance, the predicted model, and the adjusted model are 98.82%, 86.61%, and 96.18% (Table 4). Respectively, the significance analysis of the factor analysis, shown in Fig. 14, indicates that boundary conditions and lower limit conditions are significant influencing factors. This is consistent with the information presented in the factor main effect plots and interaction plots above.

Fig. 14
figure 14

Pareto chart of full factor experiment.

Full size image

Response surface analysis

The full-factor experiment clarified that the boundary range and property lower limits have a significant impact on the mean value of reserve calculations. Therefore, in the response surface analysis, boundary range and property lower limits are selected as experimental factors. The appropriate variogram range is then adjusted to minimize its impact on reserve calculations, ensuring the stability of the calculation process. By modifying the boundary conditions and property lower limits, the optimal parameter combination is ultimately determined.

The experiment follows the principles of Central Composite Design (CCD), with the possible boundary and physical property lower limit two (POR > 5.3, PERM > 0.7) serving as the center points for response surface analysis. The levels of the two significant influencing factors are shown in Table 5. According to the principles of Central Composite Design, the experiment consists of 13 runs, with each run involving the construction of a facies model and three equiprobable physical property models, followed by the calculation of three reserve results. The final evaluation is based on the average value, and the experimental results are provided in Table 6.

Table 5 Level and code of experiment variables.
Full size table
Table 6 Design scheme of response surface.
Full size table

Probabilistic reservoir simulation

The Monte Carlo simulation method was used for 1,000 reserve calculations, replacing the actual modeling and reserve calculation operations. By probability distribution graph (Fig. 15), the pessimistic reserve P10 is found to be 954.5 × 104 tons, the probable reserve P50 is 1059.7 × 10^4 tons, and the optimistic reserve P90 is 1164.9 × 104 tons.

Fig. 15
figure 15

Cumulative probability distribution curves of reserves.

Full size image

Among the 1000 Monte Carlo simulation calculations, the reserve values close to the P10, P50, and P90 probability reserves are identified, and their corresponding variable values are converted to actual values. It is found that the influencing factor values for the P10 probability reserve are similar to those under the possible boundary and effective reservoir property limit three conditions. The influencing factor values for the P50 probability reserve are similar to those under the optimistic boundary and effective reservoir property limit two conditions. The influencing factor values for the P90 probability reserve are similar to those under the optimistic boundary and effective reservoir property limit one conditions. Models are constructed under these three sets of conditions using Petrel software, and reserves are calculated. Comparing the actual reserve calculation values with the Monte Carlo reserve prediction values, the errors are all less than 2%.

By performing a cross-analysis of the 13 calculated reserves with the 3P reserves (Fig. 16), a regression equation is obtained: Y = 1.042X − 42.71 (X is the model reserve and Y is the reserve calculated by the Monte Carlo or regression formula), with an R2 value of 0.986. This indicates that the multiple regression equation derived from the response surface experiment analysis has a certain degree of predictive ability. Additionally, combined with the Monte Carlo simulation method, it is possible to predict reserves under any conditions within the response surface.

Fig. 16
figure 16

Cross plot of model-predicted versus proxy equation-calculated GIIPs.

Full size image

Conclusions

This study takes the Yong 935–936 block of the Yanjia Oilfield in the Dongying Depression as an example. Aiming at the significant uncertainties present in the construction of a single geological reservoir model, three influencing factors affecting reservoir modeling and reserve calculation uncertainty were selected: fan body boundaries, effective reservoir property thresholds, and variogram range. A full-factor experimental design method was used to design 27 different modeling schemes, ultimately achieving multiple equally probable sand-mud facies geological models. With the static reserve as the goal, the study qualitatively and quantitatively evaluates the significant level of influence of each influencing factor. Through response surface experimental analysis, parameters are optimized, and the best choices for each influencing factor are determined, ultimately constructing the target model for the study area. The main insights gained from this study are as follows:

  1. (1)

    In response to the reservoir characteristics of the study area with multiple sediment sources, a locally varying variogram was used for unified simulation of the study area, effectively solving the problem of abrupt changes in sand body thickness at partition boundaries. Seismic attributes and sand body thickness maps were comprehensively utilized to delineate pessimistic, possible, and optimistic boundary levels. Under different boundary conditions, the percentages of sand-mud under different effective reservoir property thresholds were statistically analyzed, and the obtained data were applied as conditional constraints in the reservoir modeling process, effectively reducing the uncertainties in reservoir modeling.

  2. (2)

    Combining the characteristics of the study area and the selection of variograms during the modeling process, the fan body boundary, effective property threshold, and variogram range were analyzed as the main influencing factors in reservoir modeling and reserve calculation. For different influencing factors, optimistic, possible, and pessimistic levels were set, and a full-factor experimental analysis was conducted. The significance and causes of different factors were analyzed, and it was determined that the significant influencing factors for reserve calculation were the effective property threshold and fan body boundary, with the effect of the effective property threshold being more pronounced. The non-significant influencing factor was the variogram range. In response to the actual situation of modeling and the level setting of significant influencing factors, a response surface analysis method was used to analyze the significant influencing factors and reserve results, obtaining a multivariate regression model between reserves and significant influencing factors. A Monte Carlo random simulation method was employed to obtain the cumulative probability reserve distribution, predicting the 3P reserves of the study area.

  3. (3)

    All kinds of assumptions used in this modeling process, such as boundary conditions and physical property limits, are dependent on the special geological and geophysical data of Yong935-936 block. Therefore, when applying this modeling method to other types of reservoirs, specific analysis and judgment of the corresponding reservoirs are required. Otherwise, it will lead to unnecessary bias, which will affect the result of the final reservoir modeling.