Multiple-point Statistics Research Articles

Regionalized multiple-point statistical simulation for calibrating process-based geologic models to seismic data

Calibrating process-based geologic models to seismic data is critical and has been challenging for decades. The traditional approach to data calibration involves tuning the model input parameters by trial-and-error or through an automated inverse procedure. This can improve the model calibration to data but can hardly reach a fully satisfactory result. We adopt a multiple-point statistics (MPS) approach where a process-based geologic model is used as a training image for statistical pattern recognition. First, we define a rock-physics model from the process-based geologic model and derive its seismic attributes through seismic forward modeling. Then, we use the process-based model and its seismic attributes as coupled training images for geologic pattern recognition and regeneration under seismic data constraints. The method differs from the conventional MPS method in several ways: (1) the training image is a process-based geologic model of the reservoir of interest, thus defined on the same grid of the reservoir model; (2) the training image is generally nonstationary but there is no need to partition the nonstationary training image into pseudostationary ones; (3) the geologic facies and seismic constraint are related through seismic forward modeling instead of statistical inference, thus there is no need to convert seismic data to facies proportion or probability; (4) multiple-point statistics are based on Bayes’ law and Gaussian kernel approximation of conditional probability instead of a somehow arbitrary probability combination scheme or a heuristic rule; and (5) the method does not involve an iterative optimization procedure. Therefore, it also differs from the neural-network-based machine-learning approach, where the data conditioning is achieved through an iterative optimization procedure. These differences make our method advantageous for calibrating process-based geologic models. The two examples with synthetic data illustrate the effectiveness of the method. The first one is derived from a satellite image of a river delta, and the second one from the Stanford VI reservoir model.

Pore-scale numerical investigation on spontaneous imbibition in natural fracture with heterogeneous wettability using the volume of fluid method

Understanding the spontaneous imbibition in the natural fracture with heterogeneous wettability is crucial for predicting and mitigating the impacts of unstable displacement on unconventional recovery. In this paper, the fracture structured mesh model is reconstructed based on the micro-computed tomography (micro-CT) image of naturally fractured tight sandstone. The mineralogy map-based modeling method for heterogeneous-wetting fracture is developed by combining the thin section images, X-ray diffraction (XRD) analysis, and multiple point statistics method. The simulation of the single-phase flow is performed to test the mesh independence. The effects of gravity and wettability on spontaneous imbibition in natural fracture and corresponding imbibition front dynamics are analyzed and discussed using the volume of fluid (VOF) method. The results show that (1) The structured mesh reconstruction method proposed in this paper can more effectively preserve the fracture structure compared to the unstructured mesh reconstruction method. (2) Gravity has a negligible impact on the pore-scale spontaneous imbibition in natural fracture. Under homogeneous-wetting conditions, spontaneous imbibition in natural fracture consistently exhibits stable displacement without significant residual gas formation. However, under the heterogeneous-wetting condition, the spontaneous imbibition displays typical capillary fingering, resulting in approximately 24.04% of the gas being trapped after spontaneous imbibition. The residual gas trapping mechanisms mainly include adhered, isolated, and connected gas. (3) Under both homogeneous- and heterogeneous-wetting conditions, the imbibing water saturation and the length of the imbibition front are proportional to the power of imbibition time during spontaneous imbibition in the natural fracture.

Using multiple-point geostatistics for geomodeling of a vein-type gold deposit

Geostatistical cascade modeling of Mineral Resources is challenging in vein-type gold deposits. The narrow shape and long-range features of these auriferous veins, coupled with the paucity of drill-hole data, can complicate the modeling process and make the use of two-point geostatistical algorithms impractical. Instead, multiple-point geostatistics techniques can be a suitable alternative. However, the most challenging part in implementing the MPS is to use a suitable training data set or training image (TI). In this paper, we suggest using the radial basis function algorithm to build a training image and the DeeSse algorithm, one of the multiple-point statistics (MPS) methods, to model two long-range veins in a gold deposit. It is demonstrated that DeeSse can replicate long-range vein features better than plurigaussian simulation techniques when there is a lack of conditioning data. This is shown by several validation processes, such as comparing simulation results with an interpretive geological block model and replicating geological proportions.

Borehole‐Based Interval Kriging for 3D Lithofacies Modeling

AbstractDeveloping a three‐dimensional (3D) lithofacies model from boreholes is critical for providing a coherent understanding of complex subsurface geology, which is essential for groundwater studies. This study aims to introduce a new geostatistical method—interval kriging—to efficiently conduct 3D borehole‐based lithological modeling with sand/non‐sand binary indicators. Interval kriging is a best linear unbiased estimator for irregular interval supports. Interval kriging considers 3D anisotropies between two orthogonal components—a horizontal plane and a vertical axis. A new 3D interval semivariogram is developed. To cope with the nonconvexity of estimation variance, the minimization of estimation variance is regulated with an additional regularization term. The minimization problem is solved by a global‐local genetic algorithm embedded with quadratic programming and Brent's method to obtain kriging weights and kriging length. Four numerical and real‐world case studies demonstrate that interval kriging is more computationally efficient than 3D kriging because the covariance matrix is largely reduced without sacrificing borehole data. Moreover, interval kriging produces more realistic geologic characteristics than 2.5D kriging, while conditional to spatial borehole data. Compared to the multiple‐point statistics (MPS) algorithm—SNESIM, interval kriging can reproduce the geological architecture and spatial connectivity of channel‐type features, meanwhile producing tabular‐type features with better connectivity. Because the regularization term constrains kriged value toward 0 or 1, interval kriging produces more certainty in sand/non‐sand classification than 2.5D kriging, 3D kriging, and SNESIM. In conclusion, interval kriging is an effective and efficient 3D geostatistical algorithm that can capture the 3D structural complexity while significantly reducing computational time.

Using Pattern Counts to Quantify the Difference Between a Pair of Three-Dimensional Realizations

AbstractWhen comparing different ways of modeling discrete three-dimensional realizations such as facies, it is useful to have a measure of difference (or similarity) in the geometry of these realizations. We propose a method for evaluating such difference by comparing pattern counts for a small template. Tests on synthetic datasets demonstrate that the proposed difference effectively differentiates between realizations of a Boolean model and those generated using multiple-point statistics with the Boolean realizations as training images. We also observed that multiple-point statistics realizations based on similar training images yield smaller differences to one another compared to those based on training images from dissimilar concepts. This suggests that the proposed difference is a useful tool for comparing discrete three-dimensional realizations.

Fracture density reconstruction using direct sampling multiple-point statistics and extreme value theory

The aim of this work is to present a methodology for the reconstruction of missing fracture density within highly fractured intervals, which can represent preferential fluid flow pathways. The lack of record can be very common due to the intense presence of fractures, dissolution processes, or data acquisition issues. The superposition of numerous fractures makes the definition of fracture surfaces impossible, as a consequence, modeling such zones is challenging. In order to address this issue, the usage of direct sampling multiple-point statistics to perform gap filling in well logs is demonstrated as an alternative to other techniques. It reproduces data patterns and provides several models representing uncertainty. The method was tested in intervals from a highly fractured well, by removing previously known fracture density data, and simulating different scenarios with direct sampling. Simulation results are compared to the observed data using cross-validation and continuous rank probability score. The reference scenario training data set consists in one well and two variables: fracture density and fracture occurrence. A sensitivity analysis is carried out considering additional variables, additional wells, different intervals, resampling with extremes, and other gap filling techniques. The auxiliary variable plays an important role in pattern matching, but adding wells and logs increases the complexity of the method without improving pattern retrieval. Best results are obtained applying extreme values theory for stochastic process with the enrichment of the fracture density data at the tail region, followed by resampling of the new values. The enriched data is used for the gap filling resulting in lower continuous rank probability score, and the achievement of extreme fracture density values.

Comparison of three recent discrete stochastic inversion methods and influence of the prior choice

Groundwater flow depends on subsurface heterogeneity, which often calls for categorical fields to represent different geological facies. The knowledge about subsurface is however limited and often provided indirectly by state variables, such as hydraulic heads of contaminant concentrations. In such cases, solving a categorical inverse problem is an important step in subsurface modeling. In this work, we present and compare three recent inverse frameworks: Posterior Population Expansion (PoPEx), Ensemble Smoother with Multiple Data Assimilation (ESMDA), and DREAM-ZS (a Markov chain Monte Carlo sampler). PoPEx and ESDMA are used with Multiple-point statistics (MPS) as geostatistical engines, and DREAM-ZS is used with a Wasserstein generative adversarial network (WGAN). The three inversion methods are tested on a synthetic example of a pumping test in a fluvial channelized aquifer. Moreover, the inverse problem is solved three times with each method, each time using a different training image to check the performance of the methods with different geological priors. To assess the quality of the results, we propose a framework based on continuous ranked probability score (CRPS), which compares single true values with predictive distributions. All methods performed well when using the training image used to create the reference, but their performances were degraded with the alternative training images. PoPEx produced the least geological artifacts but presented a rather slow convergence. ESMDA showed initially a very fast convergence which reaches a plateau, contrary to the remaining methods. DREAM-ZS was overly confident in placing some incorrect geological features but outperformed the other methods in terms of convergence.

Stochastic reconstruction of geological reservoir models based on a concurrent multi-stage U-Net generative adversarial network

For complex geological reservoir modeling, some numerical-simulation-based methods, such as traditional multiple-point statistics (MPS), cannot extract nonstationary patterns of training images (TIs) effectively. CPU-intensive calculations require that variables information only be stored in RAM instead of other storage mediums to avoid huge time consumption if many simulations are performed successively. Generative Adversarial Network (GAN) modeling methods are able to overcome these issues, but when the training datasets are limited, mode collapse easily occurs. Meanwhile, fixed receptive fields might cause the problem that the corresponding reconstruction cannot take local features and global features of the TIs into account simultaneously. Therefore, we propose the Concurrent Multi-Stage U-Net GAN (Con-MSUGAN) based on a single TI. The pyramid structure is used for spatial pattern multi-scale representation in Con-MSUGAN, which can capture the corresponding feature information of TIs under different scales. Meanwhile, concurrent training patterns can make network parameters in adjacent stages correlated with each other to accelerate network convergence speed, thereby realizing stochastic multi-scale reconstructions. Stationary channel and nonstationary delta facies TIs were used to verify the performance of Con-MSUGAN. Results show that simulation of different datasets using this method are more similar in terms of spatial variance, connectivity degree and facies type frequency distribution to the original TIs. Meanwhile, muti-dimensional-scaling (MDS) plots show small discrepancies between simulation results and the corresponding TI. It indicates that Con-MSUGAN can overcome the problem that complex spatial patterns are hard to reproduce. In addition, saved network parameters can be reused. A variety of equiprobability simulation results can be obtained rapidly, thereby realizing stochastic reconstruction of geological reservoir models.

Extraction of weak geochemical anomalies based on multiple-point statistics and local singularity analysis

Extraction of weak geochemical anomalies based on multiple-point statistics and local singularity analysis

Stochastic Reconstruction of 3D Heterogeneous Microstructure Using a Column-Oriented Multiple-Point Statistics Program

Abstract Three-dimensional (3D) microstructure reconstruction is a key approach to exploring the relationship between pore characteristics and physical properties. Viewing the training image as a prior model, multiple-point statistics (MPS) focus on reproducing spatial patterns in the simulation grid. However, it is challenging to efficiently generate 3D nonstationary models with varying microstructures. In this work, we propose column-oriented simulation (ColSIM) to achieve the stochastic reconstruction of 3D porous media. A heterogeneous system is understood as a spatially evolving process that consists of frequent transitions of small magnitude and abrupt changes of large magnitude. First, a training image selection step is suggested to find representative microstructures. Our program applies modified Hausdorff distance, t-distributed stochastic neighboring embedding, and spectral clustering to organize two-dimensional (2D) candidate images. The medoid of each group is applied to guide the following programs. Second, we introduce column-oriented searching into MPS. To save simulation time, a subset of conditioning points is checked to find desired instances. Our program suggests an early stopping strategy to address complex microstructures. Third, a contrastive loss term is designed to create 3D models from 2D slice. To automatically calibrate the volume fraction and simplify parameter specification, the computer consistently monitors the difference between the present model and the target. The performance of ColSIM is examined by 3D multiphase material modeling and 3D heterogeneous shale simulation. To achieve quantitative evaluation, we compute various statistical functions and physical descriptors on simulated realizations. The proposed ColSIM exhibits competitive performance in terms of calculation efficiency, microstructure reproduction, and spatial uncertainty.

Knowledge-based multiple point statistics for soil stratigraphy simulation

An accurate representation of subsurface stratigraphy in the form of geotechnical profiles is a prerequisite for optimal underground design and safe construction. This paper proposes a novel knowledge-based multiple point statistics (KMPS) method to simulate the heterogeneities and spatial trends of subsurface strata and to generate geotechnical profiles. This approach incorporates domain knowledge in the form of cone penetration test and standard penetration test data as well as results from laboratory testing to increase simulation accuracy. A novel varying search template size and similarity threshold is adopted to alleviate the computational cost of the simulations. Ablation analyses are used to quantify the improvement of the proposed KMPS algorithm over existing baseline multiple point statistics algorithms. By way of example, the proposed approach is applied to real-world exploration data from the Suzhou No.6 metroline project. The results show that the computational cost of KMPS is effectively minimized whilst the simulation accuracy can be improved by up to 13.3%. The algorithm is further validated by extending the simulation to a test image. With sparse exploration data, the KMPS method is shown to effectively and efficiently predict soil profiles based on prior experience embodied by a training image, thus enabling a reduction in the number of field exploration boreholes.

Generative Adversarial Networks to incorporate the Training Image uncertainty in multiple-point statistics simulation

Multiple-point geostatistics (MPS) algorithms have been extensively used for geostatistical reservoir modeling. These techniques are known for building geostatistical simulations with curvilinear features, which are crucial for applications in flow simulation. One challenge in using MPS algorithms is obtaining a representative Training Image (TI). The result is that often one TI is used to condition the entire set of geostatistical simulations, and the TI uncertainty is disregarded. If the TI is uncertain, this uncertainty should be considered in the geostatistical simulations. Otherwise, the spatial uncertainty represented by the geostatistical simulations is underestimated. Incorporating the TI uncertainty requires building a data set or catalog of plausible TIs. In this context, we propose to use Generative Adversarial Networks (GANs) to build this catalog of TIs prior to the MPS simulation. Each TI generated by GAN is used to condition one realization built by the chosen MPS algorithm. The proposed workflow combines the strength of GAN to generate TIs with the strength of MPS algorithms to consider several sources of information to condition the simulations. A case study is shown, where the methodology is compared against the traditional MPS workflow, which uses a single TI. The results show that the proposed method successfully incorporated the TI uncertainty into the geostatistical workflow. The proposed method built numerical models with higher uncertainty and variability and, thus, provides a more realistic quantification of the space of uncertainty.

Multiple-point statistical modeling of three-dimensional glacial aquifer heterogeneity for improved groundwater management

Quaternary glacial aquifers are important water sources for irrigation in many agricultural regions, including eastern Nebraska, USA. Quaternary glacial aquifers are heterogeneous, with juxtaposed low-permeability and high-permeability hydrofacies. Managing groundwater in such aquifers requires a realistic groundwater-flow model parameterization, and characterization of the aquifer geometry, spatial distribution of aquifer properties, and local aquifer interconnectedness. Despite its importance in considering uncertainty during decision-making, hydrofacies probabilities generated from multiple-point statistics (MPS) are not widely applied for groundwater model parameterization and groundwater management zone delineation. This study used a combination of soft data, a cognitive training image, and hard data to generate 100 three-dimensional (3D) conditional aquifer heterogeneity realizations. The most probable model (probability of hydrofacies) was then computed at node spacing of 200 × 200 × 3 m and validated using groundwater-level hydrographs. The resulting hydrofacies probability grids revealed variations in aquifer geometry, locally disconnected aquifer systems, recharge pathways, and hydrologic barriers. The profiles from hydrofacies probability at various locations show spatial variability of the streambed and aquifer connectivity. Groundwater-level hydrographs show evidence of these aquifer characteristics, verifying the general structure of the model. Using the MPS-generated 3D hydrofacies probability and hydrologic data, a novel workflow was developed in order to better define high-resolution groundwater management zones and strategies. In general, the conditional probability of hydrofacies helps improve the understanding of glacial aquifer heterogeneity, the characterization of aquifer-to-aquifer and streambed-aquifer connections, and the delineation of groundwater management zones. This MPS workflow can be adapted to other areas for modeling 3D aquifer heterogeneity using multisource data.

N-Map: High-resolution groundwater N-retention mapping and modelling by integration of geophysical, geological, geochemical, and hydrological data

A key aspect of protecting aquatic ecosystems from agricultural nitrogen (N) is to locate (i) farmlands where nitrate leaches from the bottom of the root zone and (ii) denitrifying zones in the aquifers where nitrate is removed before entering the surface water (N-retention). N-retention affects the choice of field mitigation measures to reduce delivered N to surface water. Farmland parcels associated with high N-retention gives the lowest impact of the targeted field measures and vice versa. In Denmark, a targeted N-regulation approach is currently implemented on small catchment scale (approx. 15 km2). Although this regulatory scale is much more detailed than what has been used previously, it is still so large that regulation for most individual fields will be either over- or under-regulated due to large spatial variation in the N-retention. The potential cost reduction for farmers is of up to 20–30% from detailed retention mapping at the field scale compared to the current small catchment scale.In this study, we present a mapping framework (N-Map) for differentiating farmland according to their N-retention, which can be used for improving the effectiveness of targeted N-regulation. The framework currently only includes N-retention in the groundwater. The framework benefits from the incorporation of innovative geophysics in hydrogeological and geochemical mapping and modelling. To capture and describe relevant uncertainties a large number of equally probable realizations are created through Multiple Point Statistical (MPS) methods. This allows relevant descriptions of uncertainties of parts of the model structure and includes other relevant uncertainty measures that affects the obtained N-retention. The output is data-driven high-resolution groundwater N-retention maps, to be used by the individual farmers to manage their cropping systems due to the given regulatory boundary conditions.The detailed mapping allows farmers to use this information in the farm planning in order to optimize the use of field measures to reduce delivered agricultural N to the surface water and thereby lower the costs of the field measures. From farmer interviews, however, it is clear that not all farms will have an economic gain from the detailed mapping as the mapping costs will exceed the potential economic gains for the farmers. The costs of N-Map is here estimated to 5–7 €/ha/year plus implementation costs at the farm. At the society level, the N-retention maps allow authorities to point out opportunities for a more targeted implementation of field measures to efficiently reduce the delivered N-load to surface waters.

Markov chain Monte Carlo for seismic facies classification

Seismic facies classification aims to predict a facies model, or a set of facies models, from measured seismic data. We focus on stochastic classification methods to estimate the probability distribution of facies conditioned on seismic data. Bayesian classification methods based on analytical solutions are generally applied to seismically inverted elastic attributes or petrophysical properties rather than measured seismic data, and they do not include spatial correlation models in the prior distribution. Instead, iterative stochastic methods, such as Monte Carlo acceptance-rejection sampling and approximate Bayesian computation, can be directly applied to seismic data and can be implemented with complex prior models; however, their convergence is generally slow because it requires the simulation of a large number of facies models from the prior distribution. We develop an efficient implementation of a Markov chain Monte Carlo (MCMC) approach that adopts complex prior models, such as multiple-point statistics simulations based on a training image, to generate geologically realistic facies realizations. The novelty of the approach is that the proposal distribution of the proposed MCMC method is based on a gradual deformation algorithm, in which the proposal distribution is expressed as a linear combination of the indicator variables of the previously accepted model and a random prior simulation, according to a uniform random weight. The synthetic examples indicate accurate facies predictions with a success rate of approximately 80% for well-log facies classification based on elastic attributes and approximately 60% for seismic facies classification based on seismic data. The inversion is compared to an MCMC method with prior models sampled from a first-order Markov chain and Bayesian facies classification.

Machine Learning-Based Urban Renovation Design for Improving Wind Environment: A Case Study in Xi’an, China

The high-density urban form and building arrangement of modern cities have contributed to numerous environmental problems. The calm wind area caused by inappropriate building arrangements results in pollutant accumulation. To realize a practical design and improve urban microclimate, we investigated the spatial relationship between roads, buildings, and open space using the machine learning technique. First, region growing and k-means clustering were employed to identify roads and buildings. Based on the image masking program, we selected training areas according to the land use map. Second, we used the multiple-point statistics technique to create new urban fabric images. Viewing the training image as a prior model, our program constantly reproduced morphological structures in the target area. We intensified the similarity with training areas and enriched the variability among generated images. Third, Hausdorff distance and multidimensional scaling were applied to achieve a quality examination. The proposed method was performed to fulfill an urban renovation design in Xi’an, China. Based on the historical record, we applied computational fluid dynamics to simulate air circulation and ventilation. The results indicate that the size of calm wind area is reduced. The wind environment is significantly improved due to the rising wind speed.

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

A credible large-scale three-dimensional (3D) geological model constructed from a small amount of local high-density data provides a vital data infrastructure for the design and construction of metro lines. The multiple-point statistics (MPS) method has successfully generated 3D geological models in many fields. However, the global cognition of the spatial relations between geological objects is difficult to reconstruct using a local optimization strategy with moving templates in the MPS simulation. This study presented an MPS simulation approach for constructing 3D geological structures with two-dimensional (2D) cross-sections, where the MPS process was coupled with a fully connected deep neural network (FCDNN) simulation. The core idea of the proposed method was that the global features of each geological object were generated from the trained FCDNN, the input parameter of which was the subsurface depth. Then, the entire model was constructed on the initial model using the following iterative MPS process: the semantic relationship of geological objects was extracted from 2D cross-sections and used as a constraint in the modeling process, which ensured the order of the stratigraphic sequence and rationality of the geological structure in the final model. The joint objective function of the proposed method was composed of the loss function of the FCDNN and conventional MPS objective function constrained by the spatial semantics of geological objects. The concrete example in the Jiangtai Road Metro Station of Guangzhou Line 11, China illustrated that the details of the complex geological structures could be obtained accurately by the presented method. The global features from the trained FCDNN provided the basic geometry for the final model, which overcame the limitations of the conventional MPS method.

Development of training image database for subsurface stratigraphy

ABSTRACT Image-based stochastic simulation methods, such as multiple point statistics (MPS), can be viewed as a physics-informed Bayesian learning approach, which samples typical stratigraphic patterns from a single training image for onward conditional modelling of subsurface stratigraphy. A training image is essentially a prior geological model, which comprises representative stratigraphic connectivity at the site of interest. One key difficulty hindering wide application of image-based geological modelling methods is the lack of qualified training images. In this study, a systematic framework is proposed to develop training image databases for conditional simulations of subsurface stratigraphy. Collected training images can be further categorised based on three key factors, namely, geological origin, site location and application scenario. As a pilot study, a total of 54 geological cross-sections, mainly interpreted by experienced engineering practitioners, for weathered granite and tuff slopes in Hong Kong are collected and compiled as two training image databases. To demonstrate value and application of the established training image databases, subsurface stratigraphy for real weathered granite slope examples are used as illustrative examples, and stratigraphic uncertainty is also quantified. Results indicate that training image databases are of great significance for subsurface stratigraphy and uncertainty quantification, particularly when only limited site-specific data are available.

Multiple-point statistics and non-colocational soft data integration

In geostatistics, available conditional information is typically categorized as either hard (with no uncertainty) or soft (associated with an uncertainty) data. 2-point based Gaussian (kriging) simulation methods have since their inception been able to account for hard and soft data. In contrast, multiple-point statistical (MPS) simulation methods, that allow quantifying more complex phenomena than Gaussian methods, have focused mostly on conditioning to hard data. The most widely used method for accounting for soft data is to consider only co-located soft data. Here we define conditional information data through a probability distribution and show how a general probabilistic solution, in the form of a posterior probability distribution, can be obtained, combining conditional categorical information with information from a categorical training image. We demonstrate how to compute the conditional distribution, as needed by most MPS-based sequential simulation algorithms, that allows accounting for both co- and non-co-located hard and soft data. Examples are provided comparing the direct computation of conditional statistics to using more computational demanding rejection sampling. In order to reduce the significant computational demands, a distance-based measure is proposed to increase the entropy of conditional data, based on the distance to the point being estimated and simulated.

Numerical modelling of reservoir at pore scale: A comprehensive review

Pore scale reservoir modelling is carried out to digitally build the inner structure of the reservoir samples using proper mathematical algorithms based on certain priors extracted from training image(s) or laboratory measured data. Although various imaging devices can provide the inner structure of the rock sample directly, numerical modelling methods still has its irreplaceable advantages. Firstly, due to some intrinsic limitations, no imaging device so far can provide a 3D porous structure satisfying both resolution and field of view except for a very few of homogeneous reservoir samples with relatively large pore-throat size, while numerical modelling is an effective approach to improve the resolution of the machine produced images. Secondly, in some cases, where only 2D images (such as the images of drill cuttings) are available, numerical modelling is the only option to reconstruct a 3D porous structure. Thirdly, numerical modelling is a convenient and low-cost technique that makes it widely applicable. This paper presents a comprehensive review of the pore scale reservoir modelling. Pore structure's modelling can be treated as a supervised prediction process with two steps: 1) extraction of descriptors (prior knowledge) from the training image(s) or laboratory measured data; and 2) numerical reconstruction under the supervision of the extracted descriptors. In terms of extraction of descriptors, using multiple-point statistics (MPS) instead of random function models is a great step that made the empirical multivariate distributions inferred from training images can be directly used to reconstruct porous structure. Another great progress is the application of multiple-resolution training images where the image degradation mechanism contained in these different resolution images can be used to supervise the reconstruction. In terms of supervised reconstruction, some mathematical strategies such as Gaussian Random field, simulated annealing are applied to reproduce these descriptors in the reconstructed images. Along with the flourish of deep learning, a convolutional neural networks model is introduced into reconstruction process with an impressive performance.