Per Røe

We demonstrate the use of Ensemble Kalman Inversion (EKI) for building three-dimensional, multi-fault, kinematically restorable structural geologic models, by means of a workflow in which fault geometry, the distribution of slip on a fault, and the geometry of folded horizons are all modeled. The models are constrained by observations of faults and horizons in the present deformed state together with the expectation that horizons should restore flat. Two modeling approaches are tested: restoration from the deformed state, and forward modeling from the restored state. We first test these methods on a synthetic model involving a single fault, and then apply them to a real-field example involving five faults. Models are prone to ensemble collapse, which results in underestimation of uncertainty at small ensemble sizes, but localization and covariance inflation can mitigate this issue. With these methods, EKI can recover the true parameter values in the synthetic case and produce a solution consistent with the data in the real case, as well as quantify uncertainty in both cases. EKI, therefore, shows promise as a tool for building complex, restorable structural geologic models, and it holds the potential for integration of fault kinematics with other ensemble-based subsurface modelling workflows.

The GIG, Geophysical Inversion to Geology, consortium is a research consortium run by Norwegian Computing Center, with the aim of developing new understanding, new methods and software for obtaining reservoir properties from geophysical measurements This note gives a summary of the work in 2021 and the suggested work plan for 2022. The final work plan for 2022 will be decided on the annual meeting in Jan 2022.

Structural restoration using kinematic principles of fault-related folding is a valuable tool in building realistic geological models. Models are, however, typically uncertain and non-unique. While data inversion methods can be employed to find a best-fit model and estimate uncertainty, their use is limited to relatively simple models involving a single fault in two dimensions. In this work, we employ an iterative form of the Ensemble Kalman Filter (EnKF) together with a kinematic model for deformation around normal faults to build and restore three-dimensional structural geologic models. The EnKF is ideally suited to data inversion problems that involve large numbers of model parameters and is frequently used in reservoir simulations, which often do not consider uncertainty in geologic structure. We develop a workflow in which fault geometry, the distribution of slip on a fault, and the geometry of folded horizons are all modelled using the EnKF. The models are constrained by observations of faults and horizons in the present deformed state together with the expectation that horizons should restore flat. We test two modelling approaches: a restoration-based approach in which the model is built in the deformed state and then restored, and a forward modelling-based approach in which the model is built in the restored state and then forward modelled to match present-day data. We test these methods first on a synthetic model involving a single fault and then on a real-world example involving five faults. Both the restoration- and forward modelling-based methods work well for the synthetic model, but forward modelling works best for the more complex real-world case study. The EnKF method provides not only a best-fit model but also an estimate of model uncertainty. The estimation of uncertainty is, however, hindered by the phenomenon of ensemble collapse, which results in underestimation of the uncertainty in model parameters at small ensemble sizes. We employ bootstrap-based localization and covariance inflation to help mitigate this issue and use a dummy parameter to identify whether significant ensemble collapse has occurred. While ensemble collapse remains a challenge in some cases, the EnKF nonetheless shows considerable promise as a tool for building complex many-parameter structural models that are kinematically restorable, and it holds the potential for future integration of structural modelling with other EnKF-based workflows in subsurface modelling.

Subseismic faults are small faults or fractures that may be difficult to determine but can have large consequences for fluid flow and pressure communication in the subsurface. Thus, knowing their distributions may be important in several subsurface applications, such as hydrocarbon exploration and exploitation, geothermal energy production, and subsurface CO2 injection. The aim of this work is to use a stochastic model to populate a three-dimensional structural model of the subsurface with subseismic faults. The novelty of the proposed method is the conditioning of the stochastic model to input maps describing displacement and stress orientation along subsurface horizons. Hence, the resulting structural model will be consistent with these maps. The maps can originate from a variety of sources, for example, predictions of a geomechanical model or (indirect) measurements of subsurface displacements and stresses. The model uses simulated annealing as the optimization algorithm, where the residual between the displacement of the modeled subseismic faults and the input displacement map is minimized through an iterative process. Each subseismic fault is modeled with a three-dimensional displacement field around the fault slip plane, enabling comparisons with the input displacement map along a horizon. An example of how the model distributes the subseismic faults around larger known faults, using a synthetically created displacement map, is provided. The result shows that the model quickly converges towards a set of subseismic faults, giving total displacement and strike orientation close to the input maps.

The GIG, Geophysical Inversion to Geology, consortium is a research consortium run by Norwegian Computing Center, with the aim of developing new understanding, new methods and software for obtaining reservoir properties from geophysical measurements This note gives a status of the GIG consortium after its fourth year of operation, and also gives the suggested work plan for 2021.

We have developed an efficient methodology for Bayesian prediction of lithology and pore fluid, and layer-bounding horizons, in which we include and use spatial geologic prior knowledge such as vertical ordering of stratigraphic layers, possible lithologies and fluids within each stratigraphic layer, and layer thicknesses. The solution includes probabilities for lithologies and fluids and horizons and their associated uncertainties. The computational cost related to the inversion of large-scale, spatially coupled models is a severe challenge. Our approach is to evaluate all possible lithology and fluid configurations within a local neighborhood around each sample point and combine these into a consistent result for the complete trace. We use a one-step nonstationary Markov prior model for lithology and fluid probabilities. This enables prediction of horizon times, which we couple laterally to decrease the uncertainty. We have tested the algorithm on a synthetic case, in which we compare the inverted lithology and fluid probabilities to results from other algorithms. We have also run the algorithm on a real case, in which we find that we can make high-resolution predictions of horizons, even for horizons within tuning distance from each other. The methodology gives accurate predictions and has a performance making it suitable for full-field inversions.

Sedimentary rocks often obey vertical transvers isotropy due to the nature of the sedimentation process or the inherent orientation of rock grain minerals (e.g. clay platelets). Despite several documented cases in the literature, the effect of anisotropy is often neglected in seismic AVO inversion. We study the effect of transverse isotropy in shale encasing isotropic sand using probabilistic AVO inversion. The inversion algorithm is inverting seismic pre-stack data to lithology and fluid probabilities. A synthetic case demonstrates how anisotropy in shale may lead to wrong interpretation of fluid content in underlying sand. A field case in the northern North Sea shows that accounting for anisotropy in shale may improve the discrimination between good and poor sand of an Oxfordian turbidite reservoir. The results are encouraging and consistent with observations of sandy intervals in the wells and the depositional system of the study area. Our results clearly demonstrate that transverse anisotropy in sediments may give a significant contribution to the AVO gradient that otherwise could be misinterpreted. In the study area it becomes especially important, since there is a lack of contrast in acoustic impedance between the encasing shale and the target sand.

HAVANA is a program for simulating faults in petroleum reservoirs, and for integrating the effects of these faults into the reservoir description. The HAVANA project has a long history, the original sponsors being Statoil, BP, and Norsk Hydro. Other sponsors include Conoco Norge AS, Saga Petroleum AS and Centre for Integrated Petroleum Research at University of Bergen.

Lithology and fluid prediction from seismic data is traditionally done in two steps; first an inversion of the seismic data to elastic parameters, and subsequently a prediction of lithology and fluid based on a rock physics model linking the elastic parameters to individual lithology and fluid combinations. Recently, a number of inversion algorithms have been developed that, based on Bayesian statistical methodology, estimate the probability of lithology and fluid directly from seismic data. In this paper we compare the performance of two state-of-the-art Bayesian inversion algorithms on a real data set from the Volund field in the North Sea. The first algorithm follows the traditional two-step approach and cannot take into account the stratigraphic ordering of lithology and fluid. The second algorithm, referred to as one-step, evaluates possible lithology and fluid combinations within a vertical window around each inversion point enabling correct stratigraphic ordering. We find that the one-step inversion resolves more details and honours the data more strongly than the two-step approach. The latter is more prone to return the prior model if information in the seismic data is not sufficiently strong. Both models detect hydrocarbon filled sand injectites that are typical for the field.

Ensemble collapse is one of the main problems encountered when doing history matching using ensemble based methods. In history matching results where ensemble collapse is a problem it is very hard to distinguish updates due to spurious correlation with noise in the data with actual updates attributed to information in the data. Various localization workflows have been proposed for the synthetic case REEK and the public domain dataset from NORNE. We demonstrate moderate improvements of the localized updates for the REEK case. For the NORNE case, issues around the prior distribution were encountered that prevented further investigation of the proposed localization workflow.

Conduit fault zones and fault zones that can accommodate long-distance along-fault flow are well-documented phenomena. In reservoir simulation models, flow within these features is more correctly captured using volumetric representations of fault zones instead of employing standard two-dimensional fault planes. The present study demonstrates a method for generating fault envelope grids on full-field reservoir models, within which fault cores (i.e., regions where most of fault zone displacement is accommodated) are modeled. The modeled fault core elements are lenses and slip zones. They are defined as facies units and populated in the fault envelope grids using combined object-based simulation and deterministic techniques. Using the facies property, four reservoir simulation models are generated by modulating fault core thickness and slip zone type and permeability. Membrane slip zones (slip zones that act as partial barriers to fluid flow) cause the fault cores to form baffle–conduit systems. Along-strike positioned injector–producer pairs focus flow into the fault cores, decreasing sweep efficiency. In contrast, injected fluids of injector–producer pairs positioned to drain perpendicular to the fault cores are partitioned and distributed by the fault cores and therefore increase overall sweep efficiency. In reservoir models with conduit slip zones (slip zones that enhance flow along them and act as partial barriers to flow across them), the fault cores act as thief zones. Fluids preferentially move through the fault cores toward the nearby producers instead of through sedimentary layers with high permeability. Sweep efficiency in the reservoir models with conduit fault cores has less dependency on injector–producer configuration. Our study suggests that the improved realism added by incorporating volumetrically expressed fault cores substantially influences forecasts of field behavior and consequently should be considered during oil and gas production planning.

Underestimation of posterior parameter uncertainty is one of the main problems encountered when doing history matching using ensemble based methods. In history matching results with the partial or full ensembles collapse, it is very hard to distinguish updates due to spurious correlation with noise in the data from the actual updates attributed to information in the data. History matching of porosity and permeability based on well production data using the ensemble smoother with multiple data assimilation has been performed on a synthetic data set. The presence of ensemble collapse has been evaluated by different means: by looking at the stability of the update based on the starting ensemble, by adding dummy parameters to the update which do not affect the forward model, and by examining how well the data set used to generate the production data matches the posterior distributions of the parameters. Ensemble collapse can be avoided by increasing the number of ensembles. This is however a prohibitively expensive strategy for cases with a large number of history data. Localization methods have been proposed in the literature as a way to increase the ensemble spread and hence avoid collapse, by for example limiting the analysis update to regions of influence of the data, while at the same time keeping the number of ensembles low. A local analysis was performed to reduce the problems related to ensemble collapse. The results from the localized history matching produce a posterior distribution that better matches the original data set. Since our test data set is synthetic, we may perform measures of posterior uncertainty estimation by comparing with the true solution, with and without localization.

The spatial structure of the subsurface is an important factor when interpreting seismic data. The Bayesian methodology is a valuable tool for integrating these spatial relations in the inversion process as it merges the information together and assesses the uncertainty of the model. In the everyday use of the Bayesian methodology, however, the computational cost is a challenge. We describe a new approach that utilizes a local neighborhood to include the spatial constraints and assess the uncertainties in the inversion using fast and parallelizable computations. The approach is applicable for both discrete lithology-fluid prediction and estimation of rock properties, such as porosity and saturation.

Underestimation of posterior parameter uncertainty is one of the main problems encountered when doing history matching using ensemble based methods. In history matching results with the partial or full ensembles collapse, it is very hard to distinguish updates due to spurious correlation with noise in the data from the actual updates attributed to information in the data. History matching of porosity and permeability based on well production data using the ensemble smoother with multiple data assimilation has been performed on a synthetic data set. The presence of ensemble collapse has been evaluated by different means: by looking at the stability of the update based on the starting ensemble, by adding dummy parameters to the update which do not affect the forward model, and by examining how well the data set used to generate the production data matches the posterior distributions of the parameters. Ensemble collapse can be avoided by increasing the number of ensembles. This is however a prohibitively expensive strategy for cases with a large number of history data. Localization methods have been proposed in the literature as a way to increase the ensemble spread and hence avoid collapse, by for example limiting the analysis update to regions of influence of the data, while at the same time keeping the number of ensembles low. A local analysis was performed to reduce the problems related to ensemble collapse. The results from the localized history matching produce a posterior distribution that better matches the original data set. Since our test data set is synthetic, we may perform measures of posterior uncertainty estimation by comparing with the true solution, with and without localization.

Corner point grids is currently the standard grid representation for use in reservoir simulation. The cell faces in corner point grids are traditionally represented as bilinear surfaces where the edges between the corner points all are straight lines. This representation has the disadvantage that along faults with varying dip the cell faces on either side will not precisely match, giving overlapping cells or gaps between cells. We propose an alternative representation for the cell faces. The four vertical cell faces are still represented as bilinear surfaces, but instead of having linear edges between the cell corners along the top and bottom faces, we propose a representation of the vertical cell faces where any horizontal intersection will give a straight line, giving column faces whose shape is independent of the corner point locations of the individual grid cells. This ensures that the grid columns match up and that there are no gaps or overlapping volumes between grid cells. This representation gives a local parameterization for the whole grid column, and the top and bottom grid cell surfaces are modeled as bilinear using this parameterization. A set of local coordinates for the grid cell permits all the common grid operations like volume calculation, area calculation for cell faces, and blocking of well traces.

The internal structure and petrophysical property distribution of fault zones are commonly exceedingly complex compared to the surrounding host rock from which they are derived. This in turn produces highly complex fluid flow patterns which affect petroleum migration and trapping as well as reservoir behavior during production and injection. Detailed rendering and forecasting of fluid flow inside fault zones require high-resolution, explicit models of fault zone structure and properties. A fundamental requirement for achieving this is the ability to create volumetric grids in which modeling of fault zone structures and properties can be performed. Answering this need, a method for generating volumetric fault zone grids which can be seamlessly integrated into existing standard reservoir modeling tools is presented. The algorithm has been tested on a wide range of fault configurations of varying complexity, providing flexible modeling grids which in turn can be populated with fault zone structures and properties.

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, flowproperties 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 andstochastic 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.

Corner point grids are currently the standard representation for reservoir simulation grids. Cell faces in corner point grids are traditionally represented as bilinear surfaces where the edges between the corner points all are straight lines. This representation has the disadvantage that at pillar-gridded faults with varying dip the cell faces on either side will not precisely match, giving overlapping cells or gaps between cells. We propose an alternative representation for the cell faces. The four vertical cell faces are still represented as bilinear surfaces, but instead of having linear edges between the cell corners along the top and bottom faces, we propose a representation where all horizontal intersections through the grid give straight lines between the grid pillars, giving column faces whose shape is independent of the corner point locations of the individual grid cells. This ensures that the grid columns match up and that there are no gaps or overlapping volumes between grid cells. This representation gives a local parameterization for the whole grid column, and the top and bottom grid cell surfaces are modelled as bilinear within this domain. Using this representation we get a parameterization of the grid cell which we use to calculate the cell volume.

Conduit faults and faults that can accommodate vast long-distance along-strike flow are well-documented phenomena. In reservoir simulation models, flow within these features are more correctly captured using a volumetric representation of fault zones rather than employing 2D fault planes. We here demonstrate a method for implementing fault zone grids and features on a full-field case study. The fault zone grid is populated by fault rocks and fractures. We investigate the resulting effect on the modelled forecast of field-wide reservoir flow. Membrane slip zones cause the fault zones to form barrier-conduit systems. Along-strike positioned injector-producer pairs focus flow into the fault zone, decreasing sweep efficiency. On the other hand, injector-producer pairs positioned to drain perpendicular to faults partition the injection fluids and therefore tend to increase overall sweep efficiency. In models with conduit slip zones, the fault zones act as thief zones. Fluids preferentially move through the fault zones towards the producers. Consequently, sweep efficiency is more related to injector-producer distance than the geometric relation of well pairs to the faults. Our study suggests that the improved realism added by incorporating volumetrically expressed fault zones substantially influences forecasts of field behavior, and consequently should be considered during oil/gas production planning.

Traditionally fault seal calculations take place directly within the simulation grid. This approach works well for grids where all the faults are aligned along the grid pillars, but implementing an algorithm that works with stair-stepped representation of the faults has proven to be very difficult. Especially the calculation of the displacement field used both indirectly in the fault seal parameter calculation and directly in the calculation of fault zone permeability is challenging. It is hard to find where the different grid layers intersect the fault trace, and the layers are not always completely represented on both sides of the fault. We present a novel algorithm where the calculation of the fault zone permeability is carried out on a 2D plane representing the fault surface. The input parameters needed for calculating the fault zone permeability are resampled from the simulation grid onto the 2D plane, while the resulting fault zone permeability is resampled back into the simulation grid, prior to calculation of the fault transmissibility. The new approach is shown to generate good results both for pillar-faulted grids, and for grids with stairstepped faults, and also works well near complex truncations.

In history matching and sensitivity analysis, flexibility in the structural modelling is of great importance. The ability to easily generate multiple realizations of the model will have impact both on the updating workflow in history matching and uncertainty studies based on Monte Carlo simulations. The main contribution to fault modelling by the work presented in this paper is a new algorithm for calculating a 3D displacement field applicable to a wide range of faults due to a flexible representation. This gives the possibility to apply this field to change the displacement and thereby moving horizons and fault lines. The fault is modelled by a parametric format where the fault has a reference plane defined by a centre point, dip and strike angles. The fault itself is represented as a surface defined by a function z = f(x,y), where x, y and z are coordinates local to the reference plane, with the z-axis being normal to the plane. The displacement associated with the fault outside the fault surface is described by a 3D vector field. The displacement on the fault surface can be found by identifying the intersection lines between horizons and the fault surface (fault lines), and using kriging techniques to fill in values at all points on the surface. Away from the fault surface the displacement field is defined by a monotonic decreasing function which ensures zero displacement at a specified distance from the fault. An algorithm is developed where the displacement can be increased or decreased according to user-defined parameters. This means that the whole displacement field is changed and points on horizons around the fault can be moved accordingly by applying the modified displacement field on them. The interaction between several faults influencing the same points is taken care of by truncation rules and the ordering of the faults. The method is demonstrated on a realistic synthetic case based on a real reservoir.

Traditionally fault seal calculations take place directly within the simulation grid. This approach works well for grids where all the faults are aligned along the grid pillars, but implementing an algorithm that works with stair-stepped representation of the faults has proven to be very difficult. Especially the calculation of the displacement field used both indirectly in the fault seal parameter calculation and directly in the calculation of fault zone permeability is challenging. It is hard to find where the different grid layers intersect the fault trace, and the layers are not always completely represented on both sides of the fault. We present a novel algorithm where the calculation of the fault zone permeability is carried out on a 2D plane representing the fault surface. The input parameters needed for calculating the fault zone permeability are resampled from the simulation grid onto the 2D plane, while the resulting fault zone permeability is resampled back into the simulation grid, prior to calculation of the fault transmissibility. The new approach is shown to generate good results both for pillar-faulted grids, and for grids with stair-stepped faults, and also works well near complex truncations.

Traditionally fault seal calculations take place directly within the simulation grid. This approach works well for grids where all the faults are aligned along the grid pillars, but implementing an algorithm that works with stair-stepped representation of the faults has proven to be very difficult. Especially the calculation of the displacement field used both indirectly in the fault seal parameter calculation and directly in the calculation of fault zone permeability is challenging. It is hard to find where the different grid layers intersect the fault trace, and the layers are not always completely represented on both sides of the fault. We present a novel algorithm where the calculation of the fault zone permeability is carried out on a 2D plane representing the fault surface. The input parameters needed for calculating the fault zone permeability are resampled from the simulation grid onto the 2D plane, while the resulting fault zone permeability is resampled back into the simulation grid, prior to calculation of the fault transmissibility. The new approach is shown to generate good results both for pillar-faulted grids, and for grids with stair-stepped faults, and also works well near complex truncations.

Faults are volumetric in nature and can cause complex fluid flow inside the fault zone because of its special fault zone architecture and different petrophysical properties from the host rock. Thus explicit fault zone modeling is important for capturing the fluid flow inside and through the fault zone precisely. Generation of a refined volumetric fault zone grid is the first step to perform explicit fault zone modeling. An algorithm for generating volumetric fault zones has already been implemented in Havana, however this algorithm failed to generate continuous top and bottom surfaces for the fault zone. This lead to internal discontinuities in the fault zone grid, and made it hard to run flow simulations on the grid. We now present an improved version of the algorithm that works well on complex faults and indicate the capability of explicit fault facies modelling of real field cases.

The existence, position, shape, size and displacement of faults play an important role in several aspects of reservoir modelling including flow simulations, volume calculations and well planning. To model the faults, seismic inputs with both interpretation and migration uncertainty and well observations with uncertainty in log interpretations are applied. Uncertainty in the interpretation of the horizons causes uncertainty in the fault displacement. All these different sources of uncertainty should be included in the fault modelling. The relationship between faults and their size and shape defines the compartments of the reservoir and influences the simulated flow paths. The modelling of the fault displacement can have great impact on whether the fluids flow through the fault or the fault acts as a barrier in the model. The joint uncertainty in the position and size of boundary faults plays a role in the volume estimation, and this uncertainty should be accounted for also in well planning to avoid drilling outside the desired area or directly through a fault. We will present a fault modelling tool where different realizations of the structural model related to the above mentioned uncertainties can be generated. This is done by using a flexible representation of the faults which lets us update the fault geometry, fault displacement and fault size both in deterministic and stochastic ways. We will show how uncertainty envelopes around faults can be generated. These envelopes are based on seismic uncertainty and define the space for the fault to reside in. Both stochastic and deterministic change of position and shape of the fault can be performed within the envelope. The envelope and fault position can be conditioned to well data. The displacement and the length of the faults can be modified by deterministic or stochastic processes. The displacement is modelled by elliptic trends combined with input from interpreted seismic horizons. The deterministic modification implies either a scaling of the throw, leaving the length of the fault unchanged or increasing/decreasing the throw by a constant which means that a new fault tip line is estimated. In both the case of position change and displacement modification, the deterministic workflows are typically used to establish base and min/max scenarios while stochastic realizations are applied in a Monte-Carlo setting to calculate statistical properties. Faults below seismic resolution can be modelled by a stochastic process that distributes smaller faults in the reservoir either through global trends or as secondary faults linked to seismic visible primary faults. The modelling tool allows for setting up different workflows to analyse the sources of uncertainty one at the time or combined to give a total effect. Examples from a synthetic case based on a real reservoir will be given to show different aspects of the uncertainty modelling.

Kurs i bruk av Havana for RMS utviklere.

HAVANA is a program for simulating subseismic faults in petroleum reservoirs, and for integrating the effects of these faults into the reservoir description. The HAVANA project has a long history, the original sponsors being Statoil, BP, and Norsk Hydro. Other sponsors include Conoco Norge AS, Saga Petroleum AS and Centre for Integrated Petroleum Research at University of Bergen.

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.