## Abstract

Despite their significance, structural parameters are sometimes neglected in assessments of uncertainty on connected volumes and forecast production for compartmentalized reservoirs. A workflow is proposed for modelling multiple realizations of fault geometry and properties using 3D geomodelling software. Geometrical parameters that may be simulated include fault and horizon shape and location, fault displacement and fault pattern, while property variables include fault permeability, thickness and clay smears. Realizations are ranked by estimated connected volume, with selected models being exported for numerical flow simulation. Experimental design is used to assess sensitivity of forecast production and pressure to different parameters. The workflow is illustrated using a North Sea reservoir, in which structural heterogeneities cause considerable uncertainty on connected volumes, with implications for history matching and infill well planning. Fault geometry and permeability were the most important properties for all studied responses, however their relative significance could vary between early and late field life. A number of improvements are proposed, chiefly in the areas of connected volume estimation, handling of uncertain grid geometries and calculation of stress- or saturation-dependent fault permeabilities. Finally, the method can be integrated with conventional sedimentary and petrophysical uncertainties to investigate interactions and relative sensitivities with regard to structural parameters.

Faults affect fluid flow in sedimentary basins on multiple length and time-scales. Over a time period of thousands to millions of years, faults can act as vertical hydrocarbon migration pathways (Weber *et al.* 1978). Faults may also impede lateral fluid movement, ultimately providing a trapping mechanism for significant hydrocarbon accumulations (Smith 1966, 1980; Watts 1987; Harding & Tuminas 1989; Knott 1993; Gibson 1994; Berg & Avery 1995) or acting as baffles to economic oil and gas production (Bentley & Barry 1991; Jev *et al.* 1993; Harris *et al.* 2002). Faults affect lateral fluid flow by juxtaposing different sedimentary layers (Allan 1989; Knipe 1997) or due to the petrophysical properties of the fault zone materials themselves, which result from processes such as particulate flow and grain crushing (cataclasis); entrainment of fine-grained material into the fault zone and preferential diagenesis and mineral reprecipitation (Yielding *et al.* 1997).

Sampling of subsurface information is limited to spatially extensive but low resolution seismic data and sparse, high resolution well penetrations. Large (greater than *c*. 10 m displacement) faults are usually avoided as a drilling hazard (Cerveny *et al.* 2004) and very rarely cored (Hippler 1997; Aarland & Skjerven 1998; Knai & Knipe 1998), while typical seismic surveys are designed to image sub-horizontal reflectors rather than faults (van der Poel & Cassell 1989; Worthington & Hudson 2000). This lack of knowledge must therefore be filled by making inferences from observed responses (inverse modelling) or by employing some kind of conceptual forward model to predict fault rock properties from available information.

These tasks have been greatly facilitated by methods developed in the past decade, which allow fault sealing behaviour in siliciclastic sequences to be more realistically quantified and represented in geocellular reservoir models (Lindsay *et al.* 1993; Yielding *et al.* 1997; Manzocchi *et al.* 1999). Algorithms such as Shale Gouge Ratio (SGR), Shale Smear Factor (SSF) and Clay Smear Potential (CSP) provide a quantitative prediction of the amount of fine-grained (phyllosilicate) material contained within a fault zone. Published empirical transformations (Manzocchi *et al.* 1999; Sperrevik *et al.* 2002) provide estimates of fault rock permeability, which may be readily included in numerical reservoir simulations through the use of single-phase transmissibility multipliers on cell connections. Yielding (2002) reported the successful application of this approach to modelling intra-reservoir faults in the Scott field, and similar methodologies have resulted in good history matches for the Snorre (Sverdrup *et al.* 2003) and Heidrun fields (Knai & Knipe 1998), in the latter case using fault properties derived from small cored faults. Fault properties derived from core analysis have also been successfully employed for fields in the UKCS Brent Province (Jolley *et al.* 2007); however the authors attached greater importance to rigorous structural interpretation and quality control of reservoir simulation grids to ensure correct reservoir ‘plumbing’ across and around faults.

The grid type most commonly used in both geological modelling and reservoir simulation is the stratigraphic grid, where the top and base reservoir surfaces are connected by linear pillars aligned parallel to faults. Internal reservoir layering is defined by sub-parallel intermediate surfaces; together with vertical sections between pillars, these define the boundaries of each grid cell. The popularity of these grids stems from the ability to model geological properties and perform flow simulation using a shared model. Stratigraphic grids are, however, a compromise between the needs of these different tasks, which may necessitate simplification of the structural interpretation, for example, modification or removal of some fault surfaces in order to build a coherent grid.

Within the stratigraphic grid, geostatistical algorithms are used to interpolate and/or simulate properties interpreted from well data. As wells only sample a small fraction of the reservoir volume, it is impossible to produce a ‘true’ representation of the subsurface. Instead, geostatistical simulation algorithms aim to produce multiple equally probable realizations that are globally accurate; reproducing the interpreted property continuity while honouring all available data. Advances in both Earth modelling software and computer power have led to the widespread application of such stochastic techniques. In field appraisal and development studies, it is now standard practice to generate hundreds or thousands of equally probable combinations of uncertain input parameters such as sedimentary facies, porosity and permeability. The resulting distributions of connected volumes provide estimates of uncertainty which can aid in decision making and economic forecasting. Save for a few exceptions in the literature (England & Townsend 1998; Ottesen *et al.* 2005; Rivenæs *et al.* 2005), uncertainties in fault geometry and properties are typically excluded from such modelling workflows, despite their potential significance for production forecasts (e.g. Lia *et al.* 1997; Damsleth *et al.* 1998). Overlooking these factors could lead to sub-optimal or even failed field developments and associated unnecessary costs. To avoid such pitfalls, uncertainties in fault geometry and properties should be modelled and their influence on connected volumes assessed as for sedimentary and petrophysical properties.

The aim of this paper is to outline a workflow for modelling uncertainty on connected hydrocarbon volumes and forecast production in reservoirs that are compartmentalized by faults. Two categories of uncertainty are defined: geometrical and property uncertainties. The shape and location of fault and horizon surfaces are the principal geometrical factors, which arise from uncertainties in the acquisition, processing and interpretation of reflection seismic data. Property uncertainties include: fault permeability; thickness; and clay distribution. These parameters are estimated using appropriate analogue databases and empirical functions relating them to available model properties such as fault throw and SGR. Uncertainties can be defined on the constants and coefficients used in these empirical functions and on the underlying spatial continuity of the properties. The defined geometrical and property modelling parameters are then combined using stochastic simulation to generate hundreds of equally probable realizations of reservoir geometry and properties. Realizations may be ranked according to connected hydrocarbon volumes, and selected models exported for numerical fluid flow simulation.

The interpretation and modelling of geometrical and property uncertainties will be first described. The workflow will then be illustrated using a field case from the UK sector of the North Sea. Results will be presented as distributions of connected volumes and forecast cumulative production; sensitivity of these outputs to different factors will be analysed using experimental design techniques. Finally, we will discuss the impact of these results on this field case and propose some extensions and improvements to the workflow.

## Geometrical uncertainties

The first step in the method is the interpretation of seismic data to build a coherent three dimensional network of faults and horizons. Geometrical uncertainties arising from seismic data are estimated and then modelled using a proprietary software package for stochastic simulation of reservoir structure. The processes for uncertainty estimation and stochastic simulation are described respectively by Thore *et al.* (2002) and Lecour *et al.* (2001) and summarized below.

Uncertainties can be considered on six phases of the seismic ‘chain’, namely acquisition, pre-processing, stacking, migration, interpretation and time-to-depth conversion. Each source of uncertainty is described as a correlated vector field consisting of uncertainty direction, magnitude and correlation length at each point of a given surface. All surfaces are considered to be probabilistic, in each case located somewhere within an uncertainty volume surrounding the ‘best estimate’ interpretation. Faults are parameterized as objects with a central ‘backbone’ curve connecting secondary horizontal ‘generator’ curves, along which uncertainty values are stored at each node. Uncertainty vectors are specified at a few locations along the fault and interpolated where no information is provided. Horizons are described by triangulated surfaces, again with uncertainty values stored at each node. Finally, a set of constraints describes the geometric connections between faults and horizons, and an additional scalar property describes the correlation coefficient between horizons, in order to prevent undesired horizon crossing.

Generation of surface realizations is performed using probability field (P-field) simulation, which consists of adding a spatially correlated random variable to the ‘best estimate’ interpretation. Using a constant random number for the random field will produce a bulk lateral shift of the fault; while using different random numbers at the top and base of the backbone will cause the fault's dip to vary. Fault geometries are preserved by using a long correlation length along the horizontal axis and a monotone function (i.e. no inversion of fault curvature) in the *Z*-axis. After simulation, a post-processing geometrical optimization is used to ensure each realization obeys the prior fault-horizon constraints. Stratigraphic grids are deformed by linking the stochastic horizon surfaces to their equivalent grid layers.

Common applications of this method include assessing uncertainty on estimates of gross rock volume (GRV); optimizing well locations and using structure as an uncertain parameter in history matching. In this paper, we aim to assess uncertainty on estimates of connected volume and forecast production to fault position and displacement. Modelling uncertainty on fault displacement is achieved by simulating independent random variables (P-fields) on the hanging wall and footwall horizon surfaces. In order to ensure that faults do not reverse their sense of slip, uncertainty on either side of the fault is limited to a maximum of half the initial displacement. If both surfaces move in opposite directions, the maximum possible displacement is twice the initial value, and the minimum possible displacement is zero. A long correlation length can be used to simulate variations in ‘global’ fault displacement arising from seismic interpretation uncertainty. Local heterogeneities in fault zone structure (drag, relays, segmentation) can be modelled using a shorter range property; or the two properties may both be modelled.

## Fault property uncertainties

Having defined a method for simulating geometrical uncertainties, variability in fault properties must be considered. As noted above, faults in reservoirs are rarely sampled by wells; hence their properties must be estimated using available parameters. We use an approach similar to that proposed by Manzocchi *et al.* (1999), whereby we estimate fault rock thickness from fault displacement and fault permeability from clay content, which is itself approximated using conceptual models such as SGR and/or SSF. To simulate multiple realizations of fault properties, we define the coefficients and constants used in the empirical functions as uncertain variables, whose distributions are estimated from analysis of analogue databases. Here we describe in turn the data available on fault rock permeability and thickness, the data analyses undertaken and the property modelling process.

### Fault permeability data

The permeability data used in this study are from small faults in North Sea Middle Jurassic reservoirs (Fisher & Knipe 2001; Sperrevik *et al.* 2002), supplemented with core and outcrop data from various basins (Gibson 1998) and additional unpublished core data. The data of Gibson (1998) comprise 32 outcrop and core samples from sandstones in various depositional and tectonic settings. Faults are classified according to type (normal, reverse or strike slip) and displacement magnitude: micro (<10 cm); meso (10 cm-10 m) and macro (>10 m). From microstructural examination, fault rock type is further divided into clay-matrix gouge zones and deformation bands (cataclastic, solution or complex deformation bands). Permeability is measured perpendicular to the fault surface using Klinkenberg-corrected gas permeametry under a confining pressure of 20.7 MPa. There is no apparent systematic difference between permeability values measured at different scales; however there are only 2 data points for faults with throw greater than 10 m. Nine samples were removed from the database as either host rock clay content or fault permeability was not recorded, leaving 23 samples for further analysis, of which 12 have depth data recorded.

For all other data in this study, permeability is measured perpendicular to the fault surface using water permeametry (minimum measured value *c*. 0.0001 mD) under a confining pressure of 0.8 MPa. Further details of the experimental procedure and interpretation are given in Fisher & Knipe (1998, 2001) and Sperrevik *et al.* (2002). For a total of 193 samples, data are recorded for depth, host rock clay volume and host and fault rock porosity, permeability and capillary entry pressure. Host rock clay volume (*V*_{h}) is estimated from microscopic analysis of thin sections. As displacement is small, fault rock clay content (*V*_{f}) is usually taken to equal *V*_{h}, except for clay smears which are assumed to have *V*_{f}=60%. Samples are classified based on microstructural analysis into cataclasites, clay smears and phyllosilicate framework fault rocks. Thirty-one samples described as disaggregation zones, which lack discrete slip surfaces and instead represent deformation by grain reorganization, are excluded as we consider them to be unrepresentative of the fluid flow properties of the reservoir scale fault rocks we wish to estimate. This leaves 162 samples for further analysis, making a total of 185 samples with the data of Gibson (1998) (Fig. 1).

An exponential function is applied to all samples in order to compare measured permeability values at a reference (zero) effective stress, assuming hydrostatic compaction and relatively low stress sensitivity to permeability (David *et al.* 1994). Reference permeability is hereafter referred to as log *k*_{o}. The same function can subsequently be used to correct reference values to *in situ* stress conditions for a reservoir of interest.

### Fault zone thickness

To estimate the effects of faults on reservoir fluid flow, we are primarily concerned by the central zone of more or less deformed material (fault rock) bounded by two or more discrete slip surfaces, or the thickness of the slip surface itself where no fault rock is developed. Fault displacement and thickness data are available for 933 samples of fault rock ranging from mm-scale microfaults to crustal scale strike slip faults with several hundred kilometres offset. Consistent with published observations of displacement and thickness (e.g. Hull 1988; Knott *et al.* 1996; Sperrevik *et al.* 2002) we observe an approximately power-law relationship with a scaling exponent close to 1 over this size range (Fig. 2), which is interpreted to reflect the processes (abrasion, grain fragmentation) responsible for fault rock development with progressive displacement. A linear regression between log throw (*D*) and log thickness (*t*_{f}) for all samples gives the following expression:
1
where *c*=0.866, standard error=0.014; *e*= −1.62, standard error=0.02.

### Data analysis

Through deformation-induced mixing (Fisher & Knipe 2001), fault permeability in siliciclastic rocks is commonly believed to be negatively correlated with clay content. This is the basis for the empirical transformations of Manzocchi *et al.* (1999) and Sperrevik *et al.* (2002) and is borne out by the experimental data of Crawford *et al.* (2002). Permeability reduction during faulting in impure sandstones occurs by grain rotation and repacking, hence platy clay minerals are more likely to develop interconnected microtextures that reduce permeability. Fault permeability may additionally show a negative relationship with maximum burial depth due to the effects of mechanical and thermal compaction, growth of diagenetic minerals and enhanced quartz precipitation (e.g. Fisher & Knipe 2001).

A linear regression between log permeability and *V*_{f} for all samples (*n*=185) gives the following expression:
2
where Log *k*_{o}=log permeability at reference pressure *P*_{o}; *a*=−0.046, standard error=0.005; *b*=−1, standard error=0.15.

We assume that there is no further permeability decrease at values of *V*_{f} greater than 60%; this approximates the form of functions derived from laboratory deformation experiments where clay content can be more accurately controlled (e.g. Crawford *et al.* 2002) and prevents excessively low (<10^{−5} mD) permeability values from being output (Fig. 1).

For the data with a depth measurement recorded, a weak negative correlation is observed with fault permeability. However, this adds little predictive power to the existing regression and it is unknown if samples share a common depth datum. Furthermore, including a second term in the permeability equation introduces additional complications later in the modelling process; hence only clay content is kept as a significant parameter.

This equation can be used to compute fault permeability from a value of *V*_{f} estimated using a proxy property such as shale gouge ratio (Yielding *et al.* 1997). This is typically performed on a single realization of the reservoir model, using a single SGR-fault permeability function. Instead of regarding this function as fixed, we define the slope and intercept as uncertain variables, which take Gaussian distributions defined by the mean and standard error of *a* and *b* respectively. We assume that the input data are representative of the (larger) faults we wish to model and that a geometric average is appropriate to upscale fault permeability from core to reservoir grid scale. Use of a single permeability function throughout the model implies that all the faults developed under similar physico-chemical conditions; multiple functions could be used for reservoirs thought to have undergone several deformation events.

Several authors have proposed a relationship between fault displacement and permeability. Antonellini & Aydin (1994) reported that deformation bands in aeolian sandstone displayed a large (1–2 orders of magnitude) permeability decrease associated with the development of discrete slip surfaces. As this occurs at a displacement of *c*. 1 m, faults smaller than this are unlikely to be included in earth models, however this may in part account for the large scatter seen in permeability values for small faults in clean sands.

Crawford *et al.* (2002) noted a 1–2 orders of magnitude permeability decrease in sheared mixtures of quartz and kaolinite, particularly for Vclay <40%. Although total experimental slip displacement is only 3.5 mm, maximum shear strain is *c*. 20; hence the authors argue that the gouge fabric developed is comparable to that found in natural seismic scale faults. Most of the permeability decrease takes place during initial displacement (shear strain ≤1). Small faults in the datasets described above may accommodate similar strains, although the data needed to estimate this are unavailable. In response to applied shear strains of up to *c*. 10, Morrow *et al.* (1984) reported 1 order of magnitude permeability decrease for some fault gouge samples but little change in others.

If absolute fault permeability values may be comparable between the core and reservoir grid block scales, a difference in distribution is likely due to the effect of sample support, whereby measured variability decreases with increasing sample size. Core plug permeability measurements, as shown in Figure 1, should be inherently more variable than ‘bulk’ permeability values estimated from, for example, well tests. Upscaling statistics from one scale to the other requires knowledge of the spatial variability of fault permeability; however any such information is highly speculative. Anecdotal outcrop evidence (Antonellini & Aydin 1994; Faulkner & Rutter 1998; Foxford *et al.* 1998) and the results of modelling studies (Fairley *et al.* 2003; Zhurina 2003) suggest that fault permeability has a correlation length (variogram range) that is larger than a core plug (*c*. 10 cm) but smaller than a typical simulation grid block (50–100 m). To make this assumption obviates the need to model fault permeability using a correlated random field, however this could be easily achieved using the workflow outlined below. Even if no explicit correlation model is specified, as SGR is a simple average of laterally distributed clay content over a vertical distance equal to fault throw, it already contains implied vertical (= throw) and lateral (= Vclay property) correlation structures. This further supports the method of defining uncertainties on the slope and intercept of the SGR-permeability function rather than modelling the distribution of residual values.

Upscaling of the mean and standard error of the SGR-permeability function may be required depending on the assumed flow paths (and arrangement of heterogeneities) within the fault zone as being either random or ‘in parallel’, necessitating the use of the geometric (Sperrevik *et al.* 2002) or arithmetic (Manzocchi *et al.* 1999) average respectively. We assume that the geometric average is appropriate, hence Equation 2 can be used directly in the reservoir grid, however an arithmetic (or even harmonic) average could be used instead.

## Stochastic model setup

The process of simulating multiple realizations is based on a conventional stochastic geological property modelling workflow normally used to generate parameters such as facies, porosity and permeability (Caers 2005). Each parameter can be modelled as a constant; sampled randomly from a distribution; interpolated or simulated using various geostatistical algorithms, or defined as a function of another property. In this workflow, we use the values from the regression analyses above to define probability distributions for the coefficients and constants in the equations used to estimate fault permeability and thickness. Functions are defined to compute various desired parameters such as fault throw, SGR, SSF, fault permeability and thickness. Finally, instead of representing faults using transmissibility multipliers, we compute a thickness-weighted harmonic average permeability (*K**), which is stored in cells adjacent to the fault surface:
3
where *D*_{x,y}=cell size in *x* or *y* grid axis respectively; *t*_{f}=fault thickness; *K*_{x,y}=grid cell permeability in *x* or *y* grid axis respectively; *K*_{f}=fault permeability.

Inclusion of this average permeability in a flow simulation model is equivalent to applying a multiplier to the permeability of the upstream grid cell (Manzocchi *et al.* 1999, p. 57, Eq. 7; p. 59, Fig. 7j). While a less precise solution, it avoids the need to compute non-neighbour connections for each geometrical realization, which can be excessively time-consuming. Furthermore, the averaged permeability captures the effects of both sedimentary and fault heterogeneities, so can be used to make rapid estimates of connected hydrocarbon volumes through application of a threshold value. Comparative simulations show that such an averaged permeability property gives a very similar flow simulation response to transmissibility multipliers; however care is required in models where well perforations are close to or within faulted cells.

The workflow is shown schematically in Figure 3 and summarized as follows:

Define input PDFs.

Simulate grid geometry; compute fault displacement and SGR.

Simulate required uncertain properties; compute SSF, fault permeability and thickness.

Compute average permeability, apply threshold, estimate connected volume.

Repeat steps 2–4 until desired number of realizations attained.

Export realizations for fluid flow simulation.

Perform sensitivity analysis using experimental design.

## Case-study

The workflow described above has been tested on a producing field located in the Brent province of the UK Northern North Sea. The reservoir comprises a horst block structure with dip closure at its southern margin, with a central fault extending the length of the structure and several smaller faults (denoted ‘north’, ‘east’ and ‘west’) which truncate against the bounding faults (Fig. 4). Juxtaposition against downthrown reservoir sands defines an apparent hydrocarbon spill point to the north. Gas and condensate are mainly accumulated in sands of the Tarbert and Balta formations, with lesser accumulations in the Etive and upper Ness formations.

The Etive formation is interpreted as upper shoreface deposits prograding to the North. A transition from the upper Etive to Ness is characterized by tidal and shoreface deposits mixed with bay-lagoon facies (shales and coal). The Ness lithostratigraphical unit represents the deposits of a delta plain environment of floodplain shales, channels (meandering, braided and distributary) and crevasse splays. The Tarbert formation is composed of tidal and shoreface deposits overlain by transgressive Balta sands. Together, the Balta and uppermost Tarbert make up *c*. 30 m of reservoir sand, with a further 20 m in the lower Tarbert, separated by *c*. 20 m of non-reservoir shale or siltstone. A displacement of *c*. 30 m is required to completely offset the upper reservoir units. While these layers would then be partly juxtaposed against the lower Tarbert sands, the intervening shale could conceivably form a baffle or barrier to fluid flow through smearing and/or formation of low permeability clay gouge (Fig. 5).

An almost complete Balta-Tarbert succession was cored near the field crest, permitting not only sedimentological and petrophysical studies but also the identification of a number of small deformation features similar to those found within the analogue database described above. Core taken from a well drilled in an adjacent fault block also contained numerous small faults, which were additionally interpreted on borehole image logs.

Fault and horizon surfaces interpreted from seismic data were used to build a stratigraphic grid with *c*. 1 30 000 cells (*c*. 92 000 active in the flow simulator). Average cell sizes are *c*. 50 m in the I and J directions and 3 to 20 m vertically. A combination of object modelling, sequential indicator simulation and deterministic interpretation was used to model 15 sedimentary facies. Porosity, permeability, Vclay and water saturation were modelled per facies using sequential Gaussian simulation, with NTG being defined from a cut-off applied to porosity. Enhanced diagenesis with depth observed in regional wells was modelled using decline trends applied to both porosity and permeability. Vertical permeability was modelled using different *k*_{v}/*k*_{h} transformations per facies. The different modelling steps were combined in a property modelling workflow to generate multiple realizations, with one selected as the operational model for use in numerical fluid flow simulation and reserves estimation.

The model was used in a conventional black oil numerical flow simulator without modification, save for defining a number of vertical permeability barriers interpreted from well tests (Fig. 5). Functions were also defined for relative permeability and vertical well flow performance.

At the time of this study, the field had approximately two years of production history, with a single well producing gas and condensate by natural depletion. In the simulation model, the well was constrained by historical gas production rate, with bottom hole pressure (BHP) being the main history matching parameter. Production forecasts were based on a gas rate target, with defined minima for export pressure and economic flow rate. In order to match the observed BHP trend, all intra-reservoir faults were modelled as sealing to fluid flow, except for the northern part of the central fault, which was assigned a constant low transmissibility multiplier. Production forecasts using such a model have implications for a possible infill well in the second panel. The goals of applying a structural uncertainty workflow to this field were: (1) to assess the model's sensitivity to different parameters; and (2) to search for a match to production history using a geologically plausible combination of parameters.

### Geometrical uncertainties

Due to the depth of the reservoir (>3500 m), vertical seismic resolution is limited to *c*. 30 m. Fault offsets are interpreted at the top Dunlin (base reservoir) level and on a strong positive reflector which corresponds to a coal layer near the top of the Lower Tarbert succession. The overlying Base Cretaceous Unconformity (BCU) reflector is also interpreted throughout, although most faults truncate against it (Fig. 6).

Although considerable uncertainty exists in the various seismic processing and interpretation stages described above (notably migration and depth conversion), horizons are regarded as ‘fixed’ in order to focus on the influence of fault lateral position and throw. Furthermore, the lateral position of internal faults is also fixed, as shifting these surfaces may cause changes in layer juxtapositions, rendering impossible the separation of fault position and fault offset effects on connected volumes and forecast production.

Uncertainty on lateral position was interpreted for the bounding faults from the width of the zone of seismic noise within which each fault could reasonably be picked (e.g. Fig. 6). The standard deviation for P-field simulation was defined as one half of this width in either direction. Maximum fault throw uncertainty on either side of the fault was set to half the initial displacement. If both surfaces move in opposite directions, the maximum possible displacement is twice the initial value, and the minimum possible displacement is zero.

To better assess the relative importance of fault position and throw, three geometrical simulations were run: one in which only fault throw was simulated; one in which only fault position was varied, and one in which both fault throw and fault position were varied. In each case, 200 realizations were simulated, making a total of 600 possible reservoir grid geometries.

### Property uncertainties

In order to focus on the effects of fault uncertainties, grid parameters such as NTG, porosity, permeability and *S*_{w} were regarded as ‘fixed’. Uncertainties were defined on the slope and intercept of the equations used to estimate fault permeability and thickness (‘KfSlope’; ‘KfIntercept’; ‘tfSlope’; and ‘tfIntercept’) according to the regression analyses described above (Eqs 1 & 2). In order to account for the possibility of clay smear as a sealing mechanism, uncertainties were defined on two additional parameters. The value of Vclay above which cells are interpreted as being shale (‘ShaleCutoff’) was defined as a Gaussian variable with a mean of 0.5 and standard deviation of 0.08. The value of shale smear factor (SSF) above which clay smears become discontinuous (‘SmearContinuity’) was defined as a uniform variable with a minimum of 0 and maximum of 7. All faulted cells containing shale smears were assigned a Vclay value of 60%, giving a minimum fault permeability value of between 10^{−3} and 10^{−5} mD depending on the simulated values for KfSlope and KfIntercept (Fig. 1).

To analyse the relative importance of different geometrical and property uncertainties, three parameter sets were defined for each set of geometrical realizations; one in which the only uncertainty was geometrical (sets 1, 4 and 7); one with additional uncertainty on fault permeability and thickness (sets 2, 5 and 8), and one containing uncertain clay smears (sets 3, 6 and 9). In total, 1800 realizations were simulated in this way:

Fault throw variable,

*n*=200.Fault throw, permeability, thickness variable,

*n*=200.Fault throw, permeability, thickness, clay smears variable,

*n*=200.Fault position variable,

*n*=200.Fault position, permeability, thickness variable,

*n*=200.Fault position, permeability, thickness, clay smears variable,

*n*=200.Fault throw, position variable,

*n*=200.Fault throw, position, permeability, thickness variable,

*n*=200.Fault throw, position, permeability, thickness, clay smears variable,

*n*=200.

The final parameter to be defined is the threshold value applied to the average permeability *K** in order to estimate connected volumes. As with other input parameters, this may be a constant, sampled from a distribution, interpolated, simulated or derived from another property. In this case, a constant isotropic threshold of 1 mD was used; any cells with a lower *K** value are assumed to be disconnected from the producer well. This coincides with the general 1 mD cut-off used when modelling the effects of diagenesis on grid permeability with increasing burial depth. However, it is an arbitrary value that reflects the particular grid and fault permeabilities, fault thickness and cell sizes in this model; a similar value would not necessarily be appropriate for different cases. Moreover, it takes no account that the predominant fluid type is gas; which is likely to remain mobile at cell permeability values of less than 1 mD. This should be borne in mind when comparing estimated connected volumes with production forecasts based on numerical flow simulation.

### Results

The main heterogeneity is the central fault that divides the gas-bearing portion of the model into two panels; grid permeabilities generally exceed 1 mD and do not vary between realizations. The volume of gas connected to the single producer well may comprise only the eastern panel or, depending on the central fault's throw and permeability, both panels. The gross rock volume (GRV) of these panels is modified in parameter sets 4–9 by varying the position of the bounding faults. Additional complexity is introduced by the vertical permeability barriers shown in Figure 5, which separate the perforated zones of the well from gas accumulations in the stratigraphically lower Ness and Etive formations, as well as limiting the contribution from the western panel for different geometrical and fault property combinations.

A control model was defined with ‘base case’ geometry; fault properties estimated using the mean values (slope and intercept) from the data analysis and containing no clay smears. All reported volumes are normalized as a ratio of the total forecast gas production (FGPT) for this model, obtained from numerical flow simulation. Figure 7 shows plots of computed connected gas volume (CGIIP) against realization number; summary statistics on CGIIP for each of the nine parameter sets are given in Table 1.

The differences in CGIIP between the first three parameter sets are slight, suggesting that fault throw is the most influential variable for these models (Fig. 7a). For parameter set 4, where fault throw is fixed but the position of the bounding faults may vary, none of the realizations connect across to the western panel, hence variation is small (Fig. 7b). A greater spread of results is observed when fault permeability, thickness and clay smears are considered as uncertain parameters (sets 5 and 6), in which case about a quarter of the models connect to both panels. Varying both fault position and throw (sets 7–9) produces a result that compounds the two parameters' individual effects (Fig. 7c). For all geometries, adding shale smears (blue curves) is seen to reduce the mean connected gas volume, but has a limited effect on standard deviation. Figure 8 summarizes these observations in the form of cumulative distributions of connected volume for each parameter set.

### Numerical simulations

To investigate the model's dynamic sensitivity to fault heterogeneities, and to assess the usefulness of connected gas volume as a predictor of forecast gas production, realizations were sampled from parameter set 9, which includes uncertainty on fault throw and position, fault permeability, fault thickness and clay smears. Realizations were exported per decile of connected volume from Q0 to Q100, making a total of 11 simulations (Fig. 8).

Numerical fluid flow simulations were performed on the exported realizations with the same settings as for the operational model. From a crossplot of the two responses (Fig. 9, blue crosses), it appears that CGIIP is, at best, an approximate estimator of forecast gas production. Small connected gas volumes tend to underestimate total gas production, and large connected gas volumes overestimate total gas production. This can be understood by imagining a realization with only a single cell above the 1 mD threshold. Estimated connected volume will be high (both panels), while across-fault flow rate will be comparatively low. By contrast, a realization with many cells just below the 1 mD threshold would have a lower connected volume (eastern panel only) but a much greater simulated across fault flow rate.

A plot of P/Z, (pressure divided by gas compressibility factor), against cumulative production gives an estimate of connected gas in place and ultimate recovery for reservoirs producing by natural depletion. This is shown in Figure 10, (here computed using mostly flowing, rather than shut-in, BHP values). The observed data (black crosses) plot towards the lower end of the range of simulated realizations (shaded area); the model computed using the ‘base case’ geometry and properties (blue curve), which lies at the centre of this range, overestimates true connected gas in place by approximately one quarter.

### Sensitivity analysis

To more rigorously analyse the model's sensitivity to individual parameters, a two-level fractional factorial experimental design was employed. The eight uncertain parameter values were systematically varied between ‘low’ and ‘high’ levels in a planned exploration of the uncertainty space. A full factorial design, evaluating all possible combinations, would require 2^{8}=256 simulations. In order to reduce the computational time, a resolution IV screening design was selected with the generating fraction (in Box notation) I=ABCF=ABDG=ACDEH, requiring 32 experiments. This allowed estimation of the effects of all eight principal factors and thirteen two-factor interactions.

Due to the independent simulation of fault throw and fault position used to generate grid geometries, it was not possible to sample realizations directly from parameter set 7 for use in the designed experiments. Instead, geometrical simulations of fault position were first rerun for two input grids with ‘high’ and ‘low’ throw on the central fault. From each of these models, two realizations were selected to represent ‘high’ and ‘low’ GRV, resulting in four discrete grids representing the ‘end member’ geometries. Parameter values for the six fault property factors were taken from the minimum and maximum values output from the stochastic modelling stage of the workflow; these are listed in Table 2.

Three responses were chosen to estimate sensitivity coefficients for different parameters: CGIIP (estimated from the permeability cut-off prior to simulation); FGPT and simulated well bottom hole pressure (BHP) after two years of production. The latter response is the main history matching parameter for the operational model. All responses were normalized: CGIIP and FGPT to the ‘base case’ model total production and BHP to the base case model pressure after two years. For each response, sensitivities are estimated using a linear approximation of the form:
4
where *y*=response; *a*_{0}=constant; *a*_{1−n}=coefficients 1 to *n*; *x*_{1−n}=factors 1 to *n* (and their interactions).

The calculated sensitivity coefficients for all single factors and significant (at the 95% confidence interval) two-factor interactions are given in Figure 11; positive sensitivities (factor increase leads to response increase) are shown in blue and negative sensitivities (factor increase leads to response decrease) in red.

For CGIIP, the most sensitive parameter is fault throw, followed by gross rock volume (= fault position), fault permeability slope, fault permeability intercept and various two factor interactions (Fig. 11a). A similar rank of sensitivities is observed for BHP after 2 years of production, which is also sensitive to individual clay smear parameters (Fig. 11b). GRV is the most influential factor for FGPT (Fig. 11c), which reflects both the economic gas rate limit used to stop the simulation runs and the flow of gas around and through faults late in simulated field life. The greater sensitivity to fault permeability slope can be explained by the relatively high SGR values (average *c*. 35%) occurring in this case. A model with lower average SGR would likely be more sensitive to the value of fault permeability intercept.

Flowing P/Z against FGPT for the maximum and minimum experiments is shown in Figure 10 (green and red curves respectively). As expected from a systematic, planned exploration of the uncertainty space, there is a greater spread of values than for the realizations sampled from the CDF of connected volume (shaded area). Four of the models give a pressure response within a few percent of the observed data. These models all have the ‘high throw’ geometry, while the closest model with the ‘low throw’ geometry predicts a BHP 25% greater than observed. The four best-matched models also share the low value for fault permeability slope but have variable values of fault permeability intercept. Three of the models have the ‘large GRV’ structure. Forecast total gas production for the four models varies from 71 to 95% of the forecast total production for the ‘base case’ model, with three and a half years' difference in predicted field life between the smallest and largest of the four models.

Subsequently, additional pressure data were obtained during a shutdown of the field for scheduled maintenance. BHP was seen to increase at a constant rate for two months, indicating probable gas flow through the central fault from the western panel. Only one of the four best models is able to correctly reproduce the observed pressure build-up (Fig. 10, orange curve). This result helps to constrain fault properties in this field and in adjacent untested fault blocks, and to evaluate a possible infill well in the western fault panel.

## Discussion

### Model sensitivity to fault throw and position

In this study, we have only altered the position of the boundary faults, as shifting internal reservoir faults causes changes in layer juxtapositions, particularly for dipping horizons and shallow-dipping faults, as in this case. This limits the uncertainty on the size of the first panel and removes some interactions from the model, for example, a small panel with less sealing central fault might give a similar result to a larger panel with a more sealing fault. In this case, we have taken this decision to reduce ambiguity in interpreting sensitivities; however it could be argued that changes in juxtapositions are not only acceptable but desirable, if they preserve realistic geometries that could be interpreted from seismic data. For a study in which the objective was not to separate the effects of individual geometrical parameters but simply to assess the overall uncertainty on connected volumes and reserves, the position of internal faults could be varied along with that of the boundary faults.

### Fault permeability computation

Although all permeability samples were corrected to zero effective stress prior to data analysis, no correction was applied to estimated grid fault permeability values. Under initial reservoir conditions, the maximum reduction in permeability would be less than half an order of magnitude, however this effect would likely increase as the reservoir pressure declines during production. While initial permeabilities could be simply modified to reflect the effective stress conditions, these would require to be updated several times during flow simulation to reflect the decreasing reservoir fluid pressure. Alternatively, this behaviour could be approximated using simulator keywords for rock compressibility applied to the region of faulted grid cells. An overall decrease in fault permeability would shift the simulated production and BHP profiles downwards and increase the number of models able to match production history.

All fault permeabilities used in this study are absolute single-phase values applied to all fluid phases in the simulator. Two-phase properties could be modelled using pseudo relative permeability functions (Manzocchi *et al.* 2002) modified to work with averaged grid permeabilities rather than transmissibility multipliers. If the central fault zone were water-filled (Rivenæs & Dart 2002), relative permeability to gas would be far lower than the single phase permeabilities used in this case. Even if it were mainly saturated with hydrocarbons, end point relative gas permeability for the fault zone material might be two orders of magnitude lower than that of the reservoir material (Fisher 2005). For a synthetic model with a similar configuration to this case, Al-Busafi *et al.* (2005) observed an interaction between fault thickness and sensitivity to multiphase properties. For a constant fault thickness of 0.3 ft, single phase and multiphase fault properties gave a similar cumulative gas production, but results were very different for a 3 ft thick fault zone. In this case, estimated thickness for the central fault varies over a similar range (*c*. 0.1 to 1 m), hence for higher displacement, thicker parts of the fault, modelling multiphase fault properties might cause a difference in forecast gas production.

### Geometrical uncertainties

As noted above, the seismic quality leaves considerable ambiguity in fault and horizon interpretation. As well as stochastically varying fault position and throw, alternative structural interpretations are possible. Even with the existing fault pattern, additional near-seismic resolution faults could be included in the grid. Either of these methods would require building additional parallel models and repeating the study. Addition of internal faults to the existing model would likely reduce overall recovery factors and connected volumes, even for faults with relatively high permeability, low thickness and discontinuous clay smears. Faults with a displacement below seismic resolution could be modelled using stochastic object-based (Boolean) simulation, with their flow effects being captured as modifiers to cell permeability within the existing grid geometry. Modelling of discrete fault offsets in this manner would be more complex and would probably require stochastic simulation of triangulated surfaces. For the studied field case, we do not think that sub-seismic faults would have a significant impact on reservoir fluid flow, although they might be important in other situations.

In this study we have purposely regarded the overall position of the input horizons as fixed in order to focus on fault-related parameters. In reality, with control from only a single well penetration there is obviously considerable uncertainty arising from seismic processing, depth conversion and interpretation. A more complete assessment of uncertainty could be achieved by including horizon position as a variable in the geometrical simulations.

### Connected volume estimation

One of the main drawbacks of the current method is the use of a cut-off on averaged permeability to estimate connected hydrocarbon volumes, as this has been shown to be inaccurate in this case (e.g. Fig. 9 and accompanying discussion). An improved connectivity estimate would take into account the area and properties of the juxtaposed fault ‘windows’, as well as the distance from in place volumes to injector or producer wells. However, the method must also be sufficiently fast to compute that it can be applied to all realizations generated using the workflow.

An acceptable compromise between speed and accuracy might be found with streamline simulation, where pressure is computed in the three dimensional grid, as per conventional finite difference reservoir simulation, but the fluid flow equations are solved instead on one dimensional streamlines. To estimate connected volumes, it would not be necessary to run a full forward simulation if there were no planned changes to the well pattern or operating conditions (i.e. rate and/or pressure constraints). Instead, a cut-off could be applied on the initial streamline pattern using a combination of drainage time to producer wells and time-of-flight from any injector wells. Values for these cut-offs would need to be chosen in the light of reservoir engineering experience and tuned to the results of full-field simulation for the particular case. The option of using streamlines in this way is attractive as the geomodelling tool used to implement our workflow already contains a link to a commercial streamline simulator.

An improved estimation of connected volumes would give a more accurate probability distribution from which to sample realizations for use in full dynamic simulations and would help to better quantify the effect of heterogeneities on forecast production. This is the subject of ongoing research and development.

### Integrated uncertainty modelling

Finally, as the original purpose of the stochastic modelling tool is to generate multiple realizations of grid properties, these can easily be included in our workflow. This would allow investigation of the relative sensitivities of the model to sedimentary uncertainties (facies, Vclay, porosity, NTG, *S*_{w} etc.) and of their interactions with fault properties. Such an integrated approach can also be used to generate ensembles of realizations incorporating fault properties for use in assisted history matching schemes.

## Conclusions

We have outlined a workflow for modelling multiple realizations of fault geometry and properties using 3D geomodelling software. The workflow has been illustrated using a producing North Sea field in which structural uncertainties are believed to have an effect on connected volumes. The main conclusions are as follows:

Connected volume estimated from a cut-off on averaged permeability provides an approximate prediction of recoverable reserves. Accuracy could be improved by taking account of time/distance to the producer well, for example, using streamline simulation.

For connected volume and well bottom hole pressure after two years, fault throw and fault permeability are the most important factors. For forecast total production, fault position and fault permeability are the most important factors. This difference is due to the minimum economic rate constraint applied to the producer well and to the flow of gas through and around low permeability faults late in field life.

The slope of the function relating fault permeability to clay content is more important than the intercept, due to the relatively high average SGR values in this model.

Observed well pressures are at the lower end of the range of values predicted by probabilistic modelling. A better match to the observed data could be obtained by modifying fault permeabilities as a function of effective stress or fluid saturation.

A planned exploration of the uncertainty space using experimental design covers a wider range of possible pressure profiles. Four models are able to match the observed pressure after two years, but only one correctly reproduces a subsequent pressure build-up. This constrains fault properties and can be used to evaluate a possible infill well.

## Acknowledgments

This work was performed as part of A. Irving's PhD study at University College Dublin. Total E&P UK is thanked for sponsorship and for permission to publish this work. Don Medwedeff and an anonymous reviewer are thanked for their thoughtful comments, which considerably improved the paper.

- © The Geological Society of London 2010