Semantic rule-guided three-dimensional reservoir modeling method using an improved multiple-point geostatistics simulation

Purpose: To overcome the reliance on manual parameter adjustment and the difficulty in uncertainty quantification of traditional Multiple-Point Geostatistics (MPS) methods, and to address the general lack of geological semantic rule utilization that limits the reliability of existing 3D reservoir models.

Methods: A 3D reservoir modeling method integrating semantic rule guidance and MPS simulation is proposed. Geological semantic rules for sequence stratigraphy and fault kinematics were incorporated into the structural modeling workflow. Subsequently, an MPS algorithm featuring automatic domain kernel function selection and a rapid dual-dimensional parameter optimization strategy was employed for lithofacies and property simulation.

Results: Experimental results show that the models constructed by the proposed method not only maintain geological semantic consistency and coherence but also accurately characterize the spatial distribution of lithofacies and properties in the study area.

Conclusion: This method effectively integrates semantic rules with an optimized MPS workflow, enhancing the geological reliability and representational accuracy of 3D reservoir models for improved reservoir characterization.

1 Introduction

Three-dimensional (3D) modeling is a key technology in reservoir characterization and hydrocarbon development (Liu et al., 2022; Nemes, 2023). This technology enables precise characterization of the spatial structures and properties of subsurface reservoirs. Its quality directly determines the scientific validity and accuracy of reservoir management, reserve assessment, and production planning (Cui et al., 2025). Based on the 3D digital model, researchers can systematically analyze the geometric morphology, scale characteristics, and spatial distribution patterns of hydrocarbon reservoirs, enhancing the scientific basis for decision-making (Wu et al., 2006; Grose et al., 2021).

Multiple-Point Geostatistics (MPS) is a pivotal technique in the field of 3D reservoir characterization (Zhang et al., 2010). This method utilizes high-dimensional pattern matching, to capture multiple-scale spatial features from training images as multiple-point statistics. Compared to traditional geostatistical methods, MPS demonstrates significant advantages in characterizing lithology and lithofacies associations (Straubhaar et al., 2011). Its theoretical foundation is inherently well-suited for geological scenarios with strong heterogeneity, such as fluvial depositional systems and fault intersection zones (Feng et al., 2022). However, the essence of MPS remains rooted in a data-driven statistical modeling framework (Wambeke and Benndorf, 2016). Its simulation process relies on the spatial statistical distribution of numerical data, making it difficult to integrate fundamental principles of sedimentology and stratigraphy, such as stratigraphic connectivity and fault distribution. When modeling complex faults or heterogeneous depositional systems, relying solely on statistical features may lead to the neglect of critical geological constraints, thereby compromising the scientific validity and reliability of reservoir prediction (Liu et al., 2021).

A non-iterative MPS algorithm ENESIM was first proposed by Guardiano and Srivastava (1993), which laid the foundation for subsequent research in MPS. Strebelle (2002) presented the first available MPS algorithm SNESIM, which employs a search tree to store all possible conditional probability distribution functions in a single operation. This approach significantly enhances the efficiency of MPS simulations. The SIMPAT algorithm, Arpat (2004) further advanced the development of MPS, although it remains applicable only to the simulation of discrete variables. The FILTERSIM algorithm was proposed to advance MPS simulation by using linear filters to score training images and perform pattern classification Zhang, (2006). Mariethoz et al. (2010) proposed the Direct Sampling (DS) algorithm, which significantly reduces memory costs and improves simulation efficiency by directly sampling multiple-point patterns from training images, thereby bypassing the estimation of conditional probability distributions (Oriani et al., 2014; Li et al., 2024). Maxelon et al. (2009) introduced a three-dimensional reservoir modeling workflow for structurally complex areas. This workflow ensures the internal consistency of the generated 3D model through the interactive editing of geometries in both cross-section and map views. However, MPS methods face challenges in parameterization, as the selection of parameters directly impacts the quality of the simulation results.

To overcome the above-mentioned limitation, Gravey and Mariethoz (2020) proposed the adaptively parameterized Quick Sampling (QS) algorithm, Fast Fourier Transform (FFT) was used to decompose the distance norms and compute the compatibility mapping between search patterns and each position in training images to identify optimal candidates. The method utilized a partial sorting algorithm to rank these candidates, thereby significantly enhancing computational efficiency and operability. Gravey and Mariethoz, (2022) further proposed an automated parameter calibration method for QS, by using training images as the only reference. This approach identifies optimal parameters for pixel-based MPS simulation where the parameters should not remain static but require dynamic adjustment throughout the simulation process. An automated parameter optimization method for MPS was developed by Baninajar et al. (2019) to evaluate pattern reproduction via cross-validation with training images. The optimal MPS parameters are determined by using the Simultaneous Perturbation Stochastic Approximation (SPSA) in this method.

In the field of MPS-based 3D reservoir modeling, Chen et al. (2018) proposed a MPS-based 3D geological modeling method by using 2D cross-sections to reconstruct 3D geological models. This approach uses 2D cross-sections rather than 3D prior models as training images, enabling direct extraction of spatial features from 2D training images. Hou et al. (2021) leveraged the spatial features derived from 2D training images and implemented a globally optimized, MPS approach based on a hierarchical strategy to reconstruct 3D geological structures. Wang et al. (2022) proposed the Se3DRCS algorithm based on 2D training images and seismic information. In this method, an initial model was iteratively updated until the relative error between the observed values and synthetic seismograms falls below a threshold, thereby obtaining the optimal geological model. While most of these methods primarily focus on adaptive optimization of MPS algorithms and the integration of multivariate data with training images, they exhibit a notable lack of embedding geological semantic and knowledge.

The key to constructing a high-precision reservoir model lies in integrating multiple-source data and semantic rules to achieve 3D reconstruction of reservoir facies, and properties (Kathuria et al., 2019). However, there is the challenge of achieving collaborative integration between geological semantic rules and multiple-source data. In addition, it is difficult to reproduce heterogeneous facies and properties that adhere to structural constraints. To address the above-mentioned challenges, we propose a geological semantic rule-driven multiple-point geostatistical reservoir modeling method. First, spatial patterns of lithofacies, and properties are extracted based on a structural constraint. Meanwhile, an automated parameter optimization algorithm for multiple-point geostatistical stochastic simulation is proposed to enhance the accuracy and reliability of 3D reservoir characterization.

2 Methodology

In this work, we proposed a 3D reservoir modeling method based on semantic rules and MPS, which integrates the dual advantages of knowledge-driven and data-driven approaches. This approach leverages the unique strengths of MPS in characterizing the spatial heterogeneity of complex geological bodies (Okabe and Blunt, 2007). The proposed method reduces geological uncertainty through the systematic application of semantic rules and geological constraint techniques, thereby enhancing the accuracy of the reservoir model. The overall workflow of the proposed method is shown in Figure 1.

Figure 1

Figure 1. Workflow of the 3D reservoir modeling method guided by semantic rules and MPS.

2.1 Reservoir structural modeling guided by semantic rules

In the field of 3D reservoir modeling, structural modeling is a key step in building the geological framework. Geological semantic rules are a rule system that translates geologists’ expertise, principles of geological processes, and observational data into quantifiable and computable constraints for modeling. The integration of geological rules from multi-source data is essential for producing a final model that is both data-constrained and geologically consistent (Frank et al., 2007; Foged et al., 2014). Sequence stratigraphy and fault kinematics are exemplary representatives of geological semantic rules. This study facilitates the paradigm shift in geological modeling from “data-driven” to “knowledge-data dual-driven”.

2.1.1 Sequence stratigraphy

Sequence stratigraphy is grounded in the systematic analysis of stratigraphic sequences, their thickness, and attitude. It enables precise characterization of stratigraphic contact relationships and vertical stacking patterns. Stratigraphic sequences reflect the phased characteristics and formation mechanisms of sediment accumulation. Stratigraphic thickness is critical for constructing accurate models of vertical stratigraphic structure. Stratigraphic attitude, which describes the spatial geometric properties of strata (such as dip angle, strike), enhances the depiction of contact relationships and vertical stacking patterns. In the geological modeling process, sequence stratigraphy supports the creation of reservoir models by extracting stratigraphic evolution patterns. By integrating sedimentology theories, it significantly enhances the geological realism of the models.

In implicit structural modeling, the distribution of strata is determined by the values of a scalar field (Renaudeau et al., 2019; Zhang et al., 2023). The top stratum surface is set as the reference, with its corresponding isosurface value at 0. The scalar values for the isosurfaces of other strata are defined based on either their thickness or their distance from this reference stratum surface. To incorporate stratigraphic sequencing into the modeling process, the youngest stratum is added first. The FDI algorithm is then applied to interpolate this stratum, using the geometry of the most recent event to constrain older geological events. Through progressive interpolation (Mallet J. L., 1992; Mallet J. L., 1992), a complete stratigraphic model is gradually developed, ensuring alignment with sedimentary evolution laws and the principles of sequence stratigraphy. The workflow for implementing implicit structural modeling method based on sequence stratigraphy principles is shown in Figure 2.

Figure 2

Figure 2. Implicit structural modeling based on sequence stratigraphy principles.

2.1.2 Fault kinematics

Typically, it is infeasible to directly measure the slip direction of faults, and it is necessary to define it using probabilistic models based on geological knowledge. The correct fault information can be captured by incorporating additional kinematic constraints into implicit modeling algorithms (Zhong and Lin, 2019; Irakarama et al., 2021). Fault kinematics data effectively constrains fault morphology and surrounding deformation zones in structural modeling by defining the geometric morphology of faults, displacement characteristics, and their relationships with stratum evolution. This ensures that the structural model closely aligns with actual structural distribution. The specific workflow for incorporating fault kinematics data into implicit structural modeling is shown in Figure 3.

Figure 3

Figure 3. Fault kinematics-constrained implicit structural modeling.

As shown in Figure 4, the spatial distribution of fault kinematic effects can be represented by an elliptical pattern. The major axis of the ellipse is oriented parallel to the fault strike, representing the along-strike propagation of the kinematic influence, while the minor axis, perpendicular to the strike, defines its transverse extent. The length of the major axis corresponds to the projected length of the fault plane. The length of the minor axis is a definable parameter that can be calibrated based on specific geological conditions or modeling objectives. The ellipse is typically centered on the midpoint of the fault trace, which provides a well-defined and practical boundary for the fault-affected zone in the model. This geometric simplification offers a tractable framework for integrating fault kinematics into subsequent reservoir modeling workflows (Maxelon et al., 2009; Klise et al., 2009).

Figure 4

Figure 4. Schematic diagram of fault-affected zone.

The structural model is crucial for ensuring the realism of the 3D reservoir model. Within the reservoir modeling workflow, stratigraphic sequences are used to enhance the characterization of stratigraphic contact relationships and vertical stacking patterns, while fault kinematics data are incorporated to constrain the geometric morphology and the influence range of faults. By applying the proposed structural modeling method that integrates multi-source data and geological semantic rules, an accurate reservoir structural model can be established. This reliable structural model then serves as a foundation for subsequent reservoir property characterization.

2.2 Automatic parameters optimization of QS algorithm

The automated parameter group optimization method builds upon the QS algorithm by integrating an automatic neighborhood kernel function selection strategy and a rapid dual-dimensional parameter group optimization strategy. This enables dynamic adjustment of kernel functions during the simulation process and facilitates the rapid identification of globally optimal parameters.

The selection of kernel functions directly influences the spatial structure reproduction capability, computational efficiency, and robustness of the simulation results. By weighting the matching errors of neighboring pixels, the kernel function determines the contribution of pixels at different distances to the current simulation node. Its essence lies in balancing the ability to capture local details and global patterns. The automatic selection process of domain kernel functions is illustrated in Algorithm 1.

Algorithm 1

Algorithm 1.

As shown in Figure 5, the automatic kernel selection strategy is driven by an analysis of local structural heterogeneity in the training image (Tan et al., 2014), which dynamically determines the optimal kernel for the simulation node. This effectively overcomes potential simulation biases caused by fixed kernel functions, thereby improving the accuracy of the reproduced local features.

Figure 5

Figure 5. Workflow of the automatic parameter optimization of QS algorithm.

Simultaneously, the dual-dimensional parameter rapid optimization strategy is used to identify the optimal simulation parameter combination. The QS algorithm generates random fields with complex spatial structures by matching patterns from a training image. The neighbor count $n$and candidate count $k$are its core parameters. Their selection directly determines the accuracy, diversity, and computational efficiency of the simulation results. The parameters $n$and $k$are complementary. The parameter $n$governs the complexity of local patterns and the computational load. In contrast, $k$regulates randomness and prevents overfitting by controlling the size of the candidate pool. Together, they balance pattern-matching accuracy, computational efficiency, and result diversity. The implementation workflow for the proposed rapid optimization strategy for the dual-dimensional parameter set is presented in Algorithm 2.

Algorithm 2

Algorithm 2.

This strategy establishes an $n-k$ 2D parameter grid search framework to traverse and evaluate different parameter combinations, where $n$ is the neighborhood size and $k$ is the candidate count. Thereby identifying the global optimum parameter combinations within a broader parameter space. Compared to traditional manual tuning methods, the proposed strategy enhances the global reliability of the obtained optimal parameter combination.

By adjusting parameter n, the algorithm selects a neighborhood range more appropriate for the geological context, thereby ensuring the representativeness of captured patterns. Simultaneously, adjusting parameter k retains stochasticity during pattern selection, which balances the diversity and stability of the simulation results. The automatic optimization of kernel functions further enhances the accuracy of local features. Collectively, these proposed parameter optimization methods enable the algorithm to identify optimal parameters, ensuring the model aligns more closely with the actual reservoir.

2.3 Reservoir modeling constrained by geological structure

The key to achieving high-precision 3D reservoir modeling lies in establishing a workflow capable of effectively integrating quantitative data with qualitative geological knowledge. A structural model incorporating multi-source data and geological semantic rules is used as the backbone to constrain subsequent simulations of reservoir lithofacies and petrophysical properties (Schaaf et al., 2021).

To bridge the gap between structural modeling and property simulation, well data upscaling is required. A significant resolution mismatch exists between raw well logs and the model grid. Therefore, detailed log data are upscaled along the wellbore trajectory using averaging methods to assign equivalent properties to each intersected 3D cell. This process provides grid-scale-matched hard data as input for the subsequent QS algorithm. Under the dual constraints of the structural model and these hard data, the improved QS algorithm is employed for the 3D simulation of lithofacies and petrophysical properties. The structural model influences the simulation in two primary ways, as explicit boundary conditions, it defines structural limits such as fault surfaces and formation tops, as stratigraphic trend constraints, it encodes layer thickness variations and depositional dip. The structural model provides spatial constraints in the form of structural boundaries and stratigraphic trends, under which the QS algorithm matches multiple-point patterns. Thus, by integrating structural constraints with MPS, the generated 3D reservoir model not only honors the hard data but also aligns with macroscopic geological principles. The corresponding modeling workflow is shown in Figure 6.

Figure 6

Figure 6. Modeling workflow under structural model constraints.

3 Experimental results

3.1 Research area

The northern part of the Songliao Basin, located in northeastern China, boasts abundant shale oil resources. To validate the effectiveness and practicality of the proposed 3D reservoir modeling method, a specific area within this region was selected as the research focus.

The first processing step involves the depth conversion of the seismic horizon data, which yielded the depth-domain horizon data for Q1 through Q9 and T2. In this work, the prefix “Q” refers to the Quaternary strata, with numerals 1 through 9 assigned in ascending chronological order (thus, Q1 is the most ancient and Q9 the most recent). Conversely, “T” designates the Triassic strata, and T2 specifically identifies the Middle Triassic strata. To enhance the constraints on faults, the preprocessing step retained the data points and applied normal vector constraints in regions where the seismic horizon data has a high gradient. For relatively flat areas, data filtering was performed to ensure the quality and accuracy of the input data. The scatter plots of the processed seismic horizons and the normal vector constraints are shown in Figure 7.

Figure 7

Figure 7. Seismic interpretation horizon data. (a) Seismic interpretation horizon data. (b) Seismic interpretation horizon data and normal vector.

Furthermore, the effective integration of multi-source and heterogeneous data is a critical step in the data-driven structural modeling workflow. As shown in Figure 8, common multi-source data also includes well log data, boundary data and fault data, among others. The well log data used in this study were acquired from ten horizontal wells, designated as ${QH}_0$ through ${QH}_9$. Boundary data serves as an effective means to constrain the range of constructed model.

Figure 8

Figure 8. Multi-source data in the structural modeling process. (a) Well log data. (b) Boundary data. (c) Fault data.

3.2 Experimental results

3.2.1 Preprocessing of multi-source data

In implicit structural modeling, each geological horizon has its corresponding scalar value. The uppermost horizon is designated as the initial horizon, and its scalar field value is set to zero. The scalar values of the remaining horizons are assigned based on the vertical distance to the initial horizon that is calculated from their respective depth. To ensure the rationality and accuracy of the model, a correct stratigraphic sequence must be established. Before modeling, all horizons interpreted from seismic data are arranged according to their true geological order. The complete geological model is then progressively constructed following the geological timeline. Leveraging finite difference methods and a Cartesian grid, the FDI interpolation algorithm efficiently discretizes continuous geological models into regular volumetric cells. The regularization operator within the algorithm suppresses geologically inconsistent artifacts and accurately represents fault structures and other complex geological scenarios. As shown in Figure 9, this study adopts a top-down geometric modeling approach. Starting from the youngest top horizon, and integrating seismically interpreted horizons, boundary data, and sequence-stratigraphic constraints, the initial geological surface model was constructed using the FDI algorithm.

Figure 9

Figure 9. 3D structural surface model using FDI interpolation and a Top-Down modeling approach.

Faults, as critical structural features within or between rock layers, not only serve as key indicators of geomechanically evolution but also directly influence reservoir connectivity and hydrocarbon migration pathways. In three-dimensional structural modeling, fault data are typically incorporated in the form of fault sticks, which provide essential spatial location and orientation constraints for the model. By extracting geometric attributes such as the normal vector of the fault surface, a displacement field can be quantitatively constructed during the modeling process, ensuring a reasonable strain distribution between the two fault blocks during slip. The fault data used in this study are shown in Figure 10.

Figure 10

Figure 10. Fault data in different geometric representations.

By incorporating fault kinematics data (e.g., orientation, center, slip vector, and displacement) into structural model, the geometric fidelity and practical value of the model are significantly enhanced. The spatial extent and shape are defined by the major and minor axes, and the center point acts as a geometric reference for location of the fault in the reservoir model. The slip vector describes the relative motion between the two sides of the fault, and the displacement magnitude quantifies the amplitude of this movement. The 3D structural surface model containing faults is shown in Figure 11. It can be observed that the faults have not only influenced the dip angle and thickness of the local strata but have also caused deformation of the strata near the faults. Due to their strong control over reservoir connectivity and permeability, structurally complex fault zones are a key focus in hydrocarbon accumulation studies.

Figure 11

Figure 11. 3D structural surface model with faults. (a) Top view. (b) Side view.

The geometric framework and basic spatial structure of the model are defined by the distribution characteristics of scalar values, gradient variation patterns, and spatial correlations within the scalar field. As shown in Figure 12, the complete scalar field model in the study area exhibits distinct layering. When scalar values are mapped to a continuous color spectrum, the differences in stratigraphic depths and corresponding geological structures are clearly revealed. Specifically, shallow formations are characterized by lower scalar values, while deeper formations correspond to higher values. Local anomalies, such as abrupt value changes or anomalous gradient fluctuations, typically indicate the presence of faults and their associated damage zones.

Figure 12

Figure 12. Complete scalar field model of the research area. (a) Top view. (b) Side view.

3.2.2 Modeling workflow integrating multi-source data and geological semantic rules

We used the proposed method to build a reservoir structural model for a specific area in the Songliao Basin, based on the multi-source data and geological semantic rules. Sequence stratigraphy and fault kinematics, data were used as constraint in the modeling process. As shown in Figure 13, the detailed 3D reservoir structure model was generated using FDI interpolation algorithm. The model shows that the faults accurately follow the spatial distribution of the input data, with their dip directions and angles appearing natural and consistent with the regional structural style. Across all depth intervals, the model maintains proper spatial correspondence between the strata and faults. The sequence of horizons strictly adheres to the sequence-stratigraphic constraints, and the shallow and deep structures are vertically coherent, realistically representing the stratigraphic architecture of the study area. In summary, the 3D geological structure model resulting from the integration of multi-source data and geological semantic rules exhibits a high degree of accuracy and consistency.

Figure 13

Figure 13. Refined representation of the 3D structural model. (a) 3D structural model. (b) Top-down view. (c) Due north. (d) Due west. (e) Due south. (f) Due east.

To achieve seamless integration between structural modeling and property simulation, and to employ structural constraints in the property modeling process, we first perform upscaling on data such as well logs. Log data upscaling is the process of filling well log properties into the regular 3D grid. This resulting grid with property values is the basis for subsequent 3D reservoir simulation. Based on data from ten horizontal wells, we performed upscaling for porosity and permeability, with the results shown in Figure 14.

Figure 14

Figure 14. Well log upscaling. (a) Porosity data well upscaling. (b) Permeability data well upscaling.

The QS algorithm was used to simulate the 3D property model based on the upscaled well log data. The corresponding generated porosity and permeability models are shown in Figure 15. Figure 15a shows a plan view and section views looking north and east of the porosity model, while Figure 15b presents a plan view and section views looking north and east of the permeability model. The semantic rule-guided automatic optimized QS simulation method provides a more refined characterization of the spatial distribution of porosity and permeability.

Figure 15

Figure 15. 3D Reservoir Model. (a) Porosity model. (b) Permeability model.

Porosity is defined as the ratio of the total volume of pore spaces to the bulk volume of a rock sample. A higher total porosity indicates more developed pore spaces within the rock. Permeability, in turn, quantifies a rock’s ability to allow fluids to flow through it and is often strongly correlated with porosity. The significant correlation observed between the porosity and permeability models strongly validates the effectiveness of the structural constraints. Furthermore, comparing the boundary clarity and thickness variations in the generated models with the actual geological setting allows for an assessment of their ability to characterize the dynamic evolution of sedimentary facies.

To quantitatively evaluate the modeling performance, prediction accuracy and root mean square error (RMSE) were selected as the primary metrics. Prediction accuracy, which measures the proportion of correctly predicted samples, is calculated as shown in Equation 1:

$Accuracy=\frac{\left(TP+TN\right)}{\left(TP+TN+FP+FN\right)}(1)$

where $TP$, $TN$, $FP$, and $FN$ denote true positives, true negatives, false positives, and false negatives, respectively. RMSE, which quantifies the overall prediction error, is defined by Equation 2:

$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}(2)$

where $y_i$ is the observed value, ${\widehat{y}}_i$ is the predicted value, and $n$ is the number of samples. Lower RMSE values indicate better model fit. This study adopts the blind well validation method to evaluate the reliability and accuracy of the established model in predicting geological attributes. During the modeling process, a subset of samples was selected from the well logging data to serve as an independent test set, and the ${QH}_3$ well was designated as the blind well sample. The observational data from this well was not involved in simulation process and was reserved solely for subsequent validation. The prediction accuracy could be calculated by comparing the predicted values of porosity and permeability with the observational data at the blind well. To ensure stability and robustness of the results, ten independent simulation runs were conducted, with the prediction accuracy and error calculated for each run. The results indicate that the porosity model achieved an average prediction accuracy of 80.6% with an RMSE of approximately 1.29%, while the permeability model attained an average accuracy of 80.4% with an RMSE of approximately 0.94%. The variations in prediction accuracy and RMSE throughout the simulation runs are shown in Figure 16.

Figure 16

Figure 16. Variation in model prediction accuracy. (a) Porosity model. (b) Permeability model.

Meanwhile, to validate the superiority of the proposed method, comparative experiments were conducted using the Sequential Gaussian Simulation method and the conventional MPS method. For porosity modeling, the SGS method achieved a prediction accuracy of 74.3%, while the conventional MPS method reached 77.8%. For permeability modeling, the SGS method yielded a prediction accuracy of 74.4%, compared to 77.5% obtained by the conventional MPS method. The results demonstrate that the semantic rule-driven QS simulation method achieves a high prediction accuracy.

We systematically describe the principles and workflow of an improved MPS reservoir modeling method guided by semantic rules. To address the key challenge of 3D reservoir modeling under geological semantic constraints, we propose a structural modeling approach that integrates multi-source information with geological semantic rules, based on data from the actual study area. Using the resulting detailed structural model, we constrain the MPS simulation process. An improved QS algorithm is applied to achieve efficient parameter optimization and simulation execution. Finally, we complete the integrated fine-scale characterization of the structure-reservoir system in the study area, verifying the effectiveness and applicability of the proposed method.

4 Discussion

This work systematically investigates a semantic rule-guided 3D reservoir modeling method using an improved-MPS algorithm. Experimental validation demonstrates that the proposed approach significantly enhances model accuracy and geological applicability. However, as reservoir complexity increases, several key challenges persist. Among these, computational scalability becomes critical when applying the method to large-scale, field-wide models. The computational load of the proposed method arises primarily from two stages: first is the implicit structural modeling based on sequence stratigraphy and fault kinematics, which scales approximately linearly with grid size and constraint points, and another is the improved QS algorithm. While the QS algorithm accelerates pattern matching via FFT and reduces manual parameter tuning through neighborhood kernel function automatic selection strategy and fast optimization strategy for dual-dimensional parameters, its efficiency can still be challenged when processing large training images and evaluating numerous candidate patterns, leading to increased memory usage and computation time. Furthermore, incorporating detailed semantic constraints may add numerical complexity to the structural modeling stage.

Therefore, future research will focus on further optimizing algorithm implementation in large-scale parallel computing environments and exploring approximate methods for the parameter optimization process. In addition, subsequent work should emphasize the development of semantic rule-based multi-scale collaborative modeling methods. This involves establishing relational models that connect microscopic and macroscopic scales, enabling comprehensive characterization of dynamic reservoir features. By integrating multi-scale geological and physical constraints, the accuracy of model predictions regarding reservoir dynamic responses can be substantially enhanced, thereby providing more reliable data support for production evaluation, injection-production optimization, and long-term dynamic management during oil and gas field development.

5 Conclusion

This work addresses the challenge of integrating geological semantic information into 3D reservoir modeling by proposing a novel knowledge-data dual driven modeling paradigm. By combining the advantages of MPS simulation and implicit structural modeling, a novel 3D reservoir modeling method based on semantic rule-guided MPS simulation is proposed. The method deeply integrates multiple-source data with geological semantic rules to achieve structural modeling, and further to accomplish MPS-based 3D facies and property simulation under the constraints of the structural model. The proposed method and workflow were applied and validated through a case study conducted in a shale oil reservoir area. During the modeling process, geological semantic rules of sequence stratigraphy and fault kinematics were incorporated as knowledge constraints, thereby enhancing geological rationality. In addition, an improved MPS algorithm incorporating automatic domain kernel function selection and rapid dual-dimensional parameters optimization strategies were employed. This algorithm dynamically adjusts the kernel function and efficiently identifies the globally optimal parameter combination, thereby achieving a more realistic characterization of complex reservoirs. The experimental results demonstrate that the proposed 3D reservoir modeling method not only exhibits strong practicality but also shows completeness and operability in the overall workflow. By unifying knowledge-driven rules with data-driven simulation, this knowledge-data dual driven approach achieves accurate and detailed characterization of reservoir structures and properties, providing reliable support for hydrocarbon development strategy design.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

QC: Conceptualization, Funding acquisition, Supervision, Writing – original draft. HL: Methodology, Validation, Visualization, Writing – original draft. LX: Data curation, Visualization, Writing – original draft. DC: Methodology, Validation, Writing – review and editing. HF: Investigation, Writing – review and editing. GL: Funding acquisition, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work is supported by the National Natural Science Foundation of China (42172333 and 42372345).

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Arpat, G. B. (2004). “SIMPAT: stochastic simulation with patterns,” in 17th SCRF meeting, stanford Center for Reservoir forecasting. Stanford University.

Google Scholar

Baninajar, E., Sharghi, Y., and Mariethoz, G. (2019). MPS-APO: a rapid and automatic parameter optimizer for multiple-point geostatistics. Stoch. Environ. Res. And Risk Assess. 33, 1969–1989. doi:10.1007/s00477-019-01742-7

CrossRef Full Text | Google Scholar

Chen, Q., Mariethoz, G., Liu, G., Comunian, A., and Ma, X. (2018). Locality-based 3-D multiple-point statistics reconstruction using 2-D geological cross sections. Hydrology Earth Syst. Sci. 22 (12), 6547–6566. doi:10.5194/hess-22-6547-2018

CrossRef Full Text | Google Scholar

Cui, Z., Jiang, S., Chen, Q., Ma, X., Chen, D., and Liu, G. (2025). Characterization of Reservoir structures with knowledge-informed neural network. Soc. Petroleum Eng. J. 30, 4469–4486. doi:10.2118/228279-PA

CrossRef Full Text | Google Scholar

Feng, R., Grana, D., Mukerji, T., and Mosegaard, K. (2022). Application of Bayesian generative adversarial networks to geological facies modeling. Math. Geosci. 54 (5), 831–855. doi:10.1007/s11004-022-09994-w

CrossRef Full Text | Google Scholar

Foged, N., Marker, P. A., Christansen, A. V., Bauer-Gottwein, P., Jørgensen, F., Høyer, A. S., et al. (2014). Large-scale 3-D modeling by integration of resistivity models and borehole data through inversion, Hydrol. Earth Syst. Sci. 18, 4349–4362. doi:10.5194/hess-18-4349-2014

CrossRef Full Text | Google Scholar

Frank, T., Tertois, A. L., and Mallet, J. L. (2007). 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data. Comput. and Geosciences 33 (7), 932–943. doi:10.1016/j.cageo.2006.11.014

CrossRef Full Text | Google Scholar

Gravey, M., and Mariethoz, G. (2020). QuickSampling v1. 0: a robust and simplified pixel-based multiple-point simulation approach. Geosci. Model Dev. 13 (6), 2611–2630. doi:10.5194/gmd-13-2611-2020

CrossRef Full Text | Google Scholar

Gravey, M., and Mariethoz, G. (2022). AutoQS v1: automatic parameterization of QuickSampling based on training images analysis. Geosci. Model Dev. Discuss. 2022, 1–21. doi:10.5194/gmd-16-5265-2023

CrossRef Full Text | Google Scholar

Grose, L., Ailleres, L., Laurent, G., and Jessell, M. (2021). LoopStructural 1.0: time-aware geological modelling. Geosci. Model Dev. 14 (6), 3915–3937. doi:10.5194/gmd-14-3915-2021

CrossRef Full Text | Google Scholar

Guardiano, F. B., and Srivastava, R. M. (1993). “Multivariate geostatistics: beyond bivariate moments,” in Geostatistics Tróia ’92. Quantitative geology and geostatistics. Editor A. Soares (Dordrecht: Springer), 5, 133–144. doi:10.1007/978-94-011-1739-5_12

CrossRef Full Text | Google Scholar

Hou, W., Liu, H., Zheng, T., Shen, W., and Xiao, F. (2021). Hierarchical MPS-Based three-dimensional geological structure reconstruction with two-dimensional image(s). Earth Sci. 32, 455–467. doi:10.1007/s12583-021-1443-x

CrossRef Full Text | Google Scholar

Irakarama, M., Laurent, G., Renaudeau, J., and Caumon, G. (2021). Finite difference implicit structural modeling of geological structures. Math. Geosci. 53, 785–808. doi:10.1007/s11004-020-09887-w

CrossRef Full Text | Google Scholar

Kathuria, D., Mohanty, B. P., and Katzfuss, M. (2019). Multiscale data fusion for surface soil moisture estimation: a spatial hierarchical approach. Water Resour. Res. 55 (12), 10443–10465. doi:10.1029/2018wr024581

CrossRef Full Text | Google Scholar

Klise, K. A., Weissmann, G. S., McKenna, S. A., Nichols, E. M., Frechette, J. D., Wawrzyniec, T. F., et al. (2009). Exploring solute transport and streamline connectivity using lidar-based outcrop images and geostatistical representations of heterogeneity. Water Resour. Res. 45, W05413. doi:10.1029/2008WR007500

CrossRef Full Text | Google Scholar

Li, Z., Yi, S., Wang, N., Zhang, X., Chen, Q., and Liu, G. (2024). Adaptive direct sampling-based approach to ore grade modeling. Earth Sci. Inf. 17, 2537–2554. doi:10.1007/s12145-024-01297-4

CrossRef Full Text | Google Scholar

Liu, Z., Yu, X., Jia, L., Wang, Y. S., Song, Y. C., and Meng, H. D. (2021). The influence of distance weight on the inverse distance weighted method for ore-grade estimation. Sci. Rep. 11, 2689. doi:10.1038/s41598-021-82227-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, G., Fang, H., Chen, Q., Cui, Z., and Zeng, M. (2022). A feature-enhanced MPS approach to reconstruct 3D deposit models using 2D Geological Cross sections: a case study in the luodang Cu deposit, Southwestern China. Nat. Resour. Res. 31, 3101–3120. doi:10.1007/s11053-022-10113-z

CrossRef Full Text | Google Scholar

Mallet, J. L. (1992a). Discrete smooth interpolation in geometric modelling. Computer-aided Des. 24 (4), 178–191. doi:10.1016/0010-4485(92)90054-E

CrossRef Full Text | Google Scholar

Mallet, J. L. (1992b). “GOCAD: a computer Aided design program for geological applications,” in Three-Dimensional modeling with geoscientific information systems. NATO advanced Study institutes series. Editor A. K. Turner (Dordrecht: Springer), 354, 123–141. doi:10.1007/978-94-011-2556-7_11

CrossRef Full Text | Google Scholar

Mariethoz, G., Renard, P., and Froidevaux, R. (2009). Integrating collocated auxiliary parameters in geostatistical simulations using joint probability distributions and probability aggregation. Water Resour. Res. 45, W08421. doi:10.1029/2008WR007408

CrossRef Full Text | Google Scholar

Mariethoz, G., Renard, P., and Straubhaar, J. (2010). The direct sampling method to perform multiple-point geostatistical simulations. Water Resour. Res. 46 (11). doi:10.1029/2008WR007621

CrossRef Full Text | Google Scholar

Maxelon, M., Renard, P., Courrioux, G., Brändli, M., and Mancktelow, N. (2009). A workflow to facilitate three-dimensional geometrical modelling of complex poly-deformed geological units. Comput. and Geosciences 35 (3), 644–658. doi:10.1016/j.cageo.2008.06.005

CrossRef Full Text | Google Scholar

Nemes, I. (2023). Applications of automated petrel workflows in 3D reservoir geologic modelling–A case study. Cent. Eur. Geol. 65 (2), 122–136. doi:10.1556/24.2022.00129

CrossRef Full Text | Google Scholar

Okabe, H., and Blunt, M. J. (2007). Pore space reconstruction of vuggy carbonates using microtomography and multiple-point statistics. Water Resour. Res. 43, W12S02. doi:10.1029/2006WR005680

CrossRef Full Text | Google Scholar

Oriani, F., Straubhaar, J., Renard, P., and Mariethoz, G. (2014). Simulation of rainfall time series from different climatic regions using the direct sampling technique, hydrol. Earth Syst. Sci. 18, 3015–3031. doi:10.5194/hess-18-3015-2014

CrossRef Full Text | Google Scholar

Renaudeau, J., Malvesin, E., Maerten, F., and Caumon, G. (2019). Implicit structural modeling by minimization of the bending energy with moving least squares functions. Math. Geosci. 51, 693–724. doi:10.1007/s11004-019-09789-6

CrossRef Full Text | Google Scholar

Schaaf, A., de la Varga, M., Wellmann, F., and Bond, C. E. (2021). Constraining stochastic 3-D structural geological models with topology information using approximate Bayesian computation in GemPy 2.1. Geosci. Model Dev. 14 (6), 3899–3913. doi:10.5194/gmd-14-3899-2021

CrossRef Full Text | Google Scholar

Straubhaar, J., Renard, P., Mariethoz, G., Froidevaux, R., and Besson, O. (2011). An improved parallel multiple-point algorithm using a list approach. Math. Geosci. 43, 305–328. doi:10.1007/s11004-011-9328-7

CrossRef Full Text | Google Scholar

Strebelle, S. (2002). Conditional simulation of complex geological structures using multiple-point statistics. Math. Geol. 34, 1–21. doi:10.1023/A:1014009426274

CrossRef Full Text | Google Scholar

Tan, X., Tahmasebi, P., and Caers, J. (2014). Comparing training-image based algorithms using an analysis of distance. Math. Geosci. 46, 149–169. doi:10.1007/s11004-013-9482-1

CrossRef Full Text | Google Scholar

Wambeke, T., and Benndorf, J. (2016). An integrated approach to simulate and validate orebody realizations with complex trends: a case study in heavy mineral sands. Math. Geosci. 48, 767–789. doi:10.1007/s11004-016-9639-9

CrossRef Full Text | Google Scholar

Wang, L., Yin, Y., Zhang, C., Feng, W., Li, G., Chen, Q., et al. (2022). A MPS-Based novel method of reconstructing 3D reservoir models from 2D images using seismic constraints. Petroleum Sci. Eng. 209, 109974. doi:10.1016/j.petrol.2021.109974

CrossRef Full Text | Google Scholar

Wu, K., Van Dijke, M. I. J., Couples, G. D., Jiang, Z., Ma, J., Sorbie, K. S., et al. (2006). 3D stochastic modelling of heterogeneous porous media – applications to reservoir rocks. Transp. Porous Media 65, 443–467. doi:10.1007/s11242-006-0006-z

CrossRef Full Text | Google Scholar

Zhang, T. (2006). Multiple-point pattern simulation using filter scores: the filtersim algorithm and theory.

Google Scholar

Zhang, T., Li, D., Lu, D., and Yang, J. (2010). Research on the reconstruction method of porous media using multiple-point geostatistics. Sci. China Phys. Mech. and Astronomy 53, 122–134. doi:10.1007/s11433-009-0257-0

CrossRef Full Text | Google Scholar

Zhang, B., Du, L., Khan, U., Tong, Y., Wang, L., and Deng, H. (2023). AdaHRBF v1. 0: gradient-adaptive hermite–Birkhoff radial basis function interpolants for three-dimensional stratigraphic implicit modeling. Geosci. Model Dev. 16 (13), 3651–3674. doi:10.5194/gmd-16-3651-2023

CrossRef Full Text | Google Scholar

Zhong, D., Lin, B., Bi, L., and Jia, M. t. (2019). Implicit modeling of complex orebody with constraints of geological rules. Trans. Nonferrous Metals Soc. China 29 (11), 2392–2399. doi:10.1016/S1003-6326(19)65145-9

CrossRef Full Text | Google Scholar