We propose a smooth stochastic process for modeling the vertical well path uncertainty. This process describes the accumulation of measurement errors along the well path. We combine the stochastic process with a stochastic model for surfaces into a consistent framework for simultaneous prediction of well paths and surfaces. We show properties of the proposed stochastic process and provide examples of interaction between wells and surfaces.

When 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.

Kriging is a standard method for conditioning surfaces to observations. Kriging works for vertical wells, but may produce surfaces that cross horizontal wells between surface observations. We establish an approach that also works for horizontal wells, where surfaces are modeled as a set of correlated Gaussian random fields. The constraints imposed by the horizontal wells makes the conditional surfaces non-Gaussian. We present a method for exact conditional simulation and an approximation for prediction and prediction uncertainty. Thousands of constraints can be handled efficiently without numerical instabilities. The approach is illustrated with synthetic and real examples that show how the constraints influence the surfaces and reduce uncertainty.

We investigate how synthetic deposits generated by process-based models can be utilized in classical facies modeling workflows. As an initial approach we focus on using the synthetic deposits as training images for the Multipoint facies simulation algorithm (MPS). We enhance the MPS algorithm to accept multiple training images, and use varied deposits to represent some of the variability one can expect in natural systems.

We generate independent Gaussian random variables on a regular grid and use a spatial filter to smooth the independent random variables to obtain a spatially correlated Gaussian random field. The FFT is used to speed up the smoothing since convolution is a simple cell by-cell multiplication in the Fourier domain. A representation of the spatial convolution filter in the Fourier domain is efficiently obtained from the FFT of any stationary correlation function. Since FFT is cyclic, the grid must be padded to ensure that opposite sides are uncorrelated. The size of the padding is discussed in detail. Most standard covariance functions fail to be positive definite on finite cyclic domains. This causes striping artifacts in the final simulated realizations and failure to meet statistical properties such as variogram reproduction in the simulated realizations. These problems are addressed and solutions are provided to ensure near perfect statistical properties of the generated realizations. The method is fast and can generate a hundred million grid cell realization in approximately 1.5 minutes on a standard laptop PC. The method scales approximately linearly in the number of grid cells.

Unknown values of a random field can be predicted from observed data using kriging. As data sets grow in size, the computation times become large. To facilitate kriging with large data sets, an approximation where the kriging is performed in sub-segments with common data neighborhoods has been developed. It is shown how the accuracy of the approximation can be controlled by increasing the common data neighborhood. For four different variograms, it is shown how large the data neighborhoods must be to get an accuracy below a chosen threshold, and how much faster these calculations are compared to the kriging where all data are used. Provided that variogram ranges are small compared to the domain of interest, kriging with common data neighborhoods provides excellent speed-ups (2–40) while maintaining high numerical accuracy. Results are presented both for data neighborhoods where the neighborhoods are the same for all sub-segments, and data neighborhoods where the neighborhoods are adapted to fit the data densities around the sub-segments. Kriging in sub-segments with common data neighborhoods is well suited for parallelization and the speed-up is almost linear in the number of threads. A comparison is made to the widely used moving neighborhood approach. It is demonstrated that the accuracy of the moving neighborhood approach can be poor and that computational speed can be slow compared to kriging with common data neighborhoods.

Kriging is a data interpolation method that can be used to populate regular grids from data scattered in space, and requires the solution of a linear equation system the size of the number of data. When the data is numerous the speed of the calculation is slow. In this paper we propose to divide the regular grid into rectangular sub-segments and let all the grid cells in each sub-segment share a common data neighborhood. The advantage of this approach is that the number of data in the neighborhoods can be small compared to the complete dataset and it is possible to reuse some of the computations for all grid cells in each sub-segment. We show that the precision can be controlled through selection of neighbourhood size, and that the speed of the calculations can be optimized through selection of sub-segment size. We show that this is an efficient method for kriging when number of data is huge, giving a significant speed-up even for high data densities and precisions.

We demonstrate accurate surface predictions by imposing consistent physical and stochastic relationships between surfaces. The accuracy is improved by using all relevant information collected in wells: well markers, zone logs in horizontal sections, and gas/fluid content along wells. The conditioned surfaces are used to provide estimates of gross rock volumes of oil and gas reservoirs. In particular, we show how spill point and zone log information affect trapped volumes. We apply plain rejection sampling techniques to deal with the highly non-linear relationships between a surface and its spill point. For well path conditioning we build upon an extension of kriging to treat inequality constraints, based on an efficient rejection sampling from a high dimensional truncated multivariate Gaussian distribution. A fast approximate approach to simulating surfaces is presented and successfully applied to estimate volumes. The impact on gross rock volume distributions from different assumptions and data types is demonstrated by several examples and the uncertainties in all the involved data types are consistently handled and quantified.

Fault models are often based on interpretations of seismic data that are constrained by observations of faults and associated strata in wells. Because of uncertainties in depth migration, seismic interpretations and well data, there often is significant uncertainty in the geometry and position of the faults. Fault uncertainty impacts determinations of reservoir volume, flow properties and well planning. Stochastic simulation of the faults is important for quantifying the uncertainties and minimizing the impacts. In this paper, a framework for representing and modeling uncertainty in fault location and geometry is presented. This framework can be used for prediction and stochastic simulation of fault surfaces, visualization of fault location uncertainty, and assessments of the sensitivity of fault location on reservoir performance. The uncertainty in fault location is represented by a fault uncertainty envelope and a marginal probability distribution. To be able to use standard geostatistical methods, quantile mapping is employed to construct a transformation from the fault surface domain to a transformed domain. Well conditioning is undertaken in the transformed domain using kriging or conditional simulations. The final fault surface is obtained by transforming back to the fault surface domain. Fault location uncertainty can be visualized by transforming the surfaces associated with a given quantile back to the fault surface domain.

The gross-rock volume often accounts for the largest uncertainty in reserves. It is therefore important to obtain a correct gross-rock volume distribution and to reduce the uncertainty by using all available data. We demonstrate a way of obtaining accurate volume estimates by imposing realistic and consistent physical and stochastic relationships between the surfaces and hydrocarbon contacts that define the reservoir rock volume. The uncertainty is reduced by using all relevant information collected in wells; well markers, zone logs in horizontal sections, and gas/fluid content along wells. Uncertainties in all these data types are handled. The impact on volume distributions from different assumptions and data types are demonstrated by several examples. We will in particular demonstrate how restrictions on the possible spill point depth have impact on the potential trap size and the trapped volume. Some of the results are obtained using standard stochastic simulation (Monte Carlo) techniques but in particular the highly non-linear relationship between a surface and its spill point requires rejection sampling techniques. Rejection sampling is simple but very inefficient so a fast approximate approach to simulating surfaces is investigated. The conclusion is that the approximation works for calculating volumes but individual surface realizations have unacceptable artefacts.

A multigrid Markov mesh model for geological facies is formulated by defining a hierarchy of nested grids and defining a Markov mesh model for each of these grids. The facies probabilities in the Markov mesh models are formulated as generalized linear models that combine functions of the grid values in a sequential neighborhood. The parameters in the generalized linear model for each grid are estimated from the training image. During simulation, the coarse patterns are first laid out, and by simulating increasingly finer grids we are able to recreate patterns at different scales. The method is applied to several tests cases and results are compared to the training image and the results of a commercially available snesim algorithm. In each test case, simulation results are compared qualitatively by visual inspection, and quantitatively by using volume fractions, and an upscaled permeability tensor. When compared to the training image, the method produces results that only have a few percent deviation from the values of the training image. When compared with the snesim algorithm the results in general have the same quality. The largest computational cost in the multigrid Markov mesh is the estimation of model parameters from the training image. This is of comparable CPU time to that of creating one snesim realization. The simulation of one realization is typically ten times faster than the estimation.

The use of horizontal well data in 3D reservoir modeling has become an increasingly important task as the use of horizontal wells has become common practice. Standard gridding approaches are based on the use of well picks to define the positions of stratigraphic surfaces along well bores. Horizontal wells however, are often drilled almost parallel to the stratigraphic layering so the number of horizons intersected along a horizontal well can be relatively few. Therefore, horizontal sections of the well can be used to constrain the structural position of reservoir zones. A robust, geostatistical approach has been developed to ensure consistent use of horizontal well data in the construction of 3D structural models. Kriging is used for prediction of surface location based on well picks and constraints obtained from zone logs along horizontal wells. In contrast to standard approaches, all well data (picks and constraints) from all surfaces are treated simultaneously and will have impact on all surfaces above and below. The geostatistical approach is fast and reproducible, and allows structural models to be updated continuously as new wells are drilled. The uncertainty can be evaluated by kriging error maps or by generating stochastic realizations that honor all the well data.

Models used for evaluating CO2 plume behaviour in the subsurface often employ simplified geological reservoir descriptions. Experiences from the petroleum industry show, however, that geological heterogeneities significantly influence fluid flow. The present study addresses the need for evaluating the impact of realistic geology on CO2 behaviour in the subsurface. We here demonstrate the effect of adding realistic complexity to the top reservoir morphology. A sensitivity matrix consisting of combinations of depositional and structural irregularities creating relief along the top of a reservoir was set up and the resulting models run in a fluid flow simulator monitoring CO2 plume dynamics. Results demonstrate the interaction between specific geological features and resulting plume behaviour and added retention capacity. Our study highlights the need to include realistic geology in models forecasting migration behaviour in subsurface reservoir.

The primary trapping mechanism in CO2 storage is structural trapping, which means accumulation of a CO2 column under a deformation in the caprock. We present a study on how different top-seal morphologies will influence the CO2 storage capacity and migration patterns. Alternative top-surface morphologies are created stochastically by combining different stratigraphic scenarios with different structural scenarios. Stratigraphic surfaces are generated by Gaussian random fields, while faults are generated by marked point processes. The storage capacity is calculated by a simple and fast spill-point analysis, and by a more extensive method including fluid flow simulation in which parameters such as pressure and injection rate are taken into account. Results from the two approaches are compared. Moreover, by generating multiple realizations, we quantify how uncertainty in the top-surface morphology impacts the primary storage capacity. The study shows that the morphology of the top seal is of great importance both for the primary storage capacity and for migration patterns.

Abstract Structural updates for a complex reservoir model require time-consuming manual work, therefore, updates are rarely performed. This leads to an outdated model that gradually loses its predictability. Eventually, this results in model breakdown, and a new model must be built from scratch. Continuously updatable reservoir models avoid this and increase the value of models as a tool in decision making. In addition, easily updateable structural surfaces enable several structural realizations for spanning the uncertainty. We present the use of a method for fast and robust updates of structural surfaces in reservoir models. We will focus on updates using zone data from horizontal wells (zone-log conditioning), since this traditionally has been a bottleneck that needs tedious manual work prone to error. In zone-log conditioning, we try to generate horizon surfaces that honor the geological zonation along the well paths. This is important for property modeling, and is crucial for fluid-flow simulations. Our method is robust, fully automated, and is built on a consistent mathematical framework that includes specified input-data uncertainties. It has provided satisfactory results for large real-world reservoir models where standard methods and work processes have failed. The field example presented shows a reduction from 22.9 % to 0.9 % in incorrectly honoring of the zone logs by applying this method rather than the standard approach. The remaining 0.9 % is due to conflicting data, gridding errors, and is difficult to get rid of even with manual editing. We consider this a large step forward with respect to providing an up-to-date basis for decisions that also can account for structural uncertainties.

The robust use of horizontal well data in 3D stratigraphic horizon modeling is an important challenge.In many reservoir models there are inconsistencies between the zone log information observed inhorizontal wells, and the horizons and zones of the 3D model.The two main reasons for these inconsistencies are that:1) Zone log information is not used directly in the modeling of the surfaces;2) Well picks in deviated wells do not impact the modeling of adjacent surfaces.Inconsistencies between the 3D model and the horizontal well data are typically resolved by timeconsuming, iterative manual editing of the stratigraphic horizons.A robust, geostatistical approach for ensuring the correct modelling of multiple stacked stratigraphicsurfaces constrained by long horizontal wells is presented. The approach has been developed tomaintain 3D consistency between the stratigraphic surfaces and horizontal well data. The stratigraphicsurfaces are treated simultaneously in 3D and universal or Bayesian cokriging is used for prediction ofsurface location based on a variety of constraints including: well picks; zone logs; isochores andseismic interval velocities.In contrast to standard approaches, all well markers (picks) are treated simultaneously and will haveimpact on surfaces above and below. This ensures consistent use of all available well marker data forall surfaces.Zone log constraints in long horizontal wells are handled by identifying the sections of the well pathsthat impose soft (inequality) constraints on the surfaces. The key sections of the horizontal wells are re-sampled at approximate grid resolution and used as help points in the kriging equations. Identifyingthe soft constrains and assigning correct value to the help points is the cornerstone of the method.Velocity and isochore trends are used to ensure that overall shape of the structure and thicknessvariations are preserved. Uncertainty trends for velocities and isochores are used to allow flexibility inthe relative weights assigned to the various input constraints.

Abstract Fault geometry is modelled on basis of seismic data, but restricted by fault observations in wells. Due to uncertainties in depth migration, seismic interpretation and well data, there is a significant uncertainty in the geometry and position of the faults. Fault uncertainty impact reservoir volume, flow properties and well planning, and can be studied by stochastic simulation of faults. We have developed a method for stochastic simulation of fault surfaces and fault networks using standard geostatistical methods. This is made possible by the fault parameterization used, where the faults are modelled as tilted surfaces. This new method is more flexible and efficient compared to already existing algorithms due to a simpler parameterization. Conditioning to fault observations in wells is also made simpler. The fault is defined as a two-dimensional surface on a tilted reference plane. The uncertainty for a fault surface is bounded by a volume enclosing the fault surface. The smoothness of the simulated fault surfaces is controlled by variograms. The simulation is done by adding a simulated Gaussian residual. Well conditioning is done by kriging. Using the described method we can simulate a set of fault realizations where the simulated faults look realistic, are within the defined uncertainty volumes, and honour well observations. Technical contributions compared to previous work include efficient simulation of fault geometry, a flexible uncertainty model and well conditioning with no performance impact.

Geostatistical simulation of permeabilities on a geologically detailed resolution will account for permeabilities in the cells containing wells through kriging. These permeabilities are typically based on porosity logs and core plug measurements of both porosity and permeability. No direct measurement of permeability on the geomodelling grid scale exists. However, well tests give information on the effective permeability in an area close to the well region, covering several grid cells, and are therefore data on an aggregate scale. Given an assumption of radial flow into a vertical well, the effective permeability becomes a convolution of the permeabilities in the well test region. Downscaling of the well test data is possible by co-kriging the aggregate scaled permeability field with that of the geomodelling grid scale. Thereby, the geomodelling grid permeabilities will honour data on both scales. The effective permeability is downscaled through inverse block kriging, which implies a deconvolution procedure. Keeping the computation costs low when introducing a new conditioning parameter is ensured by transforming into the Fourier domain, since the Fast Fourier Transform algorithm is an efficient method for solving the inverse block co-kriging. Initial testing of the implemented algorithm has been made on real data from a StatoilHydro operated field on the Norwegian continental shelf in a proof of concept test. Three wells with well test data were chosen, and their derived effective permeabilities were included in the permeability simulations. Cases both with and without well test conditioning were run, and compared in a well test simulation software. All three near well regions showed a significant improvement. These tests indicate that the conditioning method can be a useful contribution in bringing dynamic data into the reservoir characterization.