Geostatistical Simulation Methods Research Articles

A novel geostatistical index of uncertainty for short-term mining plan

ABSTRACT The uncertainty associated with geostatistical estimates is as important in a model as the estimated value itself. This uncertainty is commonly assessed through different scenarios obtained by geostatistical simulation methods. Although the applications of geostatistical simulation methods have grown significantly in recent decades in long-term modeling, the same improvement has not occurred in medium- and short-term modeling. This study presents the “geological-operational risk indicator,” which comprises: (1) the uncertainty of geological contacts, determined by the variance of the sequential indicator simulation; and (2) the continuity of grades, obtained by sequential Gaussian simulation for variability of grades in relation to the economic or operational limit of that variable. The proposed index is applied to five grade variables of a short-term model of the Capitão do Mato iron mine in Brazil. The results are satisfactory and advance the concept of a geological risk indicator that is appropriate for medium- and short-term operations. This novel indicator helps to predict and manage in advance the uncertainties associated with mining plans.

Stochastic simulation of facies using deep convolutional generative adversarial network and image quilting

Sedimentary facies simulation is one of the essential works in sedimentary environment analysis and reservoir characterization. The traditional facies simulation method is based on geostatistics. However, the traditional two-point geostatistics-based facies simulation method cannot characterize complex facies structures. Most multiple-point geostatistical simulation methods are unable to flexibly generate abundant geologic patterns. To address these shortcomings, we develop an intelligent method to automatically simulate sedimentary facies according to the training image provided by geologists. The method can learn efficient representations for complex facies architectures and obtain the simulation results in a larger area, rather than being limited to an area with the same size as the training image. First, we construct a deep convolutional generative adversarial network to extract the high-dimensional features of facies. Then, a large number of specific patterns are randomly generated based on these features. Thus, the diversity of geologic patterns is improved. Finally, the patterns are spliced together to obtain possible facies maps by using an improved image quilting algorithm. A model test is described and analyzed to demonstrate the effectiveness and reliability of the new method. The results are consistent with the actual situation in the aspect of variability and continuity. The method is also applied to non-stationary geological facies unconditional simulation. The successful application indicates that the method is able to learn the features of non-stationary geological phenomena, showing the practicability of the proposed methods.

Quantitative 3-D Model of Carbonate Reef and Shoal Facies Based on UAV Oblique Photogrammetry Data: A Case Study of the Jiantanba Outcrop in West Hubei, China

Aiming at the problem of insufficient data acquisition for steep carbonate outcrops, we used unmanned aerial vehicle oblique photogrammetry to quantitatively evaluate Jiantianba outcrops from a global perspective. This method can quickly and flexibly acquire, process and interpret outcrop data, establish the three-dimensional digital outcrop model and quantitative reservoir geological knowledge database. Through the fine analysis and comparative study of multiple outcrop sections, we use lithofacies combination, hierarchical interface, sequential indicator stochastic simulation and multiple-point geostatistical simulation methods to establish a three-dimensional multi-point statistical model of the outcrop area. The results show that the model of the Jiantianba carbonate outcrop has three-dimensional coordinates that correspond to oblique photograph image data, allowing for both the accurate location of any point and the measurement of the rock body, and thus providing a base for studying the stratigraphic architecture of the outcrop. Through fine-scale geological and statistical analyses of the geological parameters of three geologic sections, we established a database and model of a typical outcrop, including the types, continuity and scale of lithology, as well as different lithological combinations and the geological evolution of the region. And we established a geological model and quantitative geological database via sequential indicator simulation, using virtual well location, lithological combination, stratigraphic correlation, three-dimensional stratigraphic and lithological models of a typical outcrop. Subsequently, we analyzed the training images of different microfacies based on these models combined with the outcrop geometry. Finally, we established a three-dimensional geological model of reef–shoal facies that is more in line with our current understanding of the geology of the outcrop area by multi-point geostatistical simulations.

Geostatistical modeling of porosity and evaluating the local and global distribution

Porosity is one of the main variables needed for reservoir characterization. For this volumetric variable, there are many methods to simulate the spatial distribution. In this article, porosity was analyzed and modeled in the local and global distribution. For simulation, Sequential Gaussian simulation (SGS) and Gaussian Random Function (GRFS) were applied. Also, kriging was used to estimate the porosity at specific locations. The main purpose of this work was to investigate the porosity to compare geostatistical simulation and estimation methods in a sandstone reservoir as a real case study. First, the data sets were normalized by the Normal Scores Transformation (NST) and stratigraphic coordinate. The model of experimental variograms was fitted in the vertical and horizontal directions. For the simulation methods, 10 realizations were generated by each method. The Q-Q plots were calculated, and both sets of quintiles (Target Porosity Distribution versus Porosity realization) came from normal distributions with the following correlation coefficients: 0.93, 0.94 and 0.97 related to GRFS, SGS and Kriging, respectively. The extracted variograms from realizations showed that the kriging couldn’t reproduce the variograms with global distribution. For local validation, the cross-validation was evaluated and three wells were omitted. The re-estimation of porosity was considered at located well logs through the well sections window where the kriging had a better performance with minimum error to estimate porosity locally. Finally, the cross-sectional models were generated by each algorithm which showed that the simple kriging tries to produce smoother distribution, whereas conditional simulations (SGS and GRFS) try to represent more global-detailed sections.

Hybrid parallel framework for multiple-point geostatistics on Tianhe-2: A robust solution for large-scale simulation

Multiple-point geostatistical (MPS) simulation methods have attracted an enormous amount of attention in earth and environmental sciences due to their ability to enhance extraction and synthesis of heterogeneous patterns. To characterize the subtle features of heterogeneous structures and phenomena, large-scale and high-resolution simulations are often required. Accordingly, the size of simulation grids has increased dramatically. Since MPS is a sequential process for each grid unit along a simulation path, it results in severe computational consumption. In this work, a new hybrid parallel framework is proposed for the case of MPS simulation on large areas with enormous amount of grid cells. Both inter-node-level and intra-node-level parallel strategies are combined in this framework. To maintain the quality of the realizations, we implement a conflict control method adapting to the Monte-Carlo process. Also, an optimization method for the simulation information is embedded to reduce the inter-node communication overhead. A series of synthetic tests were used to verify the availability and performance of the proposed hybrid parallel framework. The results corroborate that the proposed framework can efficiently achieve the high-resolution reproduction and characterization of complex structures and phenomena in earth sciences.

U-net generative adversarial network for subsurface facies modeling

Subsurface models are central pieces of information in different earth-related disciplines such as groundwater management and hydrocarbon reservoir characterization. These models are normally obtained using geostatistical simulation methods. Recently, methods based on deep learning algorithms have been applied as subsurface model generators. However, there are still challenges on how to include conditioning data and ensure model variability within a set of realizations. We illustrate the potential of Generative Adversarial Networks (GANs) to create unconditional and conditional facies models. Based on a synthetic facies dataset, we first train a Deep Convolution GAN (DCGAN) to produce unconditional facies models. Then, we show how image-to-image translation based on a U-Net GAN framework, including noise-layers, content loss function and diversity loss function, is used to model conditioning geological facies. Results show that GANs are powerful models to capture complex geological facies patterns and to generate facies realizations indistinguishable from the ones comprising the training dataset. The U-Net GAN framework performs well in providing variable models while honoring conditioning data in several scenarios. The results shown herein are expected to spark a new generation of methods for subsurface geological facies with fragmentary measurements.

SPS Model: a significant algorithm to reduce the time and computer memory required in geostatistical simulations

In geochemical anomaly classification, different mathematical-statistical models have been applied. The final classified map provides only one scenario. This model is not certain enough since every model provides several thresholds which are almost different from each other meaning dissimilarity and spatial uncertainty of the classified maps. Spatial uncertainty of the models could be quantified considering the difference between the associated geochemical scenarios simulated (called: ‘realizations’) by geostatistical simulation (GS) methods. However, the main problem with GS methods is that these methods are significantly time-consuming, and CPU- and memory-demanding. To improve such problems, in this research, the method of “scaling and projecting sample-locations (SPS)” is developed. Based on the SPS theory, first of all, the whole sample-locations were projected (centralized) and scaled into a box coordinated between (0,0) to (150,0) and (0,0) to (0,100), for example (they can be equal though), with the cell-size of 1 m2. Therefore, the time consumed and the memory demanded to generate a large number of realizations, for example, 1000 realizations based on the non-scaled/non-projected (NS/NP) and scaled/projected (S/P) sample locations per case-study were quantified. In this study, the turning bands simulation (TBSIM) were applied to geochemical datasets of three different case studies to take the area scales, regularity/irregularity and density of the samples into account. The comparison between NS/NP and S/P results statistically demonstrated the same results, however, the process and outputs of the S/P samples took a significantly shorter time and consumed a remarkably lower computer-memory. Therefore, experts are able to easily run this algorithm using any normal computer.

Stochastic modeling of subglacial topography exposes uncertainty in water routing at Jakobshavn Glacier

AbstractSubglacial topography is an important feature in numerous ice-sheet analyses and can drive the routing of water at the bed. Bed topography is primarily measured with ice-penetrating radar. Significant gaps, however, remain in data coverage that require interpolation. Topographic interpolations are typically made with kriging, as well as with mass conservation, where ice flow dynamics are used to constrain bed geometry. However, these techniques generate bed topography that is unrealistically smooth at small scales, which biases subglacial water flowpath models and makes it difficult to rigorously quantify uncertainty in subglacial drainage patterns. To address this challenge, we adapt a geostatistical simulation method with probabilistic modeling to stochastically simulate bed topography such that the interpolated topography retains the spatial statistics of the ice-penetrating radar data. We use this method to simulate subglacial topography using mass conservation topography as a secondary constraint. We apply a water routing model to each of these realizations. Our results show that many of the flowpaths significantly change with each topographic realization, demonstrating that geostatistical simulation can be useful for assessing confidence in subglacial flowpaths.

Direct simulation of non-additive properties on unstructured grids

Uncertainties related to permeability heterogeneity can be estimated using geostatistical simulation methods. Usually, these methods are applied on regular grids with cells of constant size, whereas unstructured grids are more flexible to honor geological structures and offer local refinements for fluid-flow simulations. However, cells of different sizes require to account for the support dependency of permeability statistics (support effect).This paper presents a novel workflow based on the power averaging technique. The averaging exponent ω is estimated using a response surface calibrated from numerical upscaling experiments. Using spectral turning bands, permeability is simulated on points in each unstructured cell, and later averaged with a local value of ω that depends on the cell size and shape.The method is illustrated on a synthetic case. The simulation of a tracer experiment is used to compare this novel geostatistical simulation method with a conventional approach based on a fine scale Cartesian grid. The results show the consistency of both the simulated permeability fields and the tracer breakthrough curves. The computational cost is much lower than the conventional approach based on a pressure-solver upscaling.

Evaluation of liquefaction probability of earth-fill dam over next 50 years using geostatistical method based on CPT

A method for evaluating the liquefaction probability of an earth-fill dam over the next 50 years is presented through the use of a geostatistical method for the measured values from cone penetration tests (CPTs). In particular, this paper discusses a new procedure for evaluating the liquefaction probability based on CPTs. Although the fines content, Fc, and the N-value are required in the Japanese standards to evaluate the liquefaction risk, the number of test data is not enough for the statistical modeling. Herein, Fc and the N-value are derived directly from CPTs. The statistical modeling procedure for Fc and the N-value is the unique point of this study. Since CPTs can be conducted with short intervals, especially in the horizontal direction, the geostatistical parameters can be determined, and the geostatistical simulation method is applicable for evaluating the liquefaction probability. In addition, since the frequency of the seismic load at the studied site will affect the liquefaction probability, the seismic hazard should be evaluated properly. An illustrative example, assessing the liquefaction probability of an earth-fill dam in Japan, is presented to demonstrate the capability of the proposed method. Finally, the spatial average of the liquefaction probability of the dam over the next 50 years is calculated. The proposed procedure is confirmed to work well for actual design problems.

Identification of electrical and petrophysical rock types based on core and well logs: Utilizing the results to delineate prolific zones in deep water sandy packages from the Shah Deniz gas field in the south Caspian Sea basin

The construction of a geological model using geostatistical simulation methods plays a key role in reservoir characterization. One of the most important parts of the reservoir studies is the discovery of high reservoir quality areas which in turn are dependent on the depositional properties of reservoir formations. The main properties of rocks are usually defined by descriptive studies in the field or on cored intervals (lithofacies) and in more detail by laboratory studies (petrofacies and microfacies). Considering that the core data are not available in all drilled wells, the main focus is on well log data which exist in all wells, continuously. In the current study, with focus on deep sandy reservoir deposits in the southern area of the Caspian Sea, reservoir rocks were characterized based on the integration of the results from core description, thin section studies and petrophysical analysis of well logs in one of the sandstone fields with five members. Balakhany and Fasila sand members are the deepest sandstone members of the field. First petrofacies were identified based on core description and thin section studies. In the next step, the electrical facies and electro-petrophysical facies were investigated using a combination of log data and petrophysical parameters based on Multi-Resolution Graph Base Clustering (MRGC) method. The results from this stage in combination with identified petrofacies were used for the construction of a reservoir static model. Finally, a 3D model of reservoir facies was obtained in each reservoir zone by using variogram analysis and sequential indicator simulation (SIS). The constructed facies model suggests that the high reservoir quality facies are well distributed in the sandy package 3 (SP3) and Fasila B.

Closed-loop field development with multipoint geostatistics and statistical performance assessment

Closed-loop field development (CLFD) optimization is a comprehensive framework for optimal development of subsurface resources. CLFD involves three major steps: 1) optimization of full development plan based on current set of models, 2) drilling new wells and collecting new spatial and temporal (production) data, 3) model calibration based on all data. This process is repeated until the optimal number of wells is drilled. This work introduces a new CLFD implementation for complex systems described by multipoint geostatistics (MPS). Model calibration is accomplished in two steps: conditioning to spatial data by a geostatistical simulation method, and conditioning to production data by optimization-based PCA. A statistical procedure (TruMAP) is presented to assess the performance of CLFD. For performance assessment by TruMAP, the methodology is applied to an oil reservoir example for 25 different true-model cases. Application of a single-step of CLFD, improved the true NPV in 64%–80% of cases. The full CLFD procedure (with three steps) improved the true NPV in 96% of cases, with an average improvement of 37%. These results indicate that probability of improving true NPV increases with closed-loop step. This massive computational experiment involved about 9.5 million reservoir simulation runs that took about 320,000 CPU hours.

Developing a geostatistical simulation method to inform the quantity and placement of new monitors for a follow-up air sampling campaign.

Sampling campaign design is a crucial aspect of air pollution exposure studies. Selection of both monitor numbers and locations is important for maximizing measured information, while minimizing bias and costs. We developed a two-stage geostatistical-based method using pilot NO2 samples from Lanzhou, China with the goal of improving sample design decision-making, including monitor numbers and spatial pattern. In the first step, we evaluate how additional monitors change prediction precision through minimized kriging variance. This was assessed in a Monte Carlo fashion by adding up to 50 new monitors to our existing sites with assigned concentrations based on conditionally simulated NO2 surfaces. After identifying a number of additional sample sites, a second step evaluates their potential placement using a similar Monte Carlo scheme. Evaluations are based on prediction precision and accuracy. Costs are also considered in the analysis. It was determined that adding 28-locations to the existing Lanzhou NO2 sampling campaign captured 73.5% of the total kriged variance improvement and resulted in predictions that were on average within 10.9 μg/m3 of measured values, while using 56% of the potential budget. Additional monitor sites improved kriging variance in a nonlinear fashion. This method development allows for informed sampling design by quantifying prediction improvement (accuracy and precision) against the costs of monitor deployment.

Joint Bayesian inversion based on rock-physics prior modeling for the estimation of spatially correlated reservoir properties

The joint inversion of seismic data for elastic and petrophysical properties is an inverse problem with a nonunique solution. There are several factors that impact the accuracy of the results, such as the statistical rock-physics relations and observation errors. We have developed a general methodology to incorporate a linearized rock-physics model in a multivariate multimodal prior distribution for Bayesian seismic linearized inversion. The prior distribution is used to define a mixture model for elastic and petrophysical properties and introduce physics-based correlations between the properties. Using the rock-physics prior model and a convolutional seismic forward model in the Bayesian inversion framework, we obtain an analytical expression of the spatially independent conditional distributions to be used as a proposal distribution in a Gibbs sampling algorithm. We then combine the sampling algorithm with geostatistical simulation methods to compute the spatially correlated posterior distribution of the model parameters. We apply our method to a real angle-stack seismic data set to generate multiple geostatistical realizations of facies, P-velocity, S-velocity, density, porosity, and water saturation. The method is validated through a blind well test and a comparison with the standard Bayesian linearized inversion.

Synthesis of instrumentally and historically recorded earthquakes and studying their spatial statistical relationship (A case study: Dasht-e-Biaz, Eastern Iran)

Earthquake catalogues are the main source of statistical seismology for the long term studies of earthquake occurrence. Therefore, studying the spatiotemporal problems is important to reduce the related uncertainties in statistical seismology studies. A statistical tool, time normalization method, has been determined to revise time-frequency relationship in one of the most active regions of Asia, Eastern Iran and West of Afghanistan, (a and b were calculated around 8.84 and 1.99 in the exponential scale, not logarithmic scale). Geostatistical simulation method has been further utilized to reduce the uncertainties in the spatial domain. A geostatistical simulation produces a representative, synthetic catalogue with 5361 events to reduce spatial uncertainties. The synthetic database is classified using a Geographical Information System, GIS, based on simulated magnitudes to reveal the underlying seismicity patterns. Although some regions with highly seismicity correspond to known faults, significantly, as far as seismic patterns are concerned, the new method highlights possible locations of interest that have not been previously identified. It also reveals some previously unrecognized lineation and clusters in likely future strain release.

Multiple point statistical simulation using uncertain (soft) conditional data

Geostatistical simulation methods have been used to quantify spatial variability of reservoir models since the 80s. In the last two decades, state of the art simulation methods have changed from being based on covariance-based 2-point statistics to multiple-point statistics (MPS), that allow simulation of more realistic Earth-structures. In addition, increasing amounts of geo-information (geophysical, geological, etc.) from multiple sources are being collected. This pose the problem of integration of these different sources of information, such that decisions related to reservoir models can be taken on an as informed base as possible. In principle, though difficult in practice, this can be achieved using computationally expensive Monte Carlo methods. Here we investigate the use of sequential simulation based MPS simulation methods conditional to uncertain (soft) data, as a computational efficient alternative. First, it is demonstrated that current implementations of sequential simulation based on MPS (e.g. SNESIM, ENESIM and Direct Sampling) do not account properly for uncertain conditional information, due to a combination of using only co-located information, and a random simulation path. Then, we suggest two approaches that better account for the available uncertain information. The first make use of a preferential simulation path, where more informed model parameters are visited preferentially to less informed ones. The second approach involves using non co-located uncertain information. For different types of available data, these approaches are demonstrated to produce simulation results similar to those obtained by the general Monte Carlo based approach. These methods allow MPS simulation to condition properly to uncertain (soft) data, and hence provides a computationally attractive approach for integration of information about a reservoir model.

Spatial modelling of zonality elements based on compositional nature of geochemical data using geostatistical approach: a case study of Baghqloom area, Iran

Due to the existence of a constant sum of constraints, the geochemical data is presented as the compositional data that has a closed number system. A closed number system is a dataset that includes several variables. The summation value of variables is constant, being equal to one. By calculating the correlation coefficient of a closed number system and comparing it with an open number system, one can see an increase in the values of the closed number system, which is false. Such features of this data prevent the application of standard statistical techniques to process the data. Therefore, several methods have been proposed for transforming the data from closed to open number systems. There are various geostatistical methods consisting of estimation and simulation methods in order to model a deposit. Geostatistical simulations can produce various models for a deposit with different probability percentages. The most applicable geostatistical simulation method is the sequential Gaussian simulation technique, which is highly flexible. In this work, 392 Litho-geochemical data of the Baghqloom region of Kerman in Iran consisting of 20 elements were at first converted using an open number system. Afterwards, the elements that were helpful for exploring the area and were normally standard were simulated for 100 times. After the simulations, the valid output was chosen using geostatistical validation. The maps derived from the simulations revealed the enriched concentrations of mineralization elements in the central regions.

Directional hydrostratigraphic units simulation using MCP algorithm

Understanding the geological uncertainty of hydrostratigraphic models is important for risk assessment in hydrogeology. An important feature of sedimentary deposits is the directional ordering of hydrostratigraphic units (HSU). Geostatistical simulation methods propose efficient algorithm for assessing HSU uncertainty. Among different geostatistical methods to simulate categorical data, Bayesian maximum entropy method (BME) and its simplified version Markov-type categorical prediction (MCP) present interesting features. In particular, the zero-forcing property of BME and MCP can provide a valuable constrain on directional properties. We illustrate the ability of MCP to simulate vertically ordered units. A regional hydrostratigraphic system with 11 HSU and different abundances is used. The transitional deterministic model of this system presents lateral variations and vertical ordering. The set of 66 (11 × 12/2) bivariate probability functions is directly calculated on the deterministic model with fast Fourier transform. Despite the trends present in the deterministic model, MCP is unbiased for the HSU proportions in the non-conditional case. In the conditional cases, MCP proved robust to datasets over-representing some HSU. The inter-realizations variability is shown to closely follow the amount and quality of data provided. Our results with different conditioning datasets show that MCP replicates adequately the directional units arrangement. Thus, MCP appears to be a practical method for generating stochastic models in a 3D hydrostratigraphic context.

Phase annealing for the conditional simulation of spatial random fields

Simulated annealing (SA) is a popular geostatistical simulation method as it provides great flexibility. In this paper possible problems of conditioning its realizations are discussed. A statistical test to recognize whether the observations are well embedded in their simulated neighborhood or not is developed. A new simulated annealing method, phase annealing (PA), is presented which makes it possible to avoid poor embedding of observations. PA is based on the Fourier representation of the spatial field. Instead of the individual pixel values, phases corresponding to different Fourier components are modified (i.e. shifted) in order to match prescribed statistics. The method treats neighborhoods together and thus avoids singularities at observation locations. It is faster than SA and can be used for the simulation of high resolution fields. Examples demonstrate the applicability of the method.

Statistical and geostatistical approaches to study spatiotemporal uncertainties in earthquake catalog databases

Lack of data in time and location domain is a main problem in statistical seismology. Therefore, several methods have been developed to make synthetic catalogs to compensate these uncertainties. The approaches introduced in this paper try to decrease these uncertainties. A statistical tool has been determined to study the lack of data in time domain and revise a magnitude-frequency relationship in one of the very high-seismic regions in the world, northwest of Iran (44.00 to 48.00 E and 35.60 to 40.00 N). A geostatistical simulation method has also been utilized to reduce the uncertainties in the spatial domain. Variogram tool has been introduced as a spatial tool to assign a spatial regression among data in various directions. A geostatistical simulation produced a representative, synthetic catalog with 10,707 events to reduce spatial uncertainties (the number of data in the earthquake catalogs was 307 events). The synthetic database was classified by a geographical information system (GIS) based on simulated magnitudes to reveal the underlying seismicity patterns. Although some regions with high seismicity correspond to known faults, as far as seismic patterns are concerned, the new method could highlight locations of interest that have not been previously determined. It also suggests some previously unrecognized lineation and clusters in likely future strain release.