Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1756046 10.3389/feart.2025.1756046 Original Research MVarGOSIM: an MPS algorithm for characterizing complex structures with multiple variables Chen et al. 10.3389/feart.2025.1756046 Chen Yonghua 1 Data curation Resources Supervision Writing - review and editing Hou Weisheng 2 3 4 * † Conceptualization Data curation Formal Analysis Funding acquisition Validation Visualization Writing - original draft Writing - review and editing Li Yanhua 2 5 * Formal Analysis Methodology Validation Writing - original draft Writing - review and editing Ye Shuwan 1 Data curation Formal Analysis Investigation Writing - review and editing Guangzhou Metro Design & Research Institute Co. Ltd., Guangzhou, China School of Earth Sciences and Engineering, Sun Yat-sen University, Zhuhai, China Guangdong Provincial Key Laboratory of Geological Processes and Mineral Resources Exploration, Zhuhai, China Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai), Zhuhai, China Three Gorges Renewables Yangjiang Power Co., Ltd., Yangjiang, China Correspondence: Yanhua Li, 1789376011@qq.com; Weisheng Hou, houwsh@mail.sysu.edu.cn

ORCID: Weisheng Hou, orcid.org/0000-0003-2669-6686

20 02 2026 2025 13 1756046 28 11 2025 16 12 2025 22 12 2025 Copyright © 2026 Chen, Hou, Li and Ye. 2026 Chen, Hou, Li and Ye https://creativecommons.org/licenses/by/4.0/ This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Creating a highly accurate geological model at a large scale presents a considerable challenge, primarily due to constraints imposed by sparse data availability. A promising strategy to mitigate these limitations involves the integration of multiple variables. Nevertheless, the effective incorporation and amalgamation of patterns derived from diverse variables during the simulation process remains a significant obstacle in pattern-based methodologies. This study presents a novel iterative multiple-point statistics (MPS) algorithm to construct complex geological structures by integrating co-located multiple variables in conjunction with fully connected deep artificial neural networks (FCNs). The algorithm operates under the assumption that profiles of different data types are co-located, allowing patterns from various variables to be converted into probabilities and combined using a cross-entropy-weighted pooling method. The proposed approach consists of three main components: generating geological subsurfaces using FCNs, creating an initial model incorporating multiple variables, and refining this initial model through an iterative process. The trained FCNs generate the top and bottom surfaces of a geological object, with a loss function defined by geological contact elevation. In the initial model construction, patterns from co-located data are integrated with pattern probabilities, using lithology cross-sections as the primary variable and velocity and density profiles as auxiliary variables. Geological constraints, such as stratigraphic sequences and the thickness of geological objects, are applied in a post-processing phase to adjust the relationships in the initial model. An expectation–maximization-like (EM-like) optimization technique is used to rectify artifacts present in the initial model. The efficacy of the proposed algorithm is demonstrated by reconstructing the overthrust model developed by SEG/EAGE. Comparative analyses between the reference model and the results obtained with and without multiple variables indicate that the proposed algorithm achieves a more accurate representation of geological objects while also better preserving their geometry and interrelationships.

 cross-entropy expectation–maximization-like algorithm geological constraints multiple-point statistics multiple variables The author(s) declared that financial support was received for this work and/or its publication. This research was jointly supported by the National Natural Science Foundation of China (NSFC) Program (42372341 and 41972302) and the Guangdong Province-Introduced Innovative R&D Team of Big Data—Mathematical Earth Sciences and Extreme Geological Events Team (2021ZT09H399). section-at-acceptance Geoinformatics Introduction

A three-dimensional (3D) realistic geological model offers a perspective view and data framework for studying the earth’s evolution and practical applications, such as oil and gas reservoirs (Feng et al., 2017), mineral exploration (Wu et al., 2015; Xiao et al., 2015), civil engineering (Guo et al., 2022; Hou et al., 2023), and flow and solute transport (Bai and Tahmasebi, 2020). Since the first algorithm was proposed in 1993 (Guardiano and Srivastava, 1993), multiple-point statistics (MPS) has developed into a prominent modeling technique across various fields as it accommodates high-order spatial characteristics (Mariethoz and Caers, 2014) and mitigates the shortcomings associated with two-point statistics-based modeling methods. MPS creates a 3D model by reproducing patterns through a stochastic process, where these patterns are derived from known data by moving templates, referred to as training images (TIs). In current MPS simulations, two-dimensional (2D) profiles are commonly used to generate 3D geological models (Comunian et al., 2012; Gueting et al., 2018; Hou et al., 2023), although certain MPS-based modeling methods require 3D TIs (Wang et al., 2022).

Recent studies have shown that reconstructing complex geological structures with integrated multiple data sources can reduce uncertainties in the final outcomes (Hansen et al., 2018; Liu et al., 2004; Zhao and Zhou, 2019). Integrating multiple TIs into a single TI is a practical and straightforward strategy in the MPS-based method. An effective approach is to transform TIs into “soft data.” The SNESIM/ENSIM and DS methods utilize a simulation path that favors the 1D marginal entropy/information content of soft data, allowing them to manage scattered soft data (Hansen et al., 2018). Two position-independent TIs with varying anisotropy structures and attributes are integrated with a two-point correlation function (Hasanabadi et al., 2019). Additionally, several candidate data that are not location-specific can be integrated by assessing candidate TIs and optimally selecting based on minimum data event distances (MDevDs) (Feng et al., 2017). Other types of information are also considered in the reconstruction of groundwater models using MPS algorithms (He et al., 2014; Koch et al., 2014). Some approaches address the issue by integrating conditional probabilities transformed from multiple variables. By generalizing normal scores and back transforming for multivariate distributions, along with decomposing the target distribution into a product of univariate marginal and conditional probability density functions, a direct simulation technique enhances efficiency when dealing with more than three variables in two dimensions (Figueiredo et al., 2021). An MPS algorithm was proposed to integrate non-colocational categorical conditioning data using a probability distribution that describes the available information for each model parameter (Jóhannsson and Hansen, 2023). In the aforementioned methods, only one TI that integrates multiple variables is used to generate realizations in the simulation process. However, geological constraints such as geo-event sequences were not incorporated into the modeling process, and the distinct characteristics of different variables were not adequately considered in the process of selecting candidate patterns.

This study presents an iterative MPS algorithm that incorporates co-located multi-variables for reconstructing complex 3D geological structures. The proposed algorithm adopts the core framework proposed by Hou et al. (2025), comprising three key components: geological subsurface reconstruction via the fully connected deep artificial neural network (FCN), generation of an initial model incorporating multiple variables, and iterative refinement of the initial model through an expectation–maximization-like (EM-like) process. A primary contribution of this work lies in the multi-variable integration strategy implemented during the initial model generation phase. Under the assumption that diverse geological and geophysical data profiles are spatially co-located, patterns derived from each variable are converted into probability distributions. These distributions are then integrated using a cross-entropy approach. The remainder of this article is structured as follows: Section 2 details the kernel function and core methodologies underpinning the proposed algorithm. Section 3 presents a case study demonstrating the reconstruction of an overthrust model using integrated density, velocity, and lithology data. Finally, Section 4 provides discussion and conclusions.

 Methods Basic idea of the proposed algorithm

In accordance with the MPS-based modeling method presented by Hou et al. (2025), the proposed algorithm comprises three distinct phases, as illustrated in Figure 1. The first phase involves data pre-processing and the construction of surfaces using 2D geological cross-sections. To obtain 3D patterns from different datasets, two-dimensional profiles with zero thickness are expanded into 3D training data (TD), with the thickness defined as the edge length of a cubic template; the expansion process is performed as described by Hou et al. (2025). Geological constraints, including stratigraphic sequences, thickness, and spatial relationships among geological objects, are extracted and used in the simulation process. A trained FCN, in which the kernel function is defined by the elevations of a geological object, is used to generate the lithological surfaces. In the second phase, an initial model at the coarsest scale is constructed, in which multiple variables and geological constraints are used to select candidate patterns. This phase addresses two primary concerns: the identification and integration of multiple patterns derived from various datasets for assigning attributes between the geological surfaces of different lithologies and the validation of geological relationships among objects in accordance with the geological constraints. The final phase involves the optimization of the initial model via an iterative MPS process, which incorporates a multi-scale EM-like approach. Consequently, the difference between the method in this study and the algorithm presented by Hou et al. (2023) and Hou et al. (2025) lies in the construction of an initial model using lithological cross-sections, density profiles, and velocity profiles.

Flowchart of the proposed algorithm.

Flowchart depicting the process of building a geological model. It begins with data sources showing lithology profiles for velocity and density. This information, combined with geological constraints and a trained Fully Convolutional Network (FCN), creates lithological surfaces. An initial model is constructed by integrating multiple variables. This model undergoes EM-like iteration, resulting in a refined final model for optimization. Kernel function of MPS simulation with multiple data sources

Assuming that the final result has patterns as similar as possible to those in the TIs, the objective function of the MPS, termed the pattern distance, is presented in Equation 1 (Hou et al., 2023): d R , T D = ∑ P R f i n a l ∈ R min P T I ∈ T I D P R , P T D , where P R , P R f i n a l , and P TD are patterns in a realization R , final realization R final , and 3D TD for modeling, respectively; and D ( P R , P TD ) is the distance between P R and P TD . Considering constraints such as the stratigraphic sequence in the simulation, Equation 1 is modified to Equation 2: R f i n a l = arg   min R ∑ D P T D , P R | S R ∈ D S , where S( R ) and D S are the geological relationships in R and the constraint database from 3D TD, respectively. In Equation 2, only P TD and D S are known, and the other two variables are unknown. Therefore, minimizing the Euclidean distance D ( P TD , P R ) in the feature space is a typical undetermined problem. To solve Equation 2, the EM-like algorithm is used with an initial model to facilitate the iterative solution.

 Key steps of MPS simulation with integrating multiple variables

The presented methodology, which incorporates multiple variables, consists of several essential steps, including sequential simulation, retention of geological constraints, integration of multiple variables, and an EM-like iterative process. A significant challenge within the proposed algorithm is the generation of coherent 3D patterns as the identification and selection of candidate patterns are based on similarity. To mitigate this issue, an intuitive expansion procedure using 2D cross-sections is implemented, as illustrated by Hou et al. (2023). This procedure is conducted for each variable, commencing from the nodes adjacent to the 3D TD, layer by layer from top to bottom. Upon completion of the expansion process, the profiles of various variables are transformed into 3D sections, with a thickness corresponding to the edge length of a cubic template. Subsequently, three-dimensional patterns can be extracted from the 3D TD and categorized based on similarity. For each variable, three-dimensional patterns are clustered.

Sequential process for building an initial lithological model

Here, the initial model is reproduced by a sequential assigning process, where each geological block is constrained with bottom and top surfaces generated by FCNs. Considering that the geometry of geological objects is usually represented by contact elevations, the objective function of the FCN is presented in Equation 3 (Hou et al., 2023): L o s s F C N = M S E E , E ′ = ∑ x = 0 x n u m ∑ y = 0 y n u m e x , y , a t t i − e ′ x , y , a t t i 2 x n u m × y n u m , where xnum and ynum are grid numbers of the simulation grid (SG) in the lateral plane; att i is the ith lithological type; e (x, y, att i ) and e ′ x , y , a t t i are the elevations of att i in 3D TD and desired elevation obtained by the FCN at coordinate (x, y), respectively. For each lithology, two FCNs are constructed to generate the top Surtop (att i ) and bottom Surbott (att i ) surfaces and are trained using geological contact elevations. The attributes between Stop (att i ) and Sbott (att i ) are assigned sequentially to obtain the initial model.

In this study, Surtop (att i ) and Surbott (att i ) of the blocks of a lithological attribute att i are generated using two networks N top (i) and N bott (i), respectively, on a back-propagation (BP) neural network architecture. This architecture comprises 8 hidden layers, 1,300 neurons, and 261,451 parameters, as provided by Hou et al. (2025). The coordinates (x, y) and elevations of geological contacts of att i are the input and output of the FCNs for attribute att i , respectively. The elevations are normalized before building the training and validating datasets. The training and validation datasets occupy 80% and 20% of the total dataset, respectively. The FCN training process is terminated when the number of epoch reaches 10,000 or the loss value is smaller than 5 × 10−6. Then, Surtop (att i ) and Surbott (att i ) at the coarsest scale are predicted using the FCNs.

The lithological attributes on the grids between Stop (att i ) and Sbott (att i ) are assigned as att i using a sequential process. Here, the surfaces are generated twice using different trained FCNs if the surfaces are not the boundary of the simulation area. Therefore, the surfaces of different objects may intersect, resulting in multiple attributes being assigned to a single grid. Then, grids with several attributes are labeled as unassigned grids after clearing the corresponding attributes. The model with unassigned grids is marked as R 0 . To assign proper values to unassigned grids in R 0 , a sequential process is implemented as follows:

An unassigned grid neighboring grids with attributes in R 0 is randomly selected and labeled as the current node u c . A 3D moving window W u with the size of W L × W L × W L is defined, in which the center of W u is u c . Moreover, the number of grids with values in W u usually should be larger than W L × W L × (W L /2–1) or a user-defined value.

A candidate pattern for u c is obtained. Patterns with the maximum similarity between the pattern P R 0 u of W u and the pattern in the pattern database of lithology are selected, where the similarity is calculated using Equation 5. Let the number of candidate patterns be N; the probability Pr P T D n of each candidate’s pattern P T D n is used as the parameter to determine the candidate pattern. Pr P T D n is calculated using Equation 4: Pr P T D n = 1 D P T D n , P R 0 u ∑ n = 1 N 1 D P T D n , P R 0 u , n = 1 , 2 , 3 , ⋯ N , where D P T D n , P R 0 u is the Hamming distance between P T D n from 3D TD and P R 0 u from the initial model. The pattern with the maximum Pr P T D n is pasted on u c in R 0 . If several patterns have the same maximum value, one of them is randomly pasted on u c . Then, u c moves to another unsigned grid in R 0 , if it exists.

Steps (1) and (2) are repeated until all unsigned grids in R 0 are assigned with values. The model obtained here is labeled with R 0 ′ .

Retaining geological spatial relationships

The initial model R 0 ′ is constructed using a stochastic process without any extra constraints, resulting in wrong relationships. To rectify errors in R 0 ′ and retain reasonable relationships among geological objects in the initial model, geological constraints, including stratigraphic sequence, thickness, and spatial relationships of geological objects, are applied in the presented method.

In 2D cross-sections, the stratigraphic sequence and spatial relationships of faults and strata are implied in the spatial structures. Here, the 2D geological cross-sections are converted into image files, such as BMP or PNG files, and imported into the 3D simulation grid. Then, a scanning process is carried out column by column along the z-direction. A sequence S (x, y) of geological objects in a column is extracted from top to bottom, and the details of the extracting process are described by Hou et al. (2023). If no strata overturn occurred in the profiles, the stratigraphic sequence would be the union set U(S (x, y)) of S (x, y). Otherwise, the stratigraphic sequence should be built manually before extracting the relationships because the extraction process cannot obtain the right deposition sequence exactly due to strata overturn. If intrusive rocks occur in the profiles, the latest active time is identified by the latest intruded stratum. Then, the spatial relationships of all geological objects can be extracted and reorganized with a digraph. For example, in Figure 2, geological objects including strata such as Sa, Sb, Sc, Sd, and Se; intrusive rock Ina; and fault F1 are shown. The black, brown, and red arrows show strata contact, intrusion contact, and fault contact, respectively. The black arrows indicate the old formations within the strata. Therefore, in Figure 2, the stratigraphic sequence is “Sa→Sb→Sc→Sd→Se.” Because Sa and Sc are connected by black arrows, Sb and Sd are pinched out in the profiles. The F1 cuts through Sb to Se. The intrusive rock Ina intrudes through Sc to Se. In the cross-sections, the strata on both sides of a fault do not contact directly, indicating that the fault cannot be represented as a surface with zero thickness.

Example of organizing the relationships of geological objects.

Diagram of geological features showing a sequence from strata Sa to Se with arrows indicating strata contact. An arrow from Ina denotes intrusion affecting Sd. Red arrows indicate fault contacts pointing to fault F1. A legend defines symbols for strata, intrusion, and fault contacts.

Upon completion of the sequential simulation process, the interrelationships among geological objects in R 0 ′ are evaluated with respect to geological constraints. At any given position (x, y) in R 0 ′ , the attribute sequence seq (x, y) is extracted from the vertical grid layers. If s e q x , y ⊄ D s and seq (x, y) is null, relationship errors might occur at position (x, y) along the vertical direction. Consequently, the attributes of the grids along the vertical direction on (x, y) are cleared and labeled as unassigned values. The sequential process described in Section 2.3.1 is repeated until all relationship errors among geological objects in R 0 ′ are resolved.

In this study, Stop (att i ) and Sbott (att i ) are also generated stochastically, which may result in an unrealistic thickness of a geological object, as shown in Figure 3a. If it is greater than 1.5 times or smaller than half the average value in cross-sections, a geological object’s thickness is modified to the average value (as shown in Figure 3b). After the thicknesses and relationships of all geological objects are retained, the corrected R 0 ′ with the coarsest scale is output as the initial model R 0.

Schematic chart of strata thickness before (a) and after (b) adjusting. The circle in the black dashed line shows the differences between (a) and (b).

Two illustrations labeled (a) and (b) depict geological layers with various colors: brown, yellow, green, blue, and pink. Both images show layered structures with alterations in the blue and pink layers within black dashed ovals, suggesting changes in subsurface geology. Integrating patterns from multiple variables in the initial model

This study synthesizes characteristics derived from geological cross-sections, seismic velocity profiles, and density sections to construct an initial model, with geological cross-sections serving as the primary variable. All known data are assumed to be co-located, ensuring that each pair of velocity, density, and geological cross-sections occupies the same spatial position and extent. The candidate pattern is identified based on the similarity between the overlapping regions of the data events and potential patterns derived from 3D TD. In this study, a distinct distance function is used to assess similarity. For categorical variables, such as lithology represented in the cross-sections, the Hamming distance function is used, as shown in Equation 5: D p i , p j = ∑ h = 1 L c u b e ∑ x = 1 L c u b e ∑ y = 1 L c u b e p i x , y , h ⊕ p j x , y , h , where L cube is the edge length of the template cube, and p i and p j are the i th and j th patterns, respectively. The Euclidean distance function is used for continuous data, such as velocity and density: D P T D , P R = ∑ h = 1 L c u b e ∑ x = 1 L c u b e ∑ y = 1 L c u b e P T D h , x , y − P R h , x , y 2 .

Because different variables have different dimensions and magnitudes, the distance calculated from Equations 5, 6 cannot be directly integrated. Therefore, all distance values are normalized before integration. Then, patterns from different data can be readily integrated using probabilities.

Let Patt main (i) be the i th pattern with the central position Ptmain (xi, yi, zi) from the geological cross-sections, Patt j (xi, yi, zi) be the pattern from the j th type of data on the same position, and the distance between Patt main (i) and Patt j (i) be D (i, j); then, the similarity probability of the i th pattern of the j th variable of the current data event is presented in Equation 7: Pr P i j = 1 / D i , j ∑ i = 1 P n u m D i , j , where Pnum is the pattern number.

In this study, a linear pooling method is used to integrate different data with the cross-entropy weight. The cross-entropy H(S, Q) can quantify the differences between two probability models S and Q based on their characteristics. In this study, the main variable S is the lithology from the geological cross-section, and the auxiliary variables Q are the density and velocity profiles. H(S, Q) is calculated as shown in Equation 8: H S , Q = ∑ i = 1 P n u m S i * log ⁡ Q i , where Pnum is the pattern number from the geological cross-sections, S (.) and Q (.) are the probabilities from geological cross-sections and from density or velocity profiles, respectively. A low H(S, Q) value indicates that the similarity between patterns from two variables is high. The corresponding weighting value of the j th type of data is presented in Equation 9: w j = 1 − H j S , Q ∑ j = 1 T n u m 1 − H j S , Q , where Tnum is the number of variables. With the selected patterns and corresponding weightings, multiple variables can be integrated using the linear pool method. Therefore, the objective function of MPS-based 3D modeling with the integration of multiple variables can be expressed as follows: P R i n i t i a l = argmax ∑ j = 1 T n u m w j × Pr P i j .

Considering the geological constraints used in the simulation, Equation 10 should be represented as Equation 11: P R i n i t i a l = argmax ∑ j = 1 J w j × Pr P i j | S R ∈ D s .

The initial lithological model presented here is developed through a two-step process. In the first step, the trained FCNs generate the lithological surfaces. The second step involves assigning attributes to the grids located between these lithological surfaces, using a sequential methodology that incorporates multiple variables.

R 0 is a model constructed only using lithological cross-sections. Here, the variables are integrated as follows (shown in Figure 4):

Setting simulation path. A random or sequential path can be chosen in this study. Using a random path, the grids to be simulated are extracted and reorganized in a random order within the grid set, with grids containing more conditional data given priority during the simulation. The grids to be simulated are scanned in sequence in the sequential path.

Obtaining the candidate patterns. For the currently visited grid u to be simulated, the data events P u s of the lithology, centered at u, are extracted. According to the similarity of the overlapped area P u s and the patterns in the pattern database of lithology, the candidate patterns Pc and P u s are chosen, and the corresponding center positions of Pc and ( P u s ) are saved as Pos(Pcand). Then, the corresponding patterns of auxiliary variables, including velocity and density, on Pos(Pcand) are extracted and marked as Pcand ( P u v ) and Pcand ( P u ρ ), respectively. All those patterns are stored in a temporary database.

Calculating the probability matrix. The distances between the data events and potential candidate patterns of each variable are calculated, in which the Hamming distance is used for lithology and the Euclidean distance is used to calculate the pattern similarity of the other two variables. Then, the distance vector ( D u s , D u v , and D u ρ ) is obtained with dimensions of (1 × n), where n is the number of candidate patterns. The probabilities Pr u s , Pr u v , and Pr u ρ of each variable are calculated using Equation 7, resulting in a probability matrix M (Pr) with dimensions of (3 × n).

Choosing the pattern for integrating probabilities with cross-entropy. The cross-entropy H ( s , s ), H ( s , v ), and H ( s , ρ ) of three variables are calculated. Then, the entropy weighting matrices [ W u s , W u v , and W u ρ ] of the three variables are obtained using Equation 9 and multiplied by the matrix M (Pr). The pattern of lithology with the maximum probability is selected and pasted on the node u.

Steps (2)–(4) are repeated until all the nodes are traversed. The result is output as the initial model R 1 at the coarsest scale.

Algorithm of the sequential procedure with multiple variables.

Flowchart detailing a sequential procedure with multiple variables in a simulation. Steps include setting the simulation path, obtaining candidate patterns, calculating a probability matrix, and choosing and pasting the candidate pattern. Specific instructions are provided at each step, involving extraction of lithology data, calculation of probabilities and distances, and cross-entropy values. The procedure concludes at step six. EM-like iterative process with a multiple-scale strategy

The structures in the initial model R 1 are reorganized using an intuitive pasting operation, in which the relationships among objects are conditioned by geological constraints. However, some unreasonable structures, including discontinuities and manual artifacts, appeared in R 1 . Here, an EM-like iterative algorithm with a multiple-scale strategy is carried out to optimize R 1 .

The optimization process is implemented starting from the coarsest-scale model using the EM-like algorithm. After optimizing the model at a given scale, it is subsequently upsampled to the next scale, and this process continues until the finest scale has been optimized. This section provides a concise overview of the EM-like algorithm, with further details available in the work by Yang et al. (2016). The EM-like process at each scale includes an E-step and an M-step. In the E-step, a pattern P R 1 u on an arbitrary grid u in R 1 is selected, in which u is labeled as the current grid. A pattern P D u in the lithological pattern database is randomly selected and pasted onto the grids, with the central grid being u. The pattern with the maximum distance D ( P D u , P R 1 u ) between P R 1 u and P D u from the lithological pattern database is selected as a candidate P D u . Then, a pattern P D u ′ for grid u is randomly chosen from a searching window, with the central grid corresponding to the center of P D u in 3D TD. If D ( P D u ′ , P R 1 u ) > D ( P D u , P R 1 u ), P D u ′ replaces P R 1 u . After all windows around u have been searched, the size of the searching window is reduced to the preset value. The process is terminated when the search window is smaller than the preset template. When the E-step is ended, the M-step is implemented once to update the model. In the E-step, only one candidate pattern is selected in each cross-section. The pattern on u is replaced with the candidate pattern, with the maximum occurrence of the attributes.

 Data and results Data source

In this study, the proposed algorithm is tested using an overthrust model developed by EAGE/SEG. The dataset consists of the 3D subsurface velocity, density models, and geological surfaces. A central part with a size of 10 km × 10 km × 4 km is extracted as the reference model. Ten mutually perpendicular cross-sections of subsurfaces, along with their corresponding density and velocity measurements, have been extracted as data sources for modeling purposes, as illustrated in Figures 5, 6.

Cross-sections extracted from the reference model. B1–B8 are virtual boreholes for model validation.

3D geological model showing a grid with colorful stratified layers representing different geological formations such as limestone, marl, sandstone, and faults labeled A, B, and C. Green lines indicate boundaries within the model, with labels B1 to B8 marking specific points. A legend on the right correlates colors to formations, and a 3D axis indicates orientation.

Cross-sections of density (a) and velocity (b) from the reference model.

3D models showing subsurface geological layers divided into grids. Panel (a) illustrates color-coded density from two to 2.8 grams per cubic centimeter using a green to brown scale. Panel (b) displays velocity data from two thousand to six thousand meters per second in a blue to red scale.

In the reference model, 12 deposits and 3 faults are developed. The lithology from top to bottom consists of Weathering zone (WZ), delta plain 1 (Dp 1), Marl 1, Shale 1, Limestone 1, Marl 2, Limestone 2, Marl 3, Sandstone 1, Delta plain 2 (Dp 2), Dolomite, and Basement. Three overthrust faults cut all deposits in the middle of the model.

 Simulation results

The final model is derived from an initial model characterized by a grid size of 200 × 200 × 80, which is subsequently refined to a grid size of 400 × 400 × 160, using a template size of 7 × 7 × 7 grids. The initial model, which uses a straightforward stochastic process constrained by surfaces generated using FCNs, exhibits chaotic stratification, as shown in Figure 7a. In contrast, the strata and faults in the subsequent model with geological constraints (Figure 7b) demonstrate greater continuity, relying solely on lithological cross-sections as data inputs. Following the rectification of relationships with geological constraints, certain anomalies are addressed. The contacts depicted in the final model (Figure 7c) exhibit a smoother appearance compared to those in Figure 7b, following an EM iterative process. Furthermore, when integrating density and velocity cross-sections, the resulting surface (shown in Figure 7d) is observed to be more continuous and smoother, despite the presence of some scattered strata.

Initial and final models with lithological cross-sections (first row) and multiple variables (second row). Image (c) is the final result from the initial model using the stochastic process result (a) and with geological constraints (b). Image (e) shows the final model optimized from the initial model (d) with multiple variables. “Null” indicates that the grid is labeled with an unassigned attribute caused by multiple values in one grid.

Five 3D geological models, labeled (a) to (e), depict layers of various rock types, including WZ, sandstone, dolomite, shale, limestone, marl, and faults. The models show a cross-section view with color-coded strata. Each model highlights different aspects and configurations of geological interfaces. A color key indicates each layer type and fault.

The discrepancies observed in the final results (Figures 8a–c) exhibit a continuous distribution, which is rigorously constrained by the cross-sectional data. For example, in the y-direction, Fault B does not reach the modeling boundary as it is not represented in the cross-section at y = 400. Similarly, Fault C is only present in the model for positions where y > 100 (as illustrated in Figure 8c), consistent with the cross-sectional representations (Figure 8). Nevertheless, artifacts associated with the faults persist in the profile positions, as depicted in Figure 8d. This occurrence is attributed to the intrusive expansion of geological object geometries, which remain unrefined in subsequent iterative processes. Consequently, the strata adjacent to the faults are not directly connected when constrained by the faults in the cross-sectional views (Figure 8d).

Modeling results of Fault A (a), B (b), and C (c) and Marble 1 (d).

Four 3D plots show geological layers in different colors. (a) Red terrain with gridlines. (b) Purple terrain with contours. (c) Pink terrain with undulating surface. (d) Combined view with red, pink, and purple layers intersecting. Each plot has x, y, z axes. Model validation

Slices were extracted along the y-direction at intervals of 20 pixels, ranging from 200 to 300 pixels, as illustrated in Figure 9. Except for the slices at y = 200 and y = 300, the remaining slices represent simulation results. The images presented in the left and right columns correspond to slices derived from models based on lithological cross-sections and multiple variables, respectively. Notable discrepancies are observed at the contacts within the simulation results that utilize different data sources. Additionally, unexpected phenomena are evident within the ellipses delineated by white dashed lines in Figure 9b; for instance, a specific error is absent in the results derived from multiple variables, as shown in Figure 9h. A similar pattern is observed in Figures 9i–k, which are also highlighted with white dashed ellipses.

Slices on different y positions from simulation results by lithological cross-sections (left column) and integrating density and velocity (right column). The ellipses in white dashed lines mark the differences between slice pairs (a–l).

A series of color-coded contour plots displaying variations across different datasets, labeled from (a) to (l). Each plot shows wavy, multicolored patterns with gamma values ranging from two hundred to three hundred, incrementing by twenty. The patterns and colors vary slightly between the plots, signifying changes in data characteristics.

To assess the validity of the constructed 3D geological model, the proportions of strata within the models and the 3D TD are presented in Table 1. Except for Mar 1, Limestone 1, and the Basement, the proportions of geological entities derived from the simulation results using multiple datasets exhibit a closer alignment with the proportions found in the 3D TD compared to those obtained solely through lithological analysis.

Proportion of geological objects in 3D TD and models with and without multiple variables.

Geological object 3D TD Model with lithology Model with multiple variables
WZ 5.28 5.14 5.23
Dp 1 2.74 3.15 2.76
Mar 1 6.18 5.72 5.70
Shale 1 10.18 9.95 11.09
Limestone 1 5.85 5.81 5.64
Mar 2 5.26 5.73 5.22
Limestone 2 4.56 4.03 4.80
Mar 3 9.26 10.50 9.22
Sandstone 1 3.88 4.52 3.58
Dp 2 9.76 9.32 9.91
Dolomite 9.92 9.33 9.87
Basement 18.54 18.54 18.52

Unit: %.

In this study, lithological sequences were extracted from both reference and reconstructed models at eight designated locations (in Figure 5), as illustrated in Figure 10. Although the lithological sequences derived from the simulation results exhibit similarities to those of the reference model, notable discrepancies are present. In particular, Faults A and B, depicted in the simulation results (left column of sub-images in Figure 10), are absent in the reference model at locations B1, B3, B5, and B7 (right columns of sub-images in Figure 10). The stratigraphic sequence observed between the two fault contacts aligns with the stratigraphic sequence presented in the cross-sections. Despite the rectification of geological relationships during the initialization phase, certain unforeseen distributions of geological objects persist in the final model.

Comparison of the attributes in the virtual boreholes from the simulation result (left column in each sub-image) and the reference model (right column in sub-image), respectively.

Vertical stratigraphic columns labeled B1 to B8 display various colored strata representing different geological layers such as WZ, Dp1, Marl1, and others. Each column is marked with color-coded squares, with a legend indicating layer types and faults, including Fault A, B, and C.

To quantitatively assess the accuracy of the geometric representation of objects in the simulation results, the dispersion metric proposed by Shi and Wang (2022) is used. This metric is predicated on calculating the proportion of corresponding attributes of virtual boreholes, evaluated from top to bottom in a vertical orientation. The dispersion can be articulated as shown in Equation 12 (Shi and Wang, 2022): D p = ∑ i = 1 N I a t t r i = a t t m i / N , where attr(i) and attm(i) are the attributes of the ith grid from the reference model and the simulation result, respectively; I (·) is an indicator function that equals 1 when attr(i) and attm(i) have the same value and 0 otherwise; N is the grid number of the virtual borehole. The dispersion value is confined to the interval [0, 1]. A higher dispersion value indicates greater precision in the simulation outcomes. As illustrated in Figure 10, the minimum and maximum dispersion values are recorded at 65.6% and 85.0%, respectively. Table 2 presents the dispersion values of eight boreholes, comparing results obtained with multiple variables to those derived from lithology alone. Notably, the average dispersion of results incorporating multiple variables is 2% greater than that of results based solely on lithological data.

Dispersion of boreholes from simulation results with lithology and those with multiple variables.

Borehole Results with lithology Results with multiple variables
B1 75.6 77.5
B2 82.5 85.0
B3 68.8 70.6
B4 67.3 66.2
B5 65.4 65.6
B6 82 84.4
B7 68.1 68.1
B8 78.5 80.0

Unit: %.

Figure 11 shows variogram functions derived from the simulation results using identical parameters across various directional orientations. This figure encompasses three results generated with multiple variables and twenty results based solely on lithology. The variograms associated with the simulations using multiple variables exhibit a closer alignment with the 3D TD compared to those derived from lithological data alone. This observation suggests that simulations incorporating multiple variables are more effective in accurately replicating the characteristics of patterns observed in the 3D TD.

The variogram functions of 3D TD and simulation results with lithology and integrating multiple variables in the x- (left), y- (middle), and z-directions (right).

Three line graphs display data trends comparing "Multiple variants" in red, "Lithology" in gray, and "3D TD" in blue, against voxel lag. The graphs show similar upward trends and variations in all panels, with specific peaks and troughs. Discussions

The proposed algorithm has been evaluated using various pattern sizes and numbers of candidate patterns, as detailed in Table 3. The results indicate that smaller pattern sizes yield more continuous representations of objects. For example, Shale 1, depicted in the results with a pattern size of 7 × 7 × 7 (illustrated in Figures 12a–d), exhibits a thinner profile than the model using a pattern size of 5 × 5 × 5 (shown in Figures 12e–h). Each segment of the strata, which is intersected by three faults, demonstrates a continuous distribution characterized by uneven contacts and pinching at the fault surfaces (as shown in Figures 12b,d,f,h). Despite the presence of certain anomalies within the model, there is no significant confusion regarding the contact relationships among the geological entities. Conversely, when the number of candidate patterns is set to 5, the uncertainty in the results appears to increase. For example, Figures 12c,d illustrate that Sandstone 1 maintains a continuous distribution and interacts with faults A and B in the model with a pattern size of 7 × 7 × 7. In the slices at y = 250, where the study area is intersected by three faults, Sandstone 1 extends in the y-direction with minimal dip angles, resulting in indented boundaries (as depicted in Figures 12a,c). The model also indicates the erosion of thrust objects, as highlighted by the white dashed lines in Figures 12b,d.

Different parameters for testing the proposed algorithm.

Parameter type Pattern size Candidate pattern number
Para 1 7 × 7 × 7 3
Para 2 7 × 7 × 7 5
Para 3 5 × 5 × 5 3
Para 4 5 × 5 × 5 5

Simulation results with different parameters. The images in different rows show the results with parameters of Para 1, Para 2, Para 3, and Para 4, respectively. Images in the left column are the slices on y = 250 of simulation results with different parameters shown in Table 3. Images in the right column are 3D models of Shale 1, Limestone 1, and three faults by various parameters (a–h).

Six side-by-side images show geological models and cross-sections. The left column displays colorful stratified layers with dashed white ovals highlighting specific areas. The right column features corresponding 3D models with similarly marked sections. Each row represents a different stage or variation in the models, labeled (a) through (h). Axes are indicated in each model.

This study identifies certain implausible scenarios within the final model. For instance, as illustrated in Figures 12a,c,e a small segment of the Weathering Zone is encased by Delta Plain 1 on the y = 250 slice, with the exception of the outcome derived from the parameter Para 4. This discrepancy arises because the stratigraphic sequence is not verified during the EM iteration process. Consequently, the unconstrained iterative stochastic process may yield erroneous results. It is imperative that future research seriously consider an appropriate methodology for preserving geological relationships throughout the iterative process.

The proposed approach emphasizes the importance of searching for and selecting candidate patterns, which is crucial for achieving favorable outcomes based on the similarities between data events and potential patterns. When an optimal solution is found, only one candidate pattern is identified. However, using multiple cross-sections significantly affects the selection of candidate patterns and, consequently, the final results.

The multiple variants are integrated in constructing the initial model. Some small differences between the initial models constructed with and without multiple variants can be observed visually in Figures 7b,d. The multiple-dimension scaling (MDS) plot illustrates that there are significant differences between those initial models (Figure 13). The initial models, obtained by integrating multiple variants, do not cluster with other initial models. Therefore, some different patterns can be obtained in the initial models constructed by integrating multiple variants.

MDS plot of initial models constructed using lithology and by multiple variables.

Scatter plot with blue and pink dots. Blue dots represent the "initial model by lithology," clustered in the bottom right corner. Pink dots represent the "initial model by multiple variants," scattered vertically on the left side.

During the simulation, the number of each pattern chosen from various cross-sections is calculated, as illustrated in Figure 14. In the lithology simulation, the selected patterns are unevenly distributed along the faults, with a higher likelihood of selection occurring at fault contacts and intersection areas of cross-sections (as shown in Figure 14a). Although the most frequently selected patterns in the simulation with multiple variables also align with faults, the preferred areas in Figure 14b differ from those in Figure 14a. Additionally, Figure 14b reveals many regions with fewer selections, indicating a broader selection of patterns from cross-sections. The presence of fault contact relationships in the study area can lead to incorrect stratigraphic sequences during the simulation. To maintain a logical relationship among geological objects in the simulation results, the search and iteration processes in fault contact areas are more frequent than in other regions, increasing the likelihood of selecting patterns near faults. According to the relationship between the number of patterns and their corresponding selection frequencies shown in Figure 15, most patterns are selected fewer than 100 times. Figure 15 indicates a negative correlation between selection frequency and the number of patterns; as the frequency of selection increases, the number of corresponding patterns decreases sharply before stabilizing. Furthermore, in the simulation with multiple variables, the number of patterns selected only once exceeds 2,000,000, which is nearly four times greater than in the previous simulation.

Number distribution of patterns selected from cross-sections for the simulation by lithology (a) and multiple variables (b).

Two 3D simulation graphics labeled (a) and (b), show gridded structures with varying color patterns ranging from light blue to red, representing different values. The accompanying color scale on the right indicates values from 0 to 100. Axes labeled x, y, and z are displayed in green, red, and blue, respectively.

Relationship of the selected frequency and the number of patterns in simulation with lithology (in blue line) and with multiple variables (in red line).

Log-log graph showing patterns number versus numbers of selection. The blue line represents lithology and the red line represents multiple variants. Both lines decrease rapidly and fluctuate after a selection number of about 1.5 on the x-axis. The y-axis ranges from ten to ten million.

This study proposed an algorithm for integrating co-located multiple variables in the framework of MPS-based 3D geological modeling. Data fusion using a linear pooling method that incorporates entropy-based weighting values is applied during the process of selecting patterns for constructing an initial model, rather than by integrating multiple variables into a single mixed variable before modeling, as was done in the previous study. Here, the similarity of different variables is measured by cross-entropy values of the probability distribution of patterns, in which geological cross-sections are treated as the primary variable. Different variables are integrated using the linear pooling method, where the weighting values are calculated using the cross-entropy. Although profiles of different data types are co-located in this study, the integrating method is actually only related to the spatial distribution. Therefore, the proposed method is independent of the spatial positions of patterns. Considering profiles of seismic velocity, density, and lithology typically have different spatial resolutions and coverage. Thus, in further study, the concept of the proposed method can be applied to the integration of location-independent data.

Benefitting from their multi-layer network structure, deep learning (DL) methods, such as generative adversarial networks (GANs) and convolutional neural networks (CNNs), have successfully been used to construct 3D geological structures in different fields due to their strong image reproduction capabilities (Chan and Elsheikh, 2020; Laloy et al., 2017; Song et al., 2021; Yang et al., 2022). Theoretically, a CNN has a much better ability to extract and reconstruct spatial structures than the network used in this study. However, when only several cross-sections are used as TIs, it becomes a difficult task to capture enough patterns for training the GAN or CNN. In addition, overfitting is likely to occur when constructing geological structures from sparse data using CNNs (Hou et al., 2023). Here, an FCN based on the BP framework with a few parameters is used to simulate geological surfaces using the elevations of geological contacts, thereby avoiding the need for a large amount of training data. Because on each cross-section, the number of the extracted geological contacts is larger than the number of patterns that can be obtained from the cross-sections. In this study, subsurfaces of all geological objects are generated using the same FCN structures. Therefore, if a geological object is too small in the cross-section, the risk of overfitting in training would possibly appear. Furthermore, two FCNs are used for each block in the proposed method, resulting in high computational cost and challenges associated with multiple attributes per grid. Thus, in further research, an adaptive FCN framework should be introduced.

 Conclusion

This study presented a novel algorithm for integrating multiple variables within the framework of iterative MPS-based three-dimensional geological modeling. The data fusion process is carried out during the construction of the initial model, in which various variables are combined using a linear pooling method that incorporates entropy-based weighting values. It is posited that profiles of different data types are collocated. The lithology cross-section serves as the primary variable, while two additional profiles function as auxiliary variables. Two FCNs, characterized by a kernel function defined by contact elevation, are developed to generate surfaces corresponding to each geological object. Following the assignment of attributes to the grids, a stochastic process is used to construct a lithological model by integrating density and velocity data. An iterative EM-like algorithm is used to refine the initial model. The outcomes of reconstructing the overthrust model demonstrate that the geometry and interrelationships of faults and strata are effectively reconstructed, with enhanced precision achieved through the integration of multiple variables in the simulation.

The proposed algorithm overcomes the limitations associated with local optimization by generating surfaces that constrain initial modeling and using a multi-scale iterative strategy during simulation. The kernel function of this method synthesizes the loss functions of the FCN and traditional MPS simulation, constrained by geological relationships. Although the relationships among geological objects are adjusted in a post-processing phase, the absence of geological constraints during the iterative process may lead to the emergence of implausible scenarios in the final results. Consequently, the incorporation of geological constraints within the EM iterative process is critical for future research.

To save simulation time, the integration process is exclusively conducted during the initial model construction in this study. As a result, the outcomes of the EM-like iterative process may be influenced by the candidate patterns of lithology. Furthermore, the colocated fusion of multiple variables presents practical challenges in obtaining such data sources. Therefore, future research should focus on developing a novel method for integrating non-colocated multiple variables.

 Data availability statement

Publicly available datasets were analyzed in this study. These data can be found at https://zenodo.org/records/4252588.

 Author contributions

YC: Data curation, Resources, Supervision, Writing – review and editing. WH: Conceptualization, Data curation, Formal analysis, Funding acquisition, Validation, Visualization, Writing – original draft, Writing – review and editing. YL: Formal analysis, Methodology, Validation, Writing – original draft, Writing – review and editing. SY: Data curation, Formal analysis, Investigation, Writing – review and editing.

 Acknowledgements

The authors greatly appreciate the valuable comments of the reviewers in refining this manuscript.

 Conflict of interest

Authors YC and SY were employed by Guangzhou Metro Design & Research Institute Co. Ltd. Author YL was employed by Three Gorges Renewables Yangjiang Power Co., Ltd.

The remaining 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 Bai T. Tahmasebi P. (2020). Hybrid geological modeling: combining machine learning and multiple-point statistics. Comput. and Geosciences 142, 104519. 10.1016/j.cageo.2020.104519 Chan S. Elsheikh A. H. (2020). Parametrization of stochastic inputs using generative adversarial networks with application in geology. Front. Water 2, 5. 10.3389/frwa.2020.00005 Chen Q. Mariethoz G. Liu G. Comunian A. 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. 10.5194/hess-22-6547-2018 Comunian A. Renard P. Straubhaar J. (2012). 3D multiple-point statistics simulation using 2D training images. Comput. and Geosciences 40, 49–65. 10.1016/j.cageo.2011.07.009 Feng W. Wu S. Yin Y. Zhang J. Zhang K. (2017). A training image evaluation and selection method based on minimum data event distance for multiple-point geostatistics. Comput. and Geosciences 104, 35–53. 10.1016/j.cageo.2017.04.004 Feng W. Yin Y. Zhang C. Duan T. Zhang W. Houc G. (2019). A training image optimal selecting method based on composite correlation coefficient ranking for multiple-point geostatistics. J. Petroleum Sci. Eng. 179, 292–311. 10.1016/j.petrol.2019.04.046 Figueiredo L. P. d. Schmitz T. Lunelli R. Roisenberg M. Freitas D. S. d. Grana D. (2021). Direct multivariate simulation - a stepwise conditional transformation for multivariate geostatistical simulation. Comput. Geosciences 147, 104659. 10.1016/j.cageo.2020.104659 Guardiano F. B. Srivastava R. M. (1993). “Multivariate geostatistics: beyond bivariate moments,” in Geostatistics-Tróia’92 (Dordrecht: Kluwer Academic Publications), 133–144. Gueting N. Caers J. Comunian A. Vanderborght J. Englert A. (2018). Reconstruction of three-dimensional aquifer heterogeneity from two-dimensional geophysical data. Math. Geosci. 50 (1), 53–75. 10.1007/s11004-017-9694-x Guo J. T. Wang Z. X. Li C. L. Jessell M. Wu L. X. Wang J. M. (2022). Multiple-point geostatistics-based three-dimensional automatic geological modeling and uncertainty analysis for borehole data. Nat. Resour. Res. 31, 2347–2367. 10.1007/s11053-022-10071-6 Hansen T. M. Vu L. T. Mosegaard K. Cordua K. S. (2018). Multiple point statistical simulation using uncertain (soft) conditional data. Comput. and Geosciences 114, 1–10. 10.1016/j.cageo.2018.01.017 Hasanabadi A. Baniassadi M. Abrinia K. Safdari M. Garmestani H. (2019). Optimal combining of microstructures using statistical correlation functions. Int. J. Solids Struct. 160, 177–186. 10.1016/j.ijsolstr.2018.10.024 He X. L. Sonnenborg T. O. Jørgensen F. Jensen K. H. (2014). The effect of training image and secondary data integration with multiple-point geostatistics in groundwater modelling. Hydrology Earth Syst. Sci. 18, 2943–2954. 10.5194/hess-18-2943-2014 Hou W. Liu H. Zheng T. Chang H. Xiao F. (2022). Extended GOSIM: MPS-driven simulation of 3D geological structure using 2D cross-sections. Earth Space Sci. 9, e2021EA001801. 10.1029/2021ea001801 Hou W. Chen Y. Liu H. Xiao F. Liu C. Wang D. (2023). Reconstructing three-dimensional geological structures by the multiple-point statistics method coupled with a deep neural network: a case study of a metro station in guangzhou, China. TUST 136, 105089. 10.1016/j.tust.2023.105089 Hou W. Li Y. Ye S. Yang S. Xiao F. (2025). Mapping 3D overthrust structures by a hybrid modeling method. Earth Space Sci. 12, e2024EA003916. 10.1029/2024ea003916 Jóhannsson Ó. D. Hansen T. M. (2023). Multiple-point statistics and non-colocational soft data integration. Comput. Geosciences 172, 105280. 10.1016/j.cageo.2022.105280 Koch J. He X. Jensen K. H. Refsgaard J. C. (2014). Challenges in conditioning a stochastic geological model of a heterogeneous glacial aquifer to a comprehensive soft data set. Hydrology Earth Syst. Sci. 18, 2907–2923. 10.5194/hess-18-2907-2014 Laloy E. Hérault R. Lee J. Jacques D. Linde N. (2017). Inversion using a new low-dimensional representation of complex binary geological media based on a deep neural network. Adv. Water Resour. 110, 387–405. 10.1016/j.advwatres.2017.09.029 Liu Y. Harding A. Abriel W. Strebelle S. (2004). Multiple-point simulation integrating wells, three-dimensional seismic data, and geology. AAPG Bull. 88 (7), 905–921. 10.1306/02170403078 Mariethoz G. Caers J. (2014). Multiple-point geostatistics: stochastic modeling with training images. Wiley and Sons. Okabe H. Blunt M. J. (2007). Pore space reconstruction of vuggy carbonates using microtomography and multiple-point statistics. Water Resour. Res. 43, W12S02. 10.1029/2006wr005680 Pérez C. Mariethoz G. Ortiz J. M. (2014). Verifying the high-order consistency of training images with data for multiple-point geostatistics. Comput. and Geosciences 70, 190–205. 10.1016/j.cageo.2014.06.001 Shi C. Wang Y. (2022). Data-driven construction of three-dimensional subsurface geological models from limited site-specific boreholes and prior geological knowledge for underground digital twin. Tunn. Undergr. Sp. Tech. 126, 104493. 10.1016/j.tust.2022.104493 Song S. Mukerji T. Hou J. (2021). Geological facies modeling based on progressive growing of generative adversarial networks (GANs). Comput. Geosci. 25, 1251–1273. 10.1007/s10596-021-10059-w Wang L. Yin Y. Zhang C. Feng W. Li G. Chen Q. (2022). A MPS-based novel method of reconstructing 3D reservoir models from 2D images using seismic constraints. J. Petroleum Sci. Eng. 209, 109974. 10.1016/j.petrol.2021.109974 Wu Q. Xu H. Zou X. Lei H. (2015). A 3D modeling approach to complex faults with multi-source data. Comput. and Geosciences 77, 126–137. 10.1016/j.cageo.2014.10.008 Xiao K. Li N. Porwal A. Holden E.-J. Bagas L. Lu Y. (2015). GIS-based 3D prospectivity mapping: a case study of jiama copper-polymetallic deposit in Tibet, China. Ore Geol. Rev. 71, 611–632. 10.1016/j.oregeorev.2015.03.001 Yang L. Hou W. S. Cui C. J. Cui J. (2016). GOSIM: a multi-scale iterative multiple point statistics algorithm with global optimization. Comput. Geosci. 89, 57–70. 10.1016/j.cageo.2015.12.020 Yang Z. X. Chen Q. Y. Cui Z. S. Liu G. Dong S. Q. Tian Y. P. (2022). Automatic reconstruction method of 3D geological models based on deep convolutional generative adversarial networks. Comput. Geosci. 26, 1135–1150. 10.1007/s10596-022-10152-8 Yin Y. Xun H. Jixin H. Wenjie F. Linxin W. Taizhong D. (2022). A three-dimensional model of deep-water turbidity channel in plutonio oilfield, Angola: from training image generation, optimization to multi-point geostatistical modelling. J. Petroleum Sci. Eng. 195, 107650. 10.1016/j.petrol.2020.107650 Zhao Z. Zhou X.-P. (2019). An integrated method for 3D reconstruction model of porous geomaterials through 2D CT images. Comput. and Geosciences 123, 83–94. 10.1016/j.cageo.2018.11.012

Edited by: Yongzhang Zhou, Sun Yat-sen University Guangzhou, China

Reviewed by: Baoyi Zhang, Central South University, China

Bingli Liu, Chengdu University of Technology, China