## Abstract

**ABSTRACT** A comprehensive simulation study tested and analysed the sensitivity of dynamic connectivity in turbidite channel reservoirs to a large number of stratigraphic and engineering parameters. The study showed that subseismic shale architecture has a significant effect on reservoir connectivity. However, representing the complete spectrum of fine-scale architectural details in full-field simulation models is beyond the limits of existing computational capabilities. Previous work demonstrated that incorporating geologically based pseudo-relative permeabilities into relatively coarse full-field reservoir models renders practically intractable simulation cases tractable. We developed a methodology for generating pseudo-relative permeabilities at multiple geological scales, incorporating the effect of channel architecture and reservoir connectivity into fast simulation models.

We describe a dynamic modelling workflow that integrates geologically based pseudo-relative permeabilities into a two-stage automatic history-matching algorithm. The history-matching problem is posed as one of data conditioning in the Bayesian framework. We show the application of the workflow to a channelized turbidite reservoir in West Africa. It is demonstrated that multiple geologically consistent models that are conditioned to production data can be generated rapidly thanks to optimally coarse simulation models that capture the effect of subseismic channel architecture on recovery behaviour, and run efficiently as the forward model within a Bayesian inference framework. Proof-of-concept tests carried out using field data indicate that the history-matched models predict well-by-well future recovery response with good accuracy.

- turbidite reservoirs
- data integration
- history matching
- uncertainty
- reservoir forecasting
- inverse problems
- reservoir simulation

## INTRODUCTION

A complex network of meander belts, channel storeys, channel infill facies and bounding shales constitute the multi-scale stratigraphic architecture of channelized turbidite reservoirs (Prather *et al*. 2000; Navarre *et al*. 2002; Abreu *et al*. 2003; Deptuck *et al*. 2003; Mayall *et al*. 2006; Barton *et al*. 2010). Representing the complete spectrum of fine-scale geological detail in full-field dynamic models is beyond the limits of existing computational capabilities. Such models result in a very large number of gridblocks, require very long simulation run times, and there often exist a very large number of geological realizations to consider. On the other hand, conventional approaches that rely on judgement-based selection of a few (typically volumetric) end-point model realizations expose the forecast statistics to undersampling and bias. Knowledge about subseismic architectural parameters stems from outcrops, shallow analogues, and process-based modelling observations. Information is also derived from core and well-log data (image and dipmeter logs), which are, by nature, areally sparse. Thus, knowledge about subseismic architecture often comes with large uncertainties. Statistically rigorous quantification of these uncertainties leads to a large number of geological realizations. On the other hand, not all of the fine-scale detail may be pertinent to dynamic performance. Alpak *et al*. (2010) document the results of a comprehensive simulation study carried out using high-resolution (decimetre-scale) sector models and their stepwise-coarsened counterparts. Stratigraphic parameters that dominate the dynamic connectivity of turbidite channel reservoirs are tabulated for water injection, dry gas injection and single-phase (gas) depletion recovery mechanisms. A methodology is described in Alpak *et al*. (2010) for generating pseudo-relative permeabilities at multiple geological scales, incorporating the effect of channel architecture and reservoir connectivity into fast simulation models. The technique is coined as the ‘3D Connectivity Method’, hereafter simply referred to as 3DCM. In this paper, we describe a multi-scale assisted history-matching and forecast uncertainty quantification workflow that utilizes 3DCM as an integral component. The history-matching problem is posed as one of data conditioning in the Bayesian framework. We demonstrate the application of the workflow to an offshore West African reservoir, hereafter referred to as Reservoir A.

Contemporary automatic/assisted history-matching algorithms can be broadly classified into three main categories: (1) deterministic, gradient or sensitivity coefficient-based methods (Gomez *et al*. 2001; Brun *et al*. 2004; Gao & Reynolds 2006); (2) stochastic global search techniques (Kitanidis 1995; Schulze-Riegert *et al*. 2002; Subbey *et al*. 2004); and (3) recursive state estimation techniques (Skjervheim *et al*. 2007; Haugen *et al*. 2008). Stochastic global search techniques are simulator-independent, relatively more straightforward to implement and tend to avoid local optima. Several stochastic search algorithms – genetic algorithms (Erbas & Christie 2007), simulated annealing (Ouenes *et al*. 1993), tabu search (Yang *et al*. 2007), neighbourhood algorithm (Sambridge 1999; Stephen *et al*. 2006), modified Markov-chain Monte Carlo algorithms (Maučec *et al*. 2007; Alpak *et al*. 2009) and agent based optimization techniques, such as particle swarm optimization (Mohamed *et al*. 2009) – have been applied to history-matching problems. Christie *et al*. (2006) used neural networks in combination with the neighbourhood algorithm for history matching. Recently, a variety of commercial algorithms have been introduced in the literature to aid engineers and geoscientists in the history-matching process (Schulze-Riegert *et al*. 2002; Bustamante *et al*. 2005; Al-Shamma & Teigland 2006; Cullick *et al*. 2006). The main theme of such algorithms is time savings through automation of tedious and repetitive tasks, e.g. model generation with different parameter values (pre-processing), running simulations (processing) and evaluating simulation results (post-processing). The goal is to shift the focus, engineering knowledge and experience to the physical problem rather than to routine time-consuming tasks of insignificant value. Ensuing time savings could be invested into devising different input strategies, interpreting simulation results and introducing improvements in reaction to intermediate history-matching outcomes. Contemporary assisted history-matching techniques are, in general, of stochastic nature and naturally lend themselves to distributed computing where multiple simulations run on multiple computer nodes, thereby reducing the project turnaround time.

We formulate the history-matching problem as one of data integration in the Bayesian framework (see, for example, Christie *et al*. 2006; Maučec *et al*. 2007; Alpak *et al*. 2009). Challenges stemming from the nonlinear nature of conditioning reservoir realizations to production data have been discussed and a non-comprehensive, yet instructive, review of alternative techniques has been reported in Oliver *et al*. (2001). Instead of developing ‘yet another’ iterative mathematical engine to address the problem of conditioning reservoir models to production data, we propose a model parametrization strategy that (a) rigorously reflects the information that is known about the reservoir prior to the incorporation of production data (Curtis & Wood 2004; Wood & Curtis 2004); (b) facilitates geologically sensible updates to model parameters over the course of iterative procedures characteristically employed by nonlinear data-integration algorithms. Our workflow is developed with a special focus on channelized turbidite reservoirs. However, the potential to extend it to other depositional environments is self-evident.

The history-matching problem in channelized turbidite reservoirs, where the dynamic response of the reservoir is dominated by stratigraphic heterogeneities rather than structural discontinuities (e.g. faults), involves model parameters that govern reservoir connectivity at stratigraphic element (e.g. valley, meander belt, channel, channel infill) scales and petrophysical property values at the gridblock scale. Characteristically, the former parameters (for instance, valley drape transmissibility multiplier, meander belt shale drape coverage, etc.) often extend well beyond gridblock scale, hence, are classified as large-scale/global parameters. On the other hand, the latter parameters represent the spatial distributions of petrophysical properties (e.g. absolute permeability), thus, are suitably referred to as local parameters. These two classes of history-matching parameters lend themselves naturally to different types of inference techniques. A global inference technique is more suitable for obtaining a rapid first-stage history match for a few potentially highly influential large-scale parameters and a few other local parameters that influence the well behaviour (e.g. skin factors). An inference technique with strong hill-climbing capability (e.g. Neigbourhood Algorithm (NA)) is desirable at this stage. On the other hand, it has been shown that Ensemble Kalman Filter (EnKF) operates very efficiently for problems where (1) the number of estimated model parameters is very large (such as local gridblock-scale parameters) and (2) the model parameters demand smooth updates, e.g. the permeability field (Skjervheim *et al*. 2007; Haugen *et al*. 2008). Thus, we employ a two-stage algorithm using existing mathematical techniques but appropriately selected for the data-conditioning problem at hand at a given stage. The first-stage iterative search is performed by a custom implementation of the NA (Sambridge 1999) on an ensemble of parametric model realizations. In the second stage, production data are incorporated into reservoir models by adjusting gridblock parameters via the EnKF technique (Haugen *et al*. 2008).

Ensemble-based techniques employed in both stages of the two-stage data-conditioning workflow permit quantification of forecast and estimated model parameter uncertainties. We show the field application of our data conditioning and reservoir forecasting workflow on a dataset from a West African channelized turbidite reservoir.

## MULTI-SCALE TURBIDITE CHANNEL ARCHITECTURE

Model objects range from meander belt (several hundreds to thousands of metres wide) to channel storeys (a few hundreds of metres wide), to channel infill beds (decimetres thick) and to petrophysical property scales. Methods have also been proposed for recognizing and quantifying the distribution of important stratigraphic parameters using core and dipmeter data (Barton *et al*. 2010).

Key stratigraphic elements in a typical turbidite channel reservoir include (1) meander belts, (2) channel storeys, (3) channel infills and (4) bed-scale petrophysical properties (Fig. 1). The largest-scale elements, informally referred to as the first-order elements, are the meander belts (Prather *et al*. 2000; Abreu *et al*. 2003). These elements are similar to channel complexes and are comprised of a set of vertically and laterally offset sinuous channel storeys and associated overbank deposits bounded at the base by an erosional disconformity. Offset of the channel storeys was adjusted to honour swing- and sweep-type relationships observed within near-seafloor analogues (Deptuck *et al*. 2003). Critical shale elements that may influence connectivity at this scale (between meander belts) are mudstone drapes deposited along the base and margins of the meander belt. Internal to the meander belts are channel storeys (also simply referred to as channels) and inter-channel elements, such as inner-levee and crevasse-splay deposits. Intra-channel shales, channel-base drapes and mud-filled channels govern connectivity between channel storey elements (Fig. 1). Within the channel storey elements, connectivity depends on the infill geometry (convergent, layered, accretionary, etc.) and the frequency and lateral extent of intra-channel shales (third-order elements). Finally, petrophysical properties of various sand, silt, and shale facies constitute the fourth-order elements.

## 3D CONNECTIVITY METHOD

Dynamic behaviour of subseismic turbidite channel architecture in coarse field-scale reservoir models can be represented by means of connectivity factors implemented in terms of effective properties (Alpak *et al*. 2010). The 3D connectivity method differs from the solution error model described in O'Sullivan & Christie (2005) where error models are introduced into the misfit calculation within a history-matching context. On the other hand, our error modelling methodology, i.e. the 3D connectivity method, represents the coarse-to-fine corrections (the 3D connectivity factors) in the form of pseudo-relative permeability functions (Alpak *et al*. 2010). Thus, it is directly applicable to dynamic forecasting studies in the absence of historical production measurements. On the other hand, pseudo-relative permeability functions can be conditioned to production data when such data are present, as described in this paper.

Connectivity factors are derived from representative full-detail and coarse-scale sector model simulations in the absence of production data. High-resolution sector model simulations are performed to capture an accurate representation of the geological detail. Coarse-scale and high-resolution sector model simulations are used together for quantifying the impact of channel architecture in terms of 3D connectivity factors. The 3D connectivity factor (*CF*) quantifies the influence of model parameters at different model complexities. It is defined as the ratio of the fine-scale recovery factor (*RF*) to the coarse-scale recovery factor at a given point in time during production,

Connectivity factors are computed for two crucial points in production time. The first type of connectivity factor is calculated using fine- and coarse-scale sector model recovery factors at their respective breakthrough times (t = bt). The resulting connectivity factor is defined as the breakthrough-time connectivity factor, *CF*_{bt}. The second type of connectivity factor is calculated using fine- and coarse-scale sector model recovery factors at the time of ultimate recovery, also called the trapped-hydrocarbon time (t = th). Typically, a sufficiently late abandonment time is selected to quantify the hydrocarbon-trapping effect of stratigraphic architecture. The ensuing connectivity factor is coined as the trapped-hydrocarbon connectivity factor, *CF*_{th}. The stronger the hydrocarbon-trapping effect of channel architecture on recovery, the closer the connectivity factor value is to 0.0. For reservoirs where the imprint of channel architecture is weak, the connectivity factor value is closer to 1.0. This observation is valid for both types of connectivity factors.

A large number of simulations (∼31 200) were carried out on channelized turbidite models at multiple geological resolutions (Alpak *et al*. 2010). Figure 2 shows cross-sections from a fine, outcrop-scale parent model and the coarse-scale offspring models. From appraisal to development phase, a typical deep-water field project would use models with a gridblock resolution from 20 000 to 1 000 000 cells, similar to the coarse-scale L5–L1 models in Figure 2. L0 models contain the full, ‘ground-truth’ architectural detail. The models were coarsened by merging the boundary layers and the architecture together through upscaling. In L1 models the layers of channel infill architecture are upscaled into a single channel infill layer. The channel infill layers are merged with their respective channel-base drapes in L2 models, the channels within a meander belt are combined into a single-unit in L3 models, the meander belt base drapes are lumped with their respective meander belts in L4 models and all meander belt layers are upscaled into a single-layer in L5 models. The simulation study described in Alpak *et al*. (2010) demonstrates that the models of the full stratigraphic detail of the turbidite channel architecture (L0-level sector models) have two main subsurface flow characteristics in reference to the coarser representations of the stratigraphic architecture (L1-, L2-, L3-, L4- and L5-level models): (1) increased trapped-hydrocarbon volume due to stratigraphic dead-ends and baffles to flow, e.g. shale drapes; (2) early injected fluid breakthrough. The latter conclusion is valid for displacement-type recovery mechanisms such as water injection and dry gas injection. Channelized turbidite reservoirs constitute a suitable and relevant first target for the 3D connectivity workflow as the architectural detail is relatively straightforward and well understood to represent with a parametric, geologically consistent integrated reservoir modelling approach.

Predictive 3D connectivity factor functions, also called proxy functions, are developed using the outcome of the largenumber of simulations carried out for water injection, dry gas injection and single-phase (gas) depletion recovery mechanisms. Proxy functions (Lechner & Zangl 2005) are closed-form equations derived by training multilayer perceptron (MLP) neural networks using a back-propagation algorithm (Haykin 1999) on connectivity factors obtained via flow simulations on fine- and coarse-scale sector models. We developed a separate proxy function for every recovery mechanism mentioned above and for every coarse modelling scale (L1–L5) in Figure 2. Proxy functions are embedded into a spreadsheet that links to a Monte Carlo engine (Palisade Corporation 2008) and calculates connectivity factor statistical distributions based on input geological parameter statistical distributions. Statistical distributions of meander belt net-to-gross, meander belt degree of amalgamation, meander belt width, meander belt shale drape coverage, channel (storey) shale drape coverage, channel infill shale drape coverage, channel infill shale frequency, and frequency of mud-plugs constitute the input. The spreadsheet also admits an average well-spacing value as input for the water/gas injection recovery mechanisms.

In the absence of production data, the procedure is to apply Monte Carlo simulation technique directly to reservoir simulation by estimating probability density functions for high-impact parameters that govern connectivity. The most influential stratigraphic architecture parameters that govern the 3D connectivity behaviour (in channelized turbidite reservoirs) are listed above and were determined in a previous study (Alpak *et al*. 2010). Findings of this study are restricted to relatively light (20–40° API) hydrocarbon-bearing turbidite channels with sand permeabilities ranging from 250 mD to 4000 mD. Properties of Reservoir A (∼30° API oil and ∼2000 mD sand permeability) conform to the 20–40° API oil and 250–4000 mD sand permeability ranges investigated in detail in Alpak *et al*. (2010). Moreover, the reservoir is undergoing water injection from the start-up date consistent with the water injection recovery mechanism simulations analysed in Alpak *et al*. (2010).

It is worthwhile emphasizing the difference between connectivity factor and recovery factor. Recovery factor behaviouris dominated by a much larger set of parameters than 3D connectivity factor, including PVT, relative permeabilites, etc. The 3D connectivity factor should be perceived as a ‘correction factor’ for the effect of unmodelled stratigraphic detail on the recovery factor. Therefore, it is sensitive to a relatively more restricted set of parameters. Effects of many parameters that influence the recovery factor cancel out in the calculation of the connectivity factor. This is because they are already accounted for in the coarse-scale model. The 3D connectivity factor depends strongly on the level of geological detail in the coarse modelling scale, and the recovery mechanism (water injection, gas injection, depletion, etc.).

When production data are present, connectivity factor distributions computed with the above procedure are imposed as prior information to our ensemble-based data-integration algorithm along with the prior distributions of other model parameters. In other words, connectivity factors are history matched along with other relevant model parameters and the history-matching process takes advantage of the availability of prior information about connectivity factors. As a matter of fact, our ensemble-based data-integration algorithm is formulated to perform a Bayesian inference. Both NA (Stage I) and EnKF (Stage II) are Bayesian data-integration techniques, not just algorithms for identifying and refining the search of good models. The data-integration process updates the prior model parameter distributions with information extracted from the production data. Typically, one obtains narrower model parameter distributions compared to the prior distributions subsequent to the data-integration process. The process incorporates dynamic reservoir connectivity information inferred from production data into updated 3D connectivity factor distributions. The next section describes how connectivity factors are translated into pseudo-relative permeability functions that can be used in the context of reservoir simulation and, hence, in history matching.

## PSEUDOIZATION ALGORITHM

Field-scale models need to encapsulate special provisions to capture the dominant effects of (unmodelled) channel architecture and shale drapes on flow. Pseudo-relative permeability functions are modified versions of rock relative permeability functions. They introduce physically-based adjustments to the relative speeds of flowing fluid phases to account for the effect of missing elements of the channel architecture. An algorithm is developed to transform rock curves into pseudo-relative permeability functions with the help of connectivity factors. Multiple reservoir regions and/or reservoirs with different connectivity characteristics (as warranted by varying geological characteristics) can be combined and simulated within a single run using region-by-region application of connectivity factors via pseudo-relative permeabilities. These region-by-region effective properties serve as localized agents of coarse-to-fine connectivity correction for stratigraphically complex reservoirs.

The pseudoization technique differs from traditional deterministic two-phase upscaling techniques that require construction of a fine-scale model. In field applications many of the key geological parameters, especially the ones that relate to subseismic architecture, can only be known statistically. The 3D connectivity workflow eliminates the need for building models at the detail-level of subseismic scale, a scale that cannot be resolved deterministically. However, the requirement is a rapid functionality that accurately predicts connectivity factors, which serve as correctors to the global flow for the effect of unmodelled elements of reservoir architecture. In our workflow, the proxy functions compute connectivity factors without resorting to case-specific simulations. The main restriction is that the stratigraphic characteristics of the reservoir of interest have to be contained within the population of generic sector models that were simulated in developing these functions. Connectivity factors can also be estimated from history matching. However, inherent non-uniqueness of the history-matching problem can induce significant uncertainty to the estimation. Therefore, the best practice is to use all available geological data to compute the prior distributions of connectivity factors following the workflow described above. Subsequently, one should use an ensemble-based production data-integration algorithm to find posterior connectivity factor distributions that honour dynamic measurements. If the production data are informative (which they almost always are), the posterior distributions are expected to be narrower.

As described earlier, simulations carried out at multiple geological resolutions demonstrate that channel architecture has two main effects on subsurface flow: (1) increased trapped-hydrocarbon volume due to stratigraphic dead-ends and baffles to flow, e.g. shale drapes; (2) early injected fluid breakthrough. The latter conclusion is valid for displacement-type recovery mechanisms, such as water injection and dry gas injection.

For displacement-type recovery mechanisms, Effect #1 is incorporated into coarse-scale models by reducing the movable oil in the reservoir (Fig. 3) through increasing the residual oil saturation via

Here, *S*_{or} and *S*_{k} represents irreducible water saturation, *S*_{wi}, for oil–water curves and residual gas saturation, *S*_{rg}, for gas–oil curves. *CF*_{th} is the trapped-hydrocarbon connectivity factor. As *CF*_{th} approaches 1.0, the amount of correction for the trapped hydrocarbons goes to 0.0. Severely disconnected reservoirs, on the other hand, require a connectivity factor significantly smaller than 1.0.

Decreasing the curvature of the water relative permeability function (*e*_{w}) accounts for the early breakthrough effect (Effect #2) (Fig. 3). The curvature reduction is linked to the breakthrough-time connectivity factor (*CF*_{bt}) by solving Buckley-Leverett waterflooding equations of motion on notionally homogeneized coarse- and fine-scales (Alpak *et al*. 2010). The same pseudoization technique is applied to gas–oil relative permeabilities by replacing the water phase with the gas phase in the workflow. *CF*_{th} acts on the depletion-type recovery mechanisms as a post-simulation tuning factor in the current implementation.

### Discussion

Previous research (see, for example, Okano *et al*. (2006)) introduced pseudoization techniques to quantify the dynamic forecast uncertainty due to simple subgrid permeability heterogeneity in reservoir models. Our approach differs from these techniques in that the primary focus is modelling the dynamic impact of the missing (unmodelled) stratigraphic architecture in coarse-scale models via pseudo-relative permeability functions. We assume that the deterministic fine-scale model and its underlying architecture (in terms of grid and petrophysical properties) is unknown but information about the dynamic response discrepancy and, hence, about the necessary correction to reconcile fine- and coarse-scale stratigraphic model responses are available in a statistical form, namely, in the form of the convolution of stratigraphic parameter statistical distributions with the 3D connectivity proxy functions. We do not attempt to model the discrepancy (error) between fine- and coarse-grid models where the coarse-grid model is derived from the direct upscaling/upgridding of a deterministic fine-grid model with specific well locations. Correction data are embedded to proxy functions that are derived from a large set of fine- and coarse-scale simulations with varying stratigraphic characteristics and well-placement configurations. We apply a methodology to convert the parameters of the 3D connectivity method, which encapsulate the effects of missing stratigraphic architecture (namely, the 3D connectivity factors) to geologically based pseudo-relative permeability functions.

The advantage of the 3D connectivity workflow is that it does not require the construction of field-specific full-detail (L0-level) full-field geomodels as the information necessary to correct coarse-scale full-field models is embedded into the proxy functions. On the other hand, our workflow demands a large initial simulation investment for a given depositional environment and recovery mechanism combination. The simulation investment is necessary to develop accurate 3D connectivity proxy functions and involves a number of generic full-detail (L0) and coarse-scale (e.g. L1–L5) sector model simulations.

As an important difference between the workflow described in Alpak *et al*. (2010) and this paper, we do not use the pure geologically based pseudo-relative permeability functions directly in the dynamic modelling workflow described here. Rather, we first condition the 3D connectivity factors (and, in turn, the pseudo-relative permeability functions) to production data. We carry out this process within the framework of a two-stage Bayesian inference algorithm that will be described in the subsequent section. We utilize the pure geologically based 3D connectivity factors to impose physical bounds to the inference process in our workflow.

## DATA-INTEGRATION ALGORITHM

The two-stage data-integration algorithm (Fig. 4a) is formulated using the Bayesian inference framework, which constitutes a statistically consistent way of updating the information about a given set of models subject to observed data. In this framework, one starts with a set of information about the details of the dynamic model. Data used to form this prior knowledge may come from the geological description of the reservoir, seismic measurements, core samples, well-log measurements and pressure-transient analysis, and other types of data used in constructing a static reservoir model. Using this prior information, one forms a way of parametrizing the reservoir description and assigns a set of probabilities to these parameters.

The next step involves sampling a number of probable reservoir descriptions from the prior model and assessing how well the predictions made using these reservoir descriptions fit the observed data. Models that render a good fit are assigned higher probabilities derived from the Bayes' rule. Bayes' theorem provides a formal methodology to update the knowledge about the model parameter probabilities, given information in the form of observed data (Tarantola 2005). This process is called sampling the parameters of the posterior distribution. Bayes' theorem represents the probability of occurrence of the model parameters, **m**, given the measured data, **d** (i.e. the posterior distribution) as proportional to the product of a likelihood function and the prior probability distribution of the reservoir description via

where *p*(**m**|**d**), *p*(**d**|**m**) and *p*(**m**) denote posterior probability density function, likelihood function and prior probability density function, respectively. The integral expression in the denominator of equation (3) represents the probability associated with the measured data. It is independent of the model parameters and, in general, treated as a normalizing term. The Bayesian approach to data integration permits the use of arbitrary probability distributions and arbitrary measures of uncertainty, not only Gaussian distribution and variance as the measure of uncertainty. Within the Bayesian framework there exists a mechanism for the rejection of unsuitable model proposals, focusing the search on appropriate models that honour prior information. This way the analysis is extended to higher levels of interpretation as opposed to a mere minimization of an objective function (Hansen *et al*. 1997).

The Bayesian formulation of the likelihood function is defined as the conditional probability distribution function of the measurement data, given the model parameters expressed as

where **C _{d}** represents the data covariance matrix,

**d**denotes the measurement data vector of size

*N*

_{d},

**m**is the vector of model parameters of size

*N*

_{m}. The nonlinear function

*f*(

**m**) represents the forward modelling operator which maps the model parameters to the data space and thereby yields the simulated model response (simulator). The posterior probability density function is proportional to

where

From equations (5) and (6) the vector of model parameters, **m**, that minimizes the objective function, *O*(**m**), maximizes the posterior probability density function, *p*(**m**|**d**), and corresponds to the most probable estimate also known as the maximum *a* *posteriori* (MAP) estimate. Bayesian inference approach updates the statistical distributions of model parameters through constraining the posterior probability density function with the prior probability density function. The severely non-unique nature of history matching as an inverse problem makes the regularization of the solution, with all available prior information about the model, essential. Appropriate selection of the prior probability density function plays a significant role in this process. On the other hand, the typical strongly nonlinear nature of the relationship between the model parameters and the dynamic response entails an iterative solution procedure.

### Stage I: Neighbourhood Algorithm (NA)

The Neighbourhood Algorithm is a stochastic technique originally developed by Sambridge (1999) for solving a nonlinear inverse problem in seismology. The algorithm first identifies regions in the parameter space that yield a good match to the observed or measured data. It then preferentially samples the parameter space in a guided way and generates new models in the good-fit regions. NA uses Voronoi cells to represent the parameter space.

Given *n* points in a plane, their Voronoi diagram divides the plane according to the nearest-neighbour rule, namely, each point is associated with the region of the plane closest to it. The distance metric is usually Euclidean, and each Voronoi polygon in the Voronoi diagram is defined by a generator point and a nearest-neighbour region. Voronoi cells are such that the sizes of the cells are inversely proportional to the density of the points in that region, and the cells are non-overlapping. Each Voronoi polygon can be regarded as a closed set where each cell is simply the nearest-neighbour region about the generator point. The Voronoi cells possess very attractive geometrical and computational properties (Okabe *et al*. 1992), which make them suitable for the task of history matching. Specifically, the neighbourhood approximation property of Voronoi diagrams provides a simple way of performing non-smooth interpolation of an irregular distribution of points in multi-dimensional space. Thus, if one represents sought model parameters by a Voronoi diagram, the problem of generating multiple model realizations translates into one of sampling a multi-dimensional parameter space defined by the Voronoi cells. By exploiting the interpolatory properties of Voronoi cells, new models can be generated and concentrated in specific regions of the parameter space using a prescribed rule. This constitutes the basic principle underlying the sampling process.

NA differs from other stochastic algorithms in the following ways. Rather than seeking the global optimum, the objective is to generate an ensemble of data-fitting models. By constructing an approximate misfit surface for which the forward problem has been solved, the algorithm samples the parameter space by performing nonlinear interpolation in multi-dimensional space, exploiting the neighbourhood property of the Voronoi cells. The sampling is done in a guided way such that new realizations of the model are concentrated in regions of parameter space that yield good fit to the observed data.

NA iteratively generates multiple data-fitting models in parameter space in the following way. First, an initial set of *n*_{si} models is generated randomly in the search space by a Sobol–Antonov–Saleev quasi-random generator (Press *et al*. 1992). For each model, the forward problem is solved and the relevant misfit value, *O*, is obtained. In the second step, *n*_{r} best models among the most recently generated *n*_{s} models are determined. The *n*_{r} models are chosen on the basis of how well their forward solutions match the measured data. Finally, a new set of *n*_{s} models are generated via uniform random walk in the Voronoi cell of each of the *n*_{r} chosen models. The algorithm returns to the second step and the process is repeated until it reaches the user-defined number of total iterations, *maxit*. At each iteration, *n*_{r} cells are resampled and in each Voronoi cell, *n*_{s}/*n*_{r} models are generated. The performance of the algorithm depends on the ratio *n*_{s}/*n*_{r} rather than on individual tuning parameters. The ratio *n*_{s}/*n*_{r} plays an analogous role to the inverse of the temperature in simulated annealing (Sen & Stoffa 1995). Although *n*_{si}, *n*_{s}, *n*_{r} and *maxit* constitute the entire set of tuning parameters of NA, *n*_{s} and *n*_{r} govern the explorative (global) vs. exploitative (local) balance of the search characteristic, respectively. A total of *n*_{si} + *n*_{s} × *maxit* models are generated by the algorithm. Generating new *n*_{s}/*n*_{r} cells within a Voronoi cell is accomplished by a Gibbs sampler (Geman & Geman 1984).

The philosophy behind the algorithm is that the misfit of each of the previous models is representative of the region of its neighbourhood, defined by its Voronoi cell. Therefore, at each iteration, new samples are concentrated in the neighbourhoods surrounding the better data-fitting models. The algorithm exploits information obtained in all previous models to selectively sample parameter space regions which yield a better fit to the observed data. This approach also overcomes one of the principal weaknesses of stochastic algorithms: slow convergence.

An in-house implementation of NA is made by closely following Sambridge (1999) (Fig. 4a, b). An automatic servo mechanism is implemented to the algorithm similar to the one described in Erbas & Christie (2006). It restarts the sampling after a given number of iterations, if no further improvement is observed in the neighbourhoods surrounding the better data-fitting models. The servo mechanism forces the sampling algorithm to jump out of local optima and to either explore uncharted sectors of the uncertainty space or generate a new variant of the uncertainty space and explore it.

Typical values for the NA parameters employed in the Stage I history matching of the field case are as follows: *n*_{si} = 100, *n*_{s} = 30, *n*_{r} = 15 and *maxit* = 8 to 15. An exploratory approach is taken in Stage I where *n*_{s}/*n*_{r} = 2. The algorithm continues until a jump is recommended by the servo or *maxit* is reached. No more than five jump recommendations by the servo are allowed. The servo uses the same objective function as the NA algorithm. If the improvement in the objective function is less than 0.1–2% (depending on the case) for three consecutive iterations, then the sampling process is restarted from a different initial ensemble of model parameters. A test process is implemented to ensure that the new ensemble is sufficiently different from all previously tested ensembles. The servo parameters documented above are rather heuristic. The sensitivity of the servo parameters to sampling efficiency is not tested extensively. Thus, there may be other tuning parameter selections that could further improve the sampling efficiency. The NA algorithm takes advantage of distributed parallel computing. Forward runs are distributed to multiple processors at a given iteration level. The algorithm is synchronized at the end of every iteration. Then, the algorithm advances to the next iteration level, where forward runs are again distributed to multiple processors.

### Stage II: Ensemble Kalman Filter (EnKF) technique

The reservoir model state and history-matching parameters are updated sequentially in time in EnKF, based on the information contained in pressure and rate measurements from production wells (Haugen *et al*. 2008). EnKF can be interpreted as the best linear unbiased estimator of a posterior model, given a prior model, measurement data and the prior and measurement statistics in the form of covariance information. In that interpretation, ‘best’ refers to a posterior estimate with minimum variance. This corresponds to the fact that an updated ensemble represents a sample of a second-order approximation of the posterior joint probability density function of the parameters. In the second stage of the data-integration algorithm, production data are incorporated to reservoir models by adjusting gridblock parameters (e.g. permeability) via the EnKF technique (Fig. 4a). Stage I history-match parameters remain unchanged in the second stage. *A priori* information is enforced in terms of gridblock property variograms. The number of estimated model parameters is typically in the order of ten to hundreds of thousands in the second stage. An ensemble size of 100 is used in the field-case applications described in this paper. Similar to NA, EnKF implementation also takes advantage of distributed parallel computing.

## FIELD APPLICATION

A turbidite channel reservoir draped across a flank of an anticline is selected for the field test. Reservoir A is saturated with light (~30° API) oil with initial reservoir pressure above bubble-point pressure. No free gas cap and no aquifer influx exist to support production. The reservoir is subject to waterflooding from the start-up date, with the primary objective of maintaining pressure such that the reservoir pressure will not decrease below bubble-point pressure, and also to displace oil to the producers. The field application involved detailed prior work in terms of quantitative interpretation of seismic, core and well-log data, leading to case-specific geological parameter statistical distributions that are fed to the 3D connectivity proxy functions.

A distinction is made in the dynamic model between major depositional zones ( Fig. 5). Phase 1 represents the external depositional zone with amalgamated channels. Phase 2 embodies the valley that cuts through Phase 1 and encompasses a number of relatively less amalgamated meander belts. Highly amalgamated channels of Phase 1 are modelled as a layered system. Interpreted meander belts have been explicitly modelled within Phase 2. The Phase 1 and Phase 2 distinction is embedded in the dynamic model with a block-by-block (depositional) phase indicator map. The distinction is necessary because different pseudo-relative permeability functions are assigned to Phase 1 and Phase 2 zones as warranted by their significantly different stratigraphic characteristics and the level of resolution used to represent them in the static model.

### Objective function

The objective function is constructed by using bottom-hole pressure (BHP) data at the injectors, and tubing-head pressure (THP) and water rate measurements at the producers. In other words, BHP data at the injectors and THP and water rate measurements at the producers are history matched, while water injection rates (injectors) and oil production rates (producers) are enforced as constraints that are fully honoured by the simulator. Since the reservoir is operated above bubble-point pressure and the PVT model is relatively simple and reliable (water and light oil system), the gas rate at surface conditions correlates with the enforced oil rate and honours the historical gas production measurements very well.

### Fault seal modelling

A functionality is developed that calculates the distribution of fault seal factors for all gridblock connections automatically based on physical fault properties. Fault properties are determined either by correlations or they can be specified based on field-specific datasets. In the case of Reservoir A, a correlation is employed similar to the one described in Jolley *et al*. (2007) but customized for turbidite reservoirs (van der Vlugt *et al*. 2006). The geological model constructed for Reservoir A is mildly upscaled in the vertical direction. We will elaborate the upscaling process in a later section. With regards to fault modelling in an upscaled model, the permeability and thickness of the fault zones are calculated pre-upscaling on the detailed grid of the static model. These fault properties are upscaled along with all other reservoir properties that are defined in the static model. The upscaled reservoir and fault properties are used to calculate transmissibility multipliers (or fault seal factors) in the simulation grid (post-upscaling). Fault seal modelling is included in the geological pre-processing step of the upscaler but its associated parameters are not modified over the course of history matching. Since the faults are predominantly outside the zone undergoing active displacement, the misfit function of the history match did not exhibit noteworthy sensitivity to the fault seal parameters.

### Shale drape modelling

Shale drapes can have a significant impact on the recovery efficiency of a channelized turbidite reservoir and, therefore, need to be modelled when present in significant quantities (Barton *et al*. 2004; 2010; Alpak *et al*. 2010). A pre-upscaling functionality was developed that incorporates bounding shales (shale drapes) for channel infills, channels, meander belts, lobe infills, lobes and lobe complexes into a geological model in terms of transmissibility surfaces (van der Vlugt *et al*. 2006). The functionality operates in the upscaler subsequent to the loading of the static model. The function introduces shale drapes into the static model with an object-based statistical process similar to the one described in Stephen *et al*. (2001). The function is designed to operate on three-dimensional cornerpoint grids (Ponting 1989) with a potentially large number of pinched-out gridlayers. The function abides by process-based erosion rules that are designed to resolve potential double-drape issues, e.g. when lobe tops and channel bases coincide at the same surface.

The shale draping functionality creates shale over the entire surface of the reservoir objects (basal surfaces of the meander belts and the Phase 2 valley in the case of Reservoir A) and then (if desired by the user) punches elliptical-shaped holes into the shale surface until the required coverage is reached. If warranted by hard data, hole locations can be guided by a probability cube. For example, they can be concentrated around the edges of basal surfaces such that the drapes appear at channel margins. Meander belt shale drape coverage in Phase 2 and the shale drape seal factor at the Phase 1–Phase 2 interface are treated as uncertain model parameters (within plausible bounds) over the course of the iterative data-integration process.

The shale draping process is part of the iterative forward modelling inside the automatic history-matching workflow. Since the shale draping functionality is implemented pre-upscaling inside the geological front-end/upscaling module of the reservoir simulator, the model parameter updates are carried out in the geological model and the shale draping and upscaling process is repeated prior to dynamic simulation for every model realization proposed by the data-integration algorithm. Since the model is upscaled only mildly and the shale draping functionality is carefully implemented to rapidly converge to the desired coverage levels, the overhead associated with the upscaling process is negligible (in the order of seconds) compared to dynamic modelling (in the order of 0.5 to 1.5 hours).

### Representing the subseismic architecture

Scoping simulations carried out on high-resolution sector models of Reservoir A indicate that the production response is affected by the subseismic channel architecture. However, the reservoir-scale static model represents the architectural detail at seismic resolution. The effect of subseismic stratigraphic architecture is modelled through connectivity factors in the dynamic model, as described earlier. Connectivity factors are handled as model uncertainties in the data-integration process.

There exists a difference between geological resolutions of Phase 1 and Phase 2 deposits in the dynamic model (Fig. 5). Amalgamated channels of Phase 1 are modelled as a layered system in the static model. Phase 2 deposits, on the other hand, are represented with explicitly modelled meander belts.

Characteristics of the fine-scale channelized turbidite architecture are described previously in the section on multi-scale turbidite channel architecture and an example vertical cross-section is shown in Figure 2 (fine-scale model [L0]). The Phase 1 zone of the model can be classified as an L4 level coarse model (Fig. 2), where a number of highly amalgamated channels are resolved with a limited number of grid-layers to satisfy the efficiency requirements of iterative automatic history matching. Although the areal size of the grid (25 × 25 m) ( Fig. 6b) is sufficiently fine to capture the channel elements with reasonable accuracy, one has to place a large number of explicit channel elements to the model to honour the highly amalgamated nature of the Phase 1 zone. It is important to note that: (1) specific locations of the channels are invisible in the seismic data; (2) time-stratigraphic surfaces between channel elements have to be preserved in the geological model, such that the shale drapes, which often appear on these surfaces, can be placed onto the geomodel in a stratigraphically consistent fashion. In our experience, stochastic channel-placement algorithms often fail to capture the time-stratigraphic surfaces/relationships with a desired level of geological realism. On the other hand, the geologically desirable but time-consuming surface-based channel modelling approach often results in very large geological models, which cannot be simulated directly. Despite the fact that such large geomodels can be constructed with considerable effort, it is important to note that conventional local single-phase upscaling techniques tend to obliterate the dynamic effects of laterally continuous shale drapes, rendering upscaled models overconnected. Here, one has to remember that shale drapes not only constitute barriers/baffles to flow in the vertical direction but also in the horizontal direction. Their extent reaches beyond the size of coarse gridblocks, into which flow properties are homogenized, violating the assumptions of effective medium theory (Renard & de Marsily 1997). Due to the above reasons, instead of modelling amalgamated channels of Phase 1 explicitly, the dynamic effects of the stratigraphic architecture within and between channel elements are represented by use of effective pseudo-relative permeability functions, which can be interpreted as a product of a global two-phase upscaling approach with a stochastic component (Alpak *et al*. 2010).

On the other hand, the Phase 2 zone of the model is at the L3 level of coarseness (Fig. 2). Phase 2 deposits are less amalgamated compared to Phase 1 deposits. Three major meander belt-scale elements (typically 500–1100 m in width and 20–50 m in max. thickness) are illuminated by seismic data. Meander belts are explicitly modelled by retaining their bounding surfaces (a product of a slice-by-slice seismic interpretation) in the geomodel. The areal grid of size 25 × 25 m is clearly sufficiently fine for representing the meander belts accurately. The vertical grid follows the meander belt surfaces and, thus, is of variable size. Shale drapes along meander belt surfaces are modelled in terms of transmissibility multipliers. Explicitly including the Phase 2 meander belt surfaces implies that the valley surface between Phase 1 and Phase 2 deposits is also explicitly modelled in the geomodel. The valley drape between Phase 1 and Phase 2 zones is modelled in terms of a transmissibility multiplier. Internal architecture of the meander belts is modelled by use of infill layers populated with facies-based effective absolute permeability values obtained through upscaling of core and well-log data. As described in the subsequent section, meander belt infill layers of the static model are slightly upscaled to fewer layers in the dynamic model to increase the computational efficiency of the forward model. In addition, effective pseudo-relative permeability functions are employed in Phase 2 deposits to represent the two-phase dynamic effects of the unmodelled stratigraphic architecture within meander belts.

Core and well-log data quantitatively indicate that subseismic geological characteristics relevant to connectivity factor calculation (such as meander belt net-to-gross, infill drape coverage, etc.) differ for Phase 1 and Phase 2 deposits. Thus, different sets of connectivity factors are employed for Phase 1 and Phase 2. In fact, two connectivity factors are necessary for a given depositional zone: one that accounts for the hydrocarbon trapping – the trapped-hydrocarbon connectivity factor, *CF*_{th} – and the other that represents the breakthrough effects also known as the breakthrough-time connectivity factor, *CF*_{bt}. Computations performed with the proxy function to determine the search ranges for connectivity factors demonstrate that *CF*_{bt} is smaller in magnitude compared to *CF*_{th} for the channel architecture in Reservoir A. In the data-integration workflow, we link *CF*_{bt} to *CF*_{th}, with a multiplier, *c* such that *CF*_{bt} = *c × CF*_{th} and 0 ≤ *c* ≤ 1. Hereafter, we will simply refer to *CF*_{th} as the connectivity factor and *c* as the connectivity factor multiplier. Probability density functions (PDFs) for connectivity factors and connectivity factor multipliers for Phase 1 and Phase 2 are determined by use of appropriate proxy functions. These PDFs are fed to the data-integration algorithm as prior PDFs.

The pseudoization algorithm of the 3D connectivity technique is implemented as a pre-processing code into the core of the simulator. This way, if desired, connectivity factors can be applied to rock relative permeability functions on a gridblock-by-gridblock basis. It is, thus, very straightforward to apply them on a regional basis as long as there exists an array that signifies different model regions in the simulation model. The depositional phase indicator map shown in Figure 5 serves for this purpose.

The gridblocks of the geological model are populated with the net-to-gross property, which is a product of seismic inversion. Net-to-gross inversion was carried out on a voxelized seismic-grid that honours the interpreted reservoir element (valley, meander belt) surfaces. Subsequently, inversion results are resampled into the cornerpoint grid of the geological model. We would like to emphasize that fact that inversion results are subject to a degree of uncertainty. Therefore, the net-to-gross gridblock properties are also conditioned to production data by slightly adjusting their values in Phase 1 and Phase 2 deposits separately through net-to-gross multipliers (adjustment factors) in Stage I.

It is important to distinguish between the net-to-gross gridblock property in the geological model and the meander belt net-to-gross parameter admitted by the proxy function that computes the connectivity factors. Meander belt net-to-gross parameter is a stratigraphic architecture characteristic of the meander belts in the subsurface. More precisely, it is the ratio of the rock volume of sands (typically dominated by the volume of sands in channel storeys) to the total volume within a meander belt. In other words, it represents the ‘geological net-to-gross’. On the other hand, the gridblock property net-to-gross is a product of the seismic inversion. It is a fraction that adjusts the gridblock volume to compensate for the existence of subgrid-scale shales. Presence of subgrid-scale shales is due to the fact that the geological model is at L3- and L4-level coarseness (Fig. 2). If the model were constructed at L0-level detail, then there would be no need for the net-to-gross gridblock property. The gridblock property net-to-gross and the meander belt (geological) net-to-gross are, by definition, different but related properties. Therefore, during the Stage I history-matching process, special attention was given to ensuring that the ranges of meander belt net-to-gross distributions, which enter into the connectivity factor calculation, span the corresponding gridblock property net-to-gross ranges that can be created during the history-matching process via net-to-gross adjustment factors.

### Upscaling

Static to dynamic model transfer involves some degree of upscaling to reduce the computational effort spent on multi-phase flow and transport calculations, in turn, the overall run time of a simulation. Upscaling involves simplification of the properties and cellular dimensions in the model. It is critical to accurately preserve the connectivity of reservoir sands across meander belts during the transfer from static to dynamic model.

The static model is built using a surface-conforming cornerpoint gridding technique (Fig. 6b) to honour stratigraphic surfaces interpreted from seismic data (Fig. 6a). The geocellular model is of size 161 × 157 × 32 in the *x*, *y* and *z* directions, respectively. The model is constructed at a lateral cell size to be directly used in the simulation model. Lateral upscaling (*x* and *y* directions) to increase cell size can damage meander belt and shale drape geometries and reservoir juxtapositions across them. This is because it causes changes to the meander belt and shale drape locations, object widths and lengths, and modifies the areas of sand-on-sand juxtaposition surfaces. Thus, the areal grid size of the static model (25 × 25 m) is retained in the dynamic model. The model is upscaled in the vertical direction to 161 × 157 × 21 with a flow-based technique (Renard & de Marsily 1997). Vertical upscaling is performed only across meander belt infill layers (Phase 2). Infill layers are separately upgridded for each meander belt. In other words, upscaling across meander belts is carefully avoided to keep intact the key stratigraphic surfaces. The upscaling workflow preserves the geometry and property heterogeneity of the static geocellular model. As will be shown later, retaining plumbing/flow characteristics of the static model while reducing the geological complexity to its salient chacteristics greatly simplifies the production data-integration effort. Average run time for the geological pre-processing (shale drape and fault seal modelling) plus flow-based upscaling plus simulation workflow in the history-match mode is 36 min. In the forecast mode, an average 20-year simulation takes 97 min. on a 3.60 GHz Intel Xeon microprocessor.

The dynamic model is quality-checked qualitatively by use of 3D views and 2D cross-sections to ensure that key geo-objects are in place and the inter-channel connectivity characteristics observed in the static model are replicated in the dynamic model. The model is also quality-checked quantitatively: (a) upscaled and un-upscaled models yield dynamic results within 1% of each other; (b) pseudo-logs of porosity and net-sand computed at well locations compare reasonably well to profiles interpreted from gamma-ray logs.

### Model uncertainties

Parameters of the dynamic model that are modified by the data-integration algorithm to make the model honour observed production data are referred to as model uncertainties or simply as model parameters. Hereafter, we will use these terms interchangeably. Model uncertainties are classified into the following categories for the Stage I data-integration effort (NA): (i) shale drape coverage values for Phase 2 meander belts, (ii) global horizontal permeability ratios for Phase 1 and Phase 2 sands, (iii) regional (localized) horizontal permeability adjustment factors for Phase 1, (iv) Phase 1 and Phase 2 *k*_{v}/*k*_{h} (permeability anisotropy) ratios, (v) Phase 1 and Phase 2 net-to-gross adjustment factors, (vi) well skin factors, (vii) Phase 1 and Phase 2 connectivity factors and connectivity factor multipliers and (viii) the seal factor (transmissibility multiplier) representing the impact of drape coverage at the Phase 1–Phase 2 interface, controlling the level of communication between Phase 1 and Phase 2 sands. In Stage I, the initial model parameter ensembles are generated using uniform distributions that enforce the parameter ranges.

Stage II data-integration effort (EnKF) adjusts gridblock permeabilities while the values of model uncertainties determined subsequent to Stage I remain unchanged. It is worthwhile noting that (1) to initialize the EnKF, an ensemble of permeability realizations are necessary; (2) Stage I assimilates permeability information into the model (from production data) in terms of permeability adjustment regions (a more global approach). In order to make good use of this rather less-localized, yet, still valuable information and to start the Stage II work from a more intelligent point, the initial ensemble of Stage II is created such that the mean permeabilities of the regions honour the ranges obtained subsequent to the Stage I effort.

### Data-integration results

The model parameters of Stage I are reported in Table 1, together with P10 and P90 values of their prior distributions, respectively referred to as Prior Min and Prior Max. In Stage I, a number of NA runs were carried out with various algorithmic parameter selections governing the explorative vs. exploitative nature of the sampling (Sambridge 1999). Posterior inference is computed from model ensembles generated during the sampling process for each NA run, separately. Although posterior distributions are computed via post-processing, best history-match models yielded by the NA runs are set aside for further analysis (described below) in Stage I for the reason of practicality. This is mainly because the best-match model parameters are reasonably close to most likely parameter combinations in posterior distributions. Four best-match models found by NA are documented in Table 1. These models are referred to as HM#1, HM#2, HM#3 and HM#4. HM#1 and HM#4 yield similar meander belt shale drape coverage values around 0.65. HM#3 exhibits major differences from the remaining models, with high meander belt shale drape coverage (∼0.80), large Phase 2 *CF*_{th} (0.85 compared to the remaining models where Phase 2 *CF*_{th} is approximately 0.70), and larger injection well skin factors. Phase 1 *CF*_{th} is consistently low and lies in a relatively tight range between 0.33 and 0.38.

Well data-match results for HM#1–HM#4 models are documented in the panels of Figures 7 and 8 for injection and production wells, respectively. Part of the production data is intentionally left out from the data integration. This record is used to test how well history-match models predict future dynamic behaviour. Results of this test are also shown in Figures 7 and 8. Prior and posterior distributions of four (example) model parameters are displayed in the panels of Figure 9. They illustrate the exploration quality of NA sampling during one of the runs that led to HM#1. It is important to emphasize that, among the model parameters shown in Figure 9, Phase 2 meander belt shale drape coverage and Phase 2 trapped-hydrocarbon connectivity factor are the ones whose uncertainty was significantly reduced by integrating production data.

Various manifestations of the objective function over the course of NA iterations are shown in Figure 10 for an example run, the one that yields HM#1. This figure displays the progress of the overall minimum misfit, minimum misfit attained at a given iteration and the average misfit at a given iteration. The latter is a relevant measure of overall progress because NA samples the misfit function for an ensemble of models at a given iteration.

HM#1 model parameters are set constant in Stage II. Two alternative gridblock permeability initial ensembles are generated. Each ensemble is of size 100. The main difference between these two initial ensembles is the variogram model used while generating the prior (logarithmic) permeability fields via two-point statistics (Deutsch & Journel 1998). Two scenarios were investigated: a long correlation-length scenario and a relatively shorter correlation-length scenario. In these initial realizations, region-by-region mean permeability values honour the values obtained post-Stage I. The resulting two initial ensembles are separately processed through EnKF, giving rise to models HM#5 and HM#6 in Figures 7 and 8. Models HM#5 and HM#6 represent post-EnKF (Stage II) ensemble means. Pre- and post-Stage II horizontal permeability fields (*k _{x}* =

*k*) are illustrated in Figure 11. By design, EnKF yields a more smoothly varying permeability field compared to NA. This study shows that NA performs the bulk of the production data-integration task using high-impact model uncertainties (Stage I). EnKF fine-tunes the dynamic models with more localized smooth updates to gridblock parameters (Stage II). EnKF renders the permeability field geologically more realistic by smoothing the sharp interfaces using variogram information. These interfaces arise from the parametrization restrictions imposed on the model for making it more amenable to NA. HM#6 yields a better quality match compared to HM#5 among models produced by EnKF (Figs 7, 8). Breakthrough-time prediction in Producer #2 is slightly more accurate in HM#6 compared to other history-match models. However, the late-time water rate behaviour (during the history-match period) is better represented by HM#1 (Fig. 12). Overall, EnKF (Stage II) does not improve the data match significantly over NA (Stage I). In fact, the quality of the BHP match in Injector #2 decreases when EnKF is applied subsequent to NA (HM#6). Having stated that, EnKF delivers a notably more geologically consistent permeability field (Fig. 11). One can argue that the NA – EnKF (Stage I– Stage II) sequence has to be iterated a few times to obtain further tuned history matches. We did not follow this route of action because history matches produced by a single (two-stage) pass were found satisfactory for the purpose of reservoir forecasting.

_{y}The effect of ignoring the subseismic channel architecture on the behaviour of the HM#6 is explored by setting connectivity factors to 1.0. Historical and history-matched water production rate responses with and without connectivity factors in Producer #2 are shown in Figure 12. The time of water breakthrough is delayed in the absence of connectivity factors, in line with the rationale of the 3D connectivity approach. As shown in Figure 12, the 3D connectivity factors have a strong imprint on the recovery behaviour. Although not demonstrated by Figure 12, in the absence of connectivity factors, it was not possible to obtain an acceptable match for the water rate behaviour in Producer #2. This was the case despite adjusting the remaining model parameters (within their physical bounds) via the Stage I algorithm. One can argue that physical bounds on the model parameters (other than connectivity factors) could have been revisited rather than introducing connectivity factors into the data-integration process. This is not a preferable approach, because bounds on model parameters reflect our prior knowledge and are tied to hard data. It is not justifiable and rigorous to permit the data-integration process to obliterate these bounds to obtain ‘just a match’ with a possible good fit to the measured data. The forecasts made with such a history-matched model would be of questionable predictive value. On the other hand, it was quantitatively demonstrated in Alpak *et al*. (2010) that (1) subseismic channel architecture matters in channelized turbidite reservoirs, and (2) it is possible to represent its effect by use of pseudo-relative permeability functions derived from connectivity factors. Therefore, it is more appropriate to include the connectivity factors in the data-integration process. Lower and upper bounds are imposed on the history-matched connectivity factors (as in the case of the remaining parameters) based on *a priori* knowledge stemming from cores, well logs, outcrop observations and high-resolution seismic data acquired in shallow-water analogues. It is worth pointing out that the history-matching study was carried out within two weeks.

### Forecast study

HM#6 model is used to test the sensitivity of the recovery factor behaviour for two alternative scenarios. The first scenario operates the field with the current well configuration and is appropriately called the ‘no further activity’ (NFA) scenario. The second scenario (IP) evaluates the updip additional recovery potential. This scenario requires a new injector–producer pair in the updip sector of the reservoir.

IP scenario wells are placed using an automatic well location optimization algorithm (Castineira *et al*. 2009). A movable oil map approach is used to derive desirable initial-guess locations that are as close as possible to the optimal locations. A movable oil map is the projection of the spatial distribution of mobile oil volume on the 2D horizontal view of the reservoir. The volume of movable oil in the reservoir varies with time. The movable oil script programmed into the simulation deck generates maps at desired periodic intervals during simulation. At a given horizontal location in the map, a movable oil map conveys the vertically integrated (stacked) mobile oil volumetric information through colours. Movable oil calculation inherently takes into account the gridblock volume, net-to-gross, porosity, oil saturation and residual oil saturation. Use of movable oil maps constitutes a practical technique for analysing reservoir sweep and connectivity. Time-dependent analysis of movable oil maps has been found useful in assessing the dynamic impact of alternative placement strategies for the updip injector–producer pair.

Initial, post-HM and post-forecast movable oil maps of Reservoir A are shown in Figure 13, together with well locations for each scenario. The NFA scenario leaves behind a substantial amount of unswept hydrocarbons, whereas the IP scenario recovers most of the updip hydrocarbons. Estimated ultimate (water injection) recovery factor is around 0.16 for the NFA scenario, while the one for the IP scenario is approximately 0.31 (Fig. 14).

For the IP scenario, the uncertainty of the recovery factor (as a function of injected pore volume) with respect to subseismic architecture is evaluated. This is accomplished by sampling the HM#1 *a posteriori* connectivity factor cumulative density functions (estimated by the Stage I algorithm) at regular percentile intervals between P10 and P90. Here, it is important to remember that HM#6 is obtained by processing the HM#1 model through EnKF (Stage II algorithm) and the connectivity distributions obtained subsequent to Stage I algorithm remain unchanged. The estimated variability is around 0.06–0.08 in terms of recovery factor units. Another interesting result is that the recovery factor function computed using the P50 connectivity factor values (estimated from the post-Stage I ensemble) is very similar to the one computed by employing the deterministic best-match connectivity factor values (Fig. 14).

For the IP scenario, other history-matched models (HM#1–HM#5) result in updip injection and production well locations in the vicinity of the ones found for HM#6. The recovery factor ranges between 0.15 and 0.20 for the NFA scenario among all history-matched models (HM#1–HM#6). The recovery factor range computed by taking into account the connectivity factor uncertainty (as shown in Fig. 14 for HM#6) is from 0.27 to 0.35 for all history-matched models.

The forecasting study described above demonstrates how connectivity factors can be utilized to rapidly evaluate the impact of subseismic architectural parameters on the recovery behaviour. Similar studies can be carried out for other model parameters. The aforementioned reservoir forecasting study was completed within one week.

## CONCLUSIONS

A dynamic modelling workflow is developed that integrates geologically based pseudo-relative permeabilities into a two-stage automatic history-matching algorithm. Large-scale reservoir parameters, such as shale drape coverage, reservoir connectivity, net-to-gross, and the amount of communication across the valley separating major reservoir compartments, are history-matched in the first stage using a novel implementation of the Neighbourhood Algorithm (NA). Ensemble Kalman Filter (EnKF) operates in the second stage to introduce smooth local changes to gridblock permeabilities. Large-scale parameters identified in the first stage remain unchanged in the second stage.

The data-integration and forecasting workflow is applied to a channelized turbidite reservoir in West Africa. We show that an ensemble of history-match models can be produced in a matter of two weeks and a simulation-based reservoir forecasting study can be completed within one week, instead of months. It is also shown that history-matched models generated with the two-stage data-integration workflow predict future recovery with reasonable accuracy. Use of 3D connectivity effective properties accelerates the history-matching process. Coarse dynamic models reduce the simulation cycle-time, enabling the history-matching procedure to operate on an ensemble of alternative subsurface realizations. A sufficiently diverse ensemble of alternative history-matched dynamic models encompasses a greater potential for accurately capturing forecast uncertainties. The connectivity effective property introduces physically based subgrid-scale prior information into the history-matching practice. The added information accounts for the effects of missing scale and, therefore, simplifies the history-matching task. It is of crucial importance to link the connectivity effective properties to geological hard data and enforce bounds on their variability in the history-matching mode. In the absence of a strong geological basis, use of connectivity factors in conjunction with coarse-scale realizations can lead to dynamic models that honour historical production measurements but lack predictive capability. From the history-matching viewpoint, 3D connectivity factors dominate the water rate response in wells Producer #1 and Producer #2. Having stated that, 3D connectivity factors by themselves are insufficient to obtain a good history match. Additional history-matched model parameters (in Stage I) include (a) shale drape coverage value for Phase 2 meander belts, (b) global horizontal permeability multipliers for Phase 1 and Phase 2 sands, (c) regional (localized) horizontal permeability adjustment factors for Phase 1, (d) Phase 1 and Phase 2 *k*_{v}/*k*_{h} (permeability anisotropy) ratios, (e) Phase 1 and Phase 2 net-to-gross adjustment factors, (f) well skin factors and (g) the seal factor (transmissibility multiplier) representing the impact of drape coverage at the Phase 1–Phase 2 interface (controlling the level of communication between Phase 1 and Phase 2 sands). Moreover, the data-integration algorithm modifies the permeability field in Stage II.

Experience indicates that for a few highly influential large-scale parameters, a global inference technique (e.g. NA) is more suitable since it allows the likelihood function to climb hills in the uncertainty space. However, for gridblock properties a relatively smoother update technique, such as EnKF, is more appropriate. There are two important observations: (1) NA carries out the bulk of the production data-integration task in a number of investigated cases using a few key high-impact uncertainties (Stage I). EnKF contributes to the data-integration process by fine-tuning the dynamic models with relatively more localized smooth updates to gridblock parameters (Stage II). (2) It is, therefore, possible to decouple one complex history-matching problem into two relatively simpler stages and to tackle these two stages with separate techniques.

## Acknowledgments

The authors would like to thank Frans van der Vlugt, Ciaran O'Byrne, Steve Naruk, Kachi Onyeagoro and Fa Dwan for their contributions to the static modelling work. The authors are grateful to Shell International Exploration and Production Inc for permission to publish this paper.

- Received November 9, 2009.
- Accepted July 12, 2010.

- © 2011 EAGE/Geological Society of London