Integrating Data of Different Types and Different Supports into Reservoir Models

Integrating Data of Different Types and Different Supports into Reservoir Models — In this paper, we focus on the joint integration of production and 4-D inverted seismic data into reservoir models. These data correspond to different types and different scales. Therefore, we developed two-scale simulation workflows making it possible to incorporate data at the right scale. This issue also emphasized the need for adapting traditional history-matching methodologies. For instance, the formulation of the objective function and the development of customized parameterization techniques turned out to be two key factors controlling the efficiency of the matching process. Two application examples are presented. The first one is a small-size synthetic field case. It aims to build a set of reservoir models respecting either production data only or both production and 4-D seismic-related data. It is shown that the incorporation of 4-D seismic-related data in addition to production data into reservoir models contributes to reduce the uncertainty in production forecasts. The second example is a field in the North Sea offshore Norway operated by Statoil. It stresses difficulties in conditioning reservoir models to both real production and 4-D inverted seismic data among the very large number of uncertain parameters to handle and the comparison of real noisy data with numerical responses. Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 67 (2012), No. 5, pp. 823-839 Copyright © 2012, IFP Energies nouvelles DOI: 10.2516/ogst/2012024 Pore2Field – Flows and Mechanics in Natural Porous Media from Pore to Field Scale Pore2Field – Physique des écoulements en milieux poreux naturels : de l'échelle du pore à l'échelle du réservoir IFP Energies nouvelles International Conference Rencontres Scientifiques d’IFP Energies nouvelles ogst110203LeRavalec 3/12/12 16:06 Page 823


INTRODUCTION
Reservoir heterogeneity exists at multiple scales and is related to complex and intricate geological formation processes such as bedform migration, erosion and deposition. Thus, channels, lobes, bars or meanders can be observed at large scale. Heterogeneity can be also evidenced inside these geological objects. For instance, lateral accretion and meter-scale crossbedding are often observed in channels. Besides, changes in the nature of contacts, grain size or fracture density are commonly detected at the much smaller scale of cores. Reservoir heterogeneity results in variations in rock physical properties. For instance, permeability of the Culebra Dolomite aquifer, approximately 450 m above the repository horizon of the Waste Isolation Pilot Plant in southeastern New Mexico, USA, was shown to vary by five orders of magnitude (Holt, 1997). On the other hand, Henriette et al. (1989) made detailed permeability measurements of porous blocks of sandstone and limestone measuring 15 × 15 × 50 cm 3 . Each block was cut into three hundred small plugs and permeability was measured for each of them. Permeability variations were shown to approximately span two orders of magnitude.
The ever-growing processing power of computers makes it possible to build numerical geological models which are used to predict reservoir performance, optimize production planning, make key economic decisions, etc. In a nutshell, a model is a numerical representation of the spatial variations in petrophysical properties inside a geological formation. It consists of a three-dimensional grid populated with porosities, permeabilities, initial oil saturations, etc. This model can be input into a flow simulator to evaluate fluid displacements: the reliability of the simulated responses depends on the reliability of the model itself. However, no model can be shown to accurately represent the processes responsible for the past and future observed fluid flows. Clearly, a model incompatible with the data collected in the field is wrong. The best we can hope for is that the more data-consistent the model, the sounder the predictions.
There are many sources of data. Measurements on core samples extracted from wells penetrating the geological formation provide very high-resolution information. Downhole logging tools are used to describe petrophysical variations along wells. However, such direct measurements of petrophysical properties are very sparse and sample only a small reservoir volume. They have to be supplemented by indirect measurements. Due to its low resolution, 3-D seismic is routinely used for defining only the structural model. However, unlike laboratory and log data, it provides information over large areas, which makes it an invaluable candidate for better characterizing petrophysical properties. Instead of focusing on raw seismic data, we can use inverted seismic data, which basically comprise velocities and impedances. In this case, seismic inversion refers to the inverse modeling of reservoir properties from raw seismic data. Inverted seismic data can help to capture reservoir property variability between wells. All of these data are considered as static because they do not change with time. Other data are collected during reservoir producing life. Because they depend on fluid flows, they are said dynamic. They mainly consist of production data, i.e., data measured at wells such as bottom hole pressures, oil production rates, gas-oil ratios, tracer concentrations, etc. Since the late nineties, they also include data derived from 3-D seismic surveys repeated at successive times, that is 4-D seismic data (Eastwood et al., 1994;Benson and Davis, 2000;Arts et al., 2002;Behrens et al., 2002;Guderian et al., 2003;Roggero et al., 2012). Ideally, a base 3-D seismic survey is acquired before starting production. Then, monitor surveys are successively acquired after a few years of production. Differences between seismic responses are caused by fluid movements, pressure changes, temperature changes, fluid property or compositional changes, or rock changes. They help to monitor fluid fronts and pressure domains between and beyond wells.
Each of the data types mentioned above is associated to its own scale of measurement, its own level of precision and a given method of measurement with its own physical principles. The major challenge in today's reservoir characterization is to integrate all different kinds of data into reservoir models for the purpose of evaluating large-scale reservoir performance.
Due to the lack of direct measurements and the complexity of geological depositional processes, a deterministic description of the spatial variations in reservoir properties is unachievable. Referring to geostatistics is a usual practice since the pioneering work of Matheron (1963Matheron ( , 1965 in mining engineering. Within this framework, petrophysical properties are considered as random variables and their spatial variability is inferred from the available static data. Thus, realizations of these variables are randomly drawn to populate the geological model. Such realizations are usually constrained to the observed static data referring to kriging techniques (Delhomme, 1978;Chilès and Delfiner, 1999). Besides, the integration of dynamic data is an inverse problem solved using optimization techniques. It is also known as history-matching since it originally focused on the integration of production history. Mathematically, history-matching relies on the minimization of an objective function. Given a geological model, the objective function quantifies the least-square mismatch between the actual dynamic data and the corresponding responses simulated for the geological model of interest. The optimization process drives the successive adjustments of this model until the data mismatch is small enough. The final so-constrained model respects both static and dynamic data. A difficulty in such inverse problems is underdetermination: we are not guaranteed to have enough information to uniquely determine the solution. It is thus essential to integrate as many data as possible into reservoir models.
The multiscale character of reservoir properties and data makes the problem of predicting reservoir performance a natural target for multiscale methods. Geological models can be considered as multiscale models spanning several orders of magnitudes from core scale to reservoir scale. The current practice in reservoir engineering for handling various scales is based upon upscaling. The fine geological model is upscaled to the scale where data have to be incorporated. This calls for accurate upscaling methods able to transfer physical processes and properties from one scale to another. A promising alternative is multiscale flow simulation, which can provide accurate resolution of both large-scale and fine-scale flow patterns. Basically, the leading idea in multiscale modeling is to take advantage of both the simplicity and efficiency of large-scale models and the accuracy of fine-scale models. What set these techniques apart is that fundamental physics principles can be rigorously applied at any scale. Then, the computed properties are passed to the next scale up. Multiscale modeling is CPU-time intensive, which limited its extended use up to now. However, the continuous advances in computational technology can give the possibility to envision to go one step ahead. This motivated recent works in flow simulation. Gautier et al. (1999) developed a nested approach to solve the pressure equation at the fine scale. Then, they integrated the saturation equations along streamlines. Referring more specifically to multiscale modeling, Jenny et al. (2006) developed a multiscale finite volume method for multiphase flow. They used a sequential scheme for dealing separately and differently with pressure and saturation. They pointed out the accuracy and computational efficiency of the method compared to standard finite volume solutions of the fine-scale problem. A multiscale finite element approach was also presented by Efendiev et al. (2006) for two-phase flow simulation. The need for nested approaches also emerged in geostatistical modeling. Weissmann and Fogg (1999) suggested to refer to sequence stratigraphic concepts to account for various scales. Reservoirs are thus divided into strata and strata into units that have similar physical properties. Journel et al. (1998) proposed to follow a nested approach with first the simulation of facies, then of porosity and permeability values populating facies. On the other hand, kriging techniques were developed for addressing scaling issues: they incorporate different data and produce estimates representative of averages over given areas (Roth et al., 1998;Chilès and Delfiner, 1999;Tran et al., 1999;Aanonsen and Eydinov, 2006). The development of multiscale techniques, although mainly restricted to flow simulation today, should favor the emergence of adaptive history-matching methodologies. Many challenging issues will have to be addressed. For instance, there will be a need for parameterization techniques able to update reservoir heterogeneity in a sequential manner, from the coarse to the fine scale. A key question is about the identification of the finest scale to be accounted for. Shall we envision multiscale flow simulation schemes with scales going from reservoir to core to better describe multiphase flow, hence to better predict reservoir performance? The finer we will go, the more expensive the flow simulations.
In this paper, we come to a balance point so that simulation is feasible and relevant information is kept. In particular, we focus on the joint integration of production and 4-D inverted seismic data into reservoir models while keeping consistency with static data. We design specific simulation workflows capable of simultaneously handling two scales: the geological and the reservoir scales. This provides a basis for building reservoir models respecting the required constraints. However, accounting for distinct scales is not enough. The efficiency of the proposed solution strongly depends on the methods used to measure the data mismatch and to adjust the initial model. Two application examples are presented in the last sections. The first one is a small-size synthetic field case inspired from the known PUNQ-S3 benchmark case. It aims to emphasize the added value of 4-D seismic for improving the reliability of reservoir models. The second example deals with the building of a reservoir model representing an off-shore North Sea field. In this case, the purpose is to stress some of the issues to be solved when studying a real field. This example focuses on the joint integration of production data and 4-D seismic impedances.

SIMULATION WORKFLOWS
The overall problem consists in determining a model providing numerical production and seismic responses similar to the data collected in the field when given as input into a simulation workflow. This simulation workflow includes the following activities: geological modeling, upscaling, flow modeling and Petro-Elastic Model (PEM) (Fig. 1  The first activity, that is geological modeling, aims to populate the grid with petrophysical properties using geostatistical simulation algorithms (Journel et al., 1998;Chilès and Delfiner, 1999). First, we generate facies realizations, which respect the observations at wells and the facies proportions derived from logs and seismic. Porosity realizations are then simulated conditionally to the available porosity data to populate the facies realizations. At this point, permeability can be computed from empirical relationships between porosity and permeability measurements collected from cores or logs. It can be also simulated conditionally to the available permeability data and the previously simulated porosity realization. The resulting model is called a geological model: it corresponds to a fine scale and usually encompasses up to a few dozens of millions of grid blocks.
The following step entails the upscaling of the geological model to reduce the number of grid blocks and make the subsequent flow simulations feasible in a reasonable amount of time. Given a fine grid, the petrophysical properties populating the fine grid and a coarse grid, upscaling boils down to the computation of the equivalent petrophysical properties to be assigned to the coarse grid blocks. This process results in a coarse model, called reservoir model, which has to reproduce as closely as possible the flow behavior of the fine geological model. This reservoir model, which is associated to a coarse scale, generally contains between 100 000 and a few millions of grid blocks.
Then, flow equations are solved to simulate fluid displacements within the reservoir model (Aziz and Settari, 1979). The computer programs, which solve flow equations are called fluid flow simulators. The simulation outputs are production responses such as bottom hole pressures, flow rates or water cuts at wells. They also include pressure, saturation and temperature grids describing the variations in these variables over the whole reservoir model at base and monitor times. The fluid flow simulator only provides approximate solutions: its accuracy depends on grid block size. However, the finer the grid blocks, the more time demanding the flow simulator.
Since seismic responses (velocities, impedances) depend on petrophysical properties, they are computed at the fine scale. However, they also depend on the pressure, saturation and temperature grids provided by the flow simulator at the coarse scale. These grids are first moved to the fine scale throughout a downscaling step. Second, a petro-elastic model is applied to convert fluid properties (i.e., pressures, saturations, temperatures, fluid densities, fluid elastic properties) and rock properties (i.e., porosity, matrix elastic properties) into simulated elastic responses, which are used to build the required seismic responses. Very often, the petroelastic model mixes theoretical equations (Gassmann, 1951;Mindlin, 1949) and empirical relationships calibrated from core and log data (Mavko et al., 1998). Then, these numerical seismic responses are transferred into the seismic bandwidth using an appropriate band pass filter.
The interested reader can refer to Le  for a detailed description of this simulation workflow.

INVERSE PROBLEM
The determination of a reservoir model consistent with the measured production and seismic data is an inverse problem. History-matching or model calibration are essentially the same as inverse modeling. This problem can be solved with optimization, meaning that an objective function, usually defined in a least-square sense, has to be minimized: (1) J is the objective function. It depends on parameters m. Vector d obs includes the data to be matched. Covariance matrix C D quantifies uncertainties in data and modeling: it is assumed to be diagonal although this hypothesis may be questionable for seismic attributes. Operator g represents the simulation workflow, which maps the parameter space to the data space. A review of the inverse problem can be found in Carrera et al. (2005) or Oliver and Chen (2011).

Objective Function
When considering both production and inverted seismic data, the objective function is written as a two-term function: one for the production data mismatch and the other one for the seismic data mismatch. The effectiveness of the joint inversion of production and seismic data hinges on the ability to properly capture differences between reference data and simulated responses. The least-square formulation is known to be appropriate for production data but production data and seismic data are very different. First of all, the support of seismic data is the whole grid while production data are localized at wells. Consequently, we may have millions of seismic data to be balanced with hundreds of production data. In addition, the seismic data considered are not raw data: they are derived from a preliminary inversion of lowresolution seismic measurements. Thus, there is more uncertainty in inverted seismic data than in production data. Several research papers showed that the least-square formulation was not relevant for quantifying the seismic mismatch (Aanonsen et al., 2003;Roggero et al., 2012). Tillier et al. (2012) suggested a different metric called Local Dissimilarity Map (LDM). The leading idea is to evidence meaningful features instead of exactly matching seismic data in every grid blocks. The approach involves three steps. The seismic data grid and the corresponding seismic response grid are first submitted to filtering and classification to enhance the important features. Second, the local dissimilarities between the two grids are derived from a local modified Hausdorff distance. For simplicity, we assume that the studied seismic attribute is split into two classes so that the data and simulated response grids are converted into binary black and white grids. The detected seismic features are identified as black. We apply the local modified Hausdorff distance and obtain a dissimilarity grid. Let us consider a given grid block. Its local dissimilarity distance is zero when its two values in the data and simulated response binary grids are identical. Otherwise, it is set to the (non zero) distance to the nearest black grid block. The last step calls for the computation of a global dissimilarity, which is the sum of the squared local dissimilarities over the whole grid. The so-obtained value is incorporated into the objective function as the seismic mismatch term. Tillier et al. (2012) proved the efficiency of this approach for conditioning a synthetic SAGD operated field case.

Parameterization
When Equation (1) is used as the objective function, the inverse problem is usually ill-posed. This implies that there is no unique solution and that the identified solution is highly sensitive to variations in the production and seismic data. To make the problem "more well-posed", we refer to regularization (Neuman, 1973;Carrera et al., 2005). A number of methods have been proposed in the literature. A possibility consists in narrowing the search space. This can be achieved by selecting an appropriate parameterization technique (Oliver and Chen, 2011) incorporating some knowledge about the properties of the searched solution. An example of restriction consists in setting upper and lower bounds to limit parameter variations. The choice of the parameterization technique is all the more decisive because it makes it possible to decrease or not the objective function. Focusing on the geological scale, a few authors developed geostatistics-based techniques to vary petrophysical properties from a limited number of parameters while preserving spatial variability. This assumes some kind of geological continuity. De Marsily (1978) and RamaRao et al. (1995) introduced the pilot point method, Hu (2000) the gradual deformation method, Caers (2003) the probability perturbation method and Le Ravalec-Dupin and Da Veiga (2011) the co-simulation perturbation method. These methods aim to perturb the realizations populating the geological model with petrophysical properties. They are used to change the spatial distribution of heterogeneities while respecting prior geological knowledge. In particular, the basic principle of the gradual deformation method is that the sum of two Gaussian random functions is a Gaussian random function. In its simplest form, it involves the linear combination of two independent Gaussian white noises z 1 and z 2 as: (2) t is a deformation parameter. Whatever its value, z is a Gaussian white noise too. On the other hand, the generation of realizations of a random function depends on a set of z t z t z t ( ) = ( ) + ( ) 1 2 cos sin π π random numbers like uniform or normal deviates. Therefore, the gradual deformation scheme makes it possible to vary the random numbers used to generate a realization. This boils down to the gradual deformation of the realization itself. Gradual deformation applies to continuous or discrete realizations as shown in Figure 2 and Figure 3, respectively. Such a deformation feature is of interest for minimizing the objective function. It makes it possible to investigate the search space and to pick up a realization, hence a geological model, leading to numerical responses closer to the target production and seismic-related data. However, the Gaussian white noise chain built from Equation (2) corresponds to a very tiny part of the search space. The probability of getting a good match by investigating this chain solely is very small. Thus, the search process must be iterated. In practice, the overall procedure is as follows (Fig. 4). A first chain is created from the combination of the initial Gaussian white noise (z 1 ) with another one (z 2 ) randomly generated. Then, a first optimization process yields an optimal deformation parameter, that is an optimal geological model for which the simulated dynamic response is as close as possible to the actual data. This optimal deformation parameter is used for updating z 1 while another random draw provides a new z 2 . Again, applying Equation (2) yields a new chain, which can be investigated again from a new optimization process and so on, until getting a reasonable data match. The parameterization methods listed above were shown to be efficient for calibrating reservoir models to production data (pilot point, gradual deformation, probability perturbation, co-simulation perturbation). However, they do not provide enough flexibility when considering also inverted seismic data (Roggero et al., 2012). One way of improving the matching of seismic data is the careful identification of the sub-domains to be adjusted. Early works considered Voronoi diagrams defined from wells (Hoffman and Caers, 2005;Gervais et al., 2007). Then, physically rooted subdomains were defined from drainage areas (Le Ravalec-Dupin and Fenwick, 2002;Hoffman and Caers, 2005).  determine influence zones by tracing the streamlines arriving to areas with significant 4-D seismic data mismatches. Da Veiga and Gervais (2011) apply filtering and classification techniques on 3-D seismic data error maps. In addition, the parameterization techniques themselves must be tailored to account for seismic data. This has steered the development of tools like the pilot block method (Le Ravalec-Dupin, 2010) able to significantly and locally vary realizations. This parameterization method reconciles a current practice in history-matching, that is the use of multipliers, with the theoretical framework of the pilot point method. It gives the possibility to change the mean of a realization over a target sub-domain while preserving spatial continuity.
Whatever the parameterization technique considered for varying petrophysical properties at the fine scale, the resulting Gradual deformation of a continuous realization with a deformation parameter varying from -1 to 1.

Figure 3
Gradual deformation of a discrete realization with a deformation parameter varying from -1 to 1.
changes are propagated to the entire simulation workflow. Thus, a consistent link is maintained between the geological model at the fine scale and the reservoir model at the coarse scale all along the optimization process.

Optimization Workflow
In the two application cases studied in the following sections, we apply the general simulation workflow introduced in Section 1 and use the gradual deformation method to adjust the petrophysical properties populating the geological model. In other words, we first generate Gaussian white noises which yield realizations of the target petrophysical properties. Then, a flow simulation and a petro-elastic model are run to obtain numerical responses, which are compared to the available production and seismic-related data. This boils down to the computation of the objective function.
The optimization techniques usually considered to minimize the objective function can be split into two main groups. The first one, which is the most commonly referred to in history-matching calls for gradients (Tarantola, 1987;Sun, 1995). Derivatives can be derived from an adjoint-state method (Chavent et al., 1975) or more often from finitedifferences (Kiefer and Wolfowitz, 1952). In this case, parameters are successively perturbed and the objective function is evaluated for each of these perturbations. The second group of optimization techniques is independent of gradients. For instance, it includes the simplex method (Nelder and Mead, 1965), simulated annealing (Yudin, 1966), genetic algorithms (Goldberg and Kuo, 1987), etc. These techniques require a huge number of iterations, which cannot be managed when running CPU-time demanding flow simulations. This is the reason why we selected a finitedifferences based approach. The gradual deformation parameters are then sequentially adjusted until convergence. At this point, the optimal values provide new starting Gaussian white noises for a subsequent optimization process as explained in Section 2.2. The optimization processes are repeated until getting a satisfactory match.

APPLICATION CASES
Two application examples are described hereafter. The first one is a small-size synthetic field case inspired from the known PUNQ-S3 case often used for benchmarking historymatching methods. It aims to emphasize the added value of 4-D seismic. The second example deals with a real off-shore North Sea field case. It emphasizes how difficult it is in practice to condition reservoir models to both production data and 4-D seismic impedances.

Value of 4-D Seismic in Improving Forecasts
In this section, we focus on the PUNQ-S3 field case defined in the "Production forecasting with UNcertainty Quantification" project (Floris et al., 2001). We aim to build 10 constrained models following the methodology described above. Two cases are investigated. First, these 10 models are constrained to production data only. Second, they are constrained to both production data and 4-D seismic-related data. The purpose is to quantify the added value of 4-D seismic.
The PUNQ-S3 model is a small-size synthetic test case defined from a real field operated by Elf Exploration and Production. It consists of 5 independent layers with distinct petrophysical properties (Tab. 1). Layers 1, 3 and 5 have linear streaks of high-porous sands (porosity > 20%). These sand streaks of about 800 m wide are encased in a low porous shale matrix (porosity < 5%). Layer 2 includes marine or lagoonal shales with distal mouthbars or distal lagoonal delta. They translate into a low-porous (porosity < 5%), shaly sediment, with some irregular patches of somewhat higher porosity (porosity > 5%). Layer 4 is a lagoonal delta within lagoonal clay. It is associated to a low porosity matrix (porosity < 5%) with intermediate porosity patches (porosity ~15%). The spatial distributions of petrophysical properties are approximated by spherical variograms with parameters defined in Table 1. Whatever the layer, variograms are the same for porosity and permeability. In addition, the correlation coefficients between porosity and horizontal permeability and between porosity and vertical permeability are 0.8.
The formation is discretized into 19 × 28 × 5 grid blocks of which 1761 are active. Grid block dimensions are 180 × 180 m 2 along axes X and Y. Their thicknesses are varying. The field is bounded to the east and south by a fault and connects to the west and north to a strong aquifer. Due to the strength of the aquifer, injection is not required for pressure maintenance. A small gas-cap is located in the middle of the dome-shaped structure. The reservoir is produced by 6 wells (PRO-1, PRO-4, PRO-5, PRO-11, PRO-12 and PRO-15) located  around the gas-oil contact (Fig. 5). The production schedule spans a period of time between 01/01/1967 and 15/01/1975: it includes one year of extended well testing with successive production rates of 100, 200, 100 and 50 m 3 /day, three years of shut-in and then four years of field production with a rate of 100 m 3 /day. Each year, wells are shut in for two weeks for testing (Fig. 6). Wells are also submitted to a pressure constraint: when bottom hole pressure goes below 120 bar, this value is used as a target pressure and the oil flow rate target is canceled. The dynamic data collected during this first 8 year period are the pressures and oil flow rates at wells and 4-D seismic-related data. For simplicity, these seismicrelated data are the oil saturation grids simulated for January 1967 and January 1975. They correspond to the base and monitor seismic-related data. A reference reservoir model populated by porosity, horizontal and vertical permeability realizations was generated on the basis of the data reported in Table 1. The soobtained porosity values assigned to the first layer are shown in Figure 11a. Then, the production scheme described above was simulated for this reference reservoir model using the internal IFP Energies nouvelles flow simulator named PumaFlow TM . The resulting numerical responses give the reference dynamic data.
At this stage, the spatial distributions of porosity, horizontal and vertical permeabilities are assumed to be unknown. The data available are listed in Table 1. They also include the reference production and seismic-related data: pressures and oil flow rates at wells from January 1967 to January 1975 (from time 0 to 2 936 days) and oil saturations over the whole reservoir grid in January 1967 and January 1975. We apply the matching methodology introduced in Section 1 to determine 10 constrained models starting from 10 distinct initial reservoir models. The gradual deformation method is used to perturb the porosity and permeability realizations populating the 5 layers of the reservoir models. Three realizations are generated for each layer: one for porosity, one for horizontal permeability and one for vertical permeability.   The gradual deformation of each petrophysical realization is driven from two deformation parameters. Therefore, the total number of parameters to be adjusted for minimizing the objective function is 30. The minimization process is driven from a gradient-based algorithm with derivatives estimated from finite-differences (Kiefer and Wolfowitz, 1952). This significantly increases the number of flow simulations to be performed: parameters are successively perturbed and a flow simulation is run for each of these perturbations. The final optimal parameter values are used to build constrained models. Two cases are considered. First, the objective function includes the production data mismatch term solely. Second, the seismic-related data mismatch term is added to the objective function. Whatever the case considered, the constrained models reproduce very well the reference pressures and oil flow rates. Figure 7 and Figure 8 show the results obtained for the cumulative volume of oil produced. Less than 100 flow simulations are needed to decrease the production data mismatch to almost zero. Matching the seismic-related data is more challenging. In the examples studied, 400 flow simulations are typically required to reduce the seismic mismatch term by about 5 to 6 (Fig. 9).
For simplicity, we focus on porosity in Layer 1. Figure 10 shows the evolution of this property along the matching process for the fifth model studied. Figure 10a presents the porosity realization populating the initial model. Then, this initial model is modified to respect production data. The updated porosity realization is shown in Figure 10b. At this stage, the changes are negligible. The final matching phase consists in incorporating both production and seismic-related data. This now results in clear modifications of the porosity field (Fig. 10c). It is worth comparing the final constrained porosity model with the reference one (Fig. 11a). They are obviously not the same but they look alike. The means of the 10 constrained porosity models are displayed in Figure 11b,c, for the two phases of the matching process. Clearly, the models constrained to production data only are very different from the reference model (Fig. 11a). On the contrary, models constrained to both production and 4-D seismic-related data capture the main trends of the reference model. It can also be checked that uncertainty in porosity decreases when incorporating 4D-seismic-related data (Fig. 11d,e). Last, we investigate the evolution of the data mismatch relative to the various phases of the matching process. In particular, this is shown in Figure 12 for the fifth model. The first matching phase, concerned with production data only, induces a clear decrease in the production data mismatch and a moderate one for the seismic-related data. On the other hand, the second matching phase, which accounts for both production and seismic-related data, leads to an obvious improvement of the seismic match. At the same time, the production match slightly deteriorates but keeps excellent anyway. Cumulative volume of oil produced against time for the 10 models constrained to the reference pressures and oil flow rates collected at the 6 producing wells. The first 2 936 days correspond to the matching period while the following ones till 6 000 define the prediction period. The red dots on the left are reference data. Cumulative volume of oil produced against time for the 10 models constrained to the reference pressures and oil flow rates collected at the 6 producing wells and the oil saturation grids associated to times 0 and 2 936 days. The first 2 936 days correspond to the matching period while the following ones till 6 000 define the prediction period. The red dots on the left are reference data.
Last, we investigate the potential of the constrained models to forecast oil production. The forecasting period is defined from January 1975 to July 1983 right after the first 8 year matching period. During this second 8.5 year period, production conditions are unchanged. Two different types of forecasts are investigated. Oil production is evaluated first for the 10 models constrained to production data only and second for the 10 models constrained to both production and 4-D seismic-related data. Results are shown in Figure 7 and Figure 8. They stress the benefit of 4-D seismic data. The variability in the production forecasts computed for the models constrained to production data only is 23%. It is reduced to 13% when considering also 4-D seismic-related data. Seismic-related data mismatch term against the number of flow simulations. Light blue bars: perturbations performed for estimating gradients. Dark blue bars: iterations due to the minimization process. Pink bars (with labels indicating the number of flow simulations performed): beginning of the investigation of a realization chain following the gradual deformation scheme.

Figure 10
Porosity in layer 1 for the fifth model at different states. a) Initial state. b) State after matching production data. c) State after matching production and seismic-related data.

Conditioning of a North Sea Field Case to Production and 4-D Inverted Seismic Data
This section presents the matching results obtained for an off-shore North Sea field operated by Statoil. These ones were originally presented as SPE 135116 at the SPE Annual Technical Conference and Exhibition, 19-22 September 2010, Florence, Italy . The field is located in the Tarbert and Upper Jurassic formations. The Upper Jurassic deposits include highly turbiditic sandstones and shale while the Tarbert environment is made of regressive/transgressive delta-plains and shoreface systems. The production of the field started in 2000 and spread over 8 years. It had been driven first by depletion for four years. Then, water was injected from 2004 to maintain pressure. Production was ensured from seven production and two injection wells (Fig. 13). The structural limits of the reservoir correspond to faults (Fig. 13).
The purpose of this study was to determine a reservoir model constrained to the available production and 4-D a) Porosity in layer 1 -reference model; b) mean for the 10 models constrained to the reference pressures and oil flow rates; c) mean for the 10 models constrained to the reference pressures, oil flow rates and oil saturation grids. d) Variance for the 10 models constrained to the reference pressures and oil flow rates and e) to the reference pressures, oil flow rates and oil saturation grids.
inverted seismic data. Production data measured at wells consisted of water, oil and gas flow rates and RFT pressures and 4-D seismic data of acoustic impedances in 1998 and ratios of acoustic impedances in 2004 to acoustic impedances in 1998. RFT stands for Repeat Formation Tester. It is an open-hole wireline device capable of providing pressure data along wells at various depths before the beginning of production/injection. Impedance data being very noisy, there were first filtered so as to make the comparison with the corresponding numerically computed responses more meaningful. The resulting filtered seismic attributes turned out to be very similar within each MT2 unit. A typical feature of the impedance ratio data, specially in unit MT2_3A, is the strong decrease around and north of injector I1 (Fig. 14a). This was assumed to be related to the pressure build-up induced by water injection. The fine geological model (Tab. 2) encompasses 72 × 150 × 170 grid blocks of dimensions 50 m along both X and Y axes and variable thickness. It is split into several independent units. We focus hereafter on the simulation of petrophysical properties for the three Mid Tarbert (MT) 2 units, which are the main reservoir units: MT2_3B, MT2_3A and MT2_1&2 (the deepest one). The properties were provided by Statoil for the other units. A preliminary step was the random generation of facies using the pluriGaussian method and the proportion maps estimated from log data. Three facies were identified: Inner Estuary sands (IE), estuarine Heterolithics Sands (HS) and Transgressive Muds (TM). Then, porosity realizations were independently simulated to populate facies IE and HS. Their spatial variations were characterized by spherical variograms with parameters reported in Table 3. Porosity values were trimmed to fall within 0 and 0.33. Linear relationships between permeability logarithm and porosity were then applied to assign horizontal permeability values to all grid blocks. Again, these ones were trimmed to lie between 0 and 3 000 mD. Production data mismatch (purple) and seismic-related data mismatch (blue) for the fifth model at different states. Left: initial state. Middle: state after matching production data. Right: state after matching production and seismic-related data. The whole simulation workflow introduced in Section 1 was run to compute numerical production and impedances responses. As computation times had to be reduced, the fine geological model was upscaled into a coarse reservoir model with 36 × 75 × 71 grid blocks (Tab. 2). The horizontal dimensions of these coarse grid blocks were about 100 m × 100 m. Equivalent porosity and horizontal permeability values were derived from arithmetic averages to populate this coarse reservoir model. A horizontal/vertical permeability ratio of 0.1 was considered to estimate vertical permeabilities. Fluid flows were simulated from a black-oil model, meaning that there was no change in hydrocarbon composition as the field was produced.
The conditioning of an initial randomly drawn geological model was performed on the basis of 3 successive steps. First, the workflow introduced in Section 1 was applied to calibrate the initial model to production data. Second, a few additional parameters were manually adjusted to improve the 4-D inverted seismic data match, which unfortunately contributed to deteriorate the production data match. Third, the full matching workflow was run again to finally determine a model respecting both production and 4-D inverted seismic data.
Step 1. A preliminary sensitivity analysis was performed to identify a set of influential parameters among 45 before running any matching process. Then, 10 parameters were kept: the mean porosity of facies IE in unit MT2_3A, gradual deformation parameters used to vary the spatial distribution of facies IE and HS in units MT2_3B and MT2_12 and gradual deformation parameters used to vary the spatial distribution of facies IE porosity within the three MT2 units. The water flow rate at well P2 a) Vertical average value of the filtered measured; b) simulated impedance ratios in unit MT2_3A after matching production data; c) after matching impedances.
simulated for the so-obtained constrained model is displayed in Figure 15a,b. The data fit was better, in particular with a breakthrough time closer to the actual one. Nonetheless, no obvious improvement in the RFT pressure match could be achieved. The final value of the objective function is reported in Figure 16 (step 1production term). Besides, acoustic impedances in 1998 and 2004 were computed for the constrained model to be compared to the reference ones (Fig. 16, step 1 -seismic term). The resulting ratios were shown to strongly depart from data (Fig. 14b).
Step 2. Based on physical knowledge, a few parameters were then manually adjusted in order to reduce the seismic data mismatch. Basically, we tried to boost pressure around well I1 combining two strategies. Three faults (nwse07, nwse08 and main07, Fig. 13) located south of injector I1,  Water production rate against time at well P2. Red dots: data. a) Blue curve: response simulated for the initial model. b) Blue curve: response simulated for the model constrained to production data only. c) Yellow curve: response simulated for the model constrained to production data only; blue curve: response simulated for the model constrained to production and inverted seismic data.
were closed. In addition, the mean porosity in the unit MT2_3A was reduced from 0.21 to 0.15 in a local area around well I1. The porosity realization was then updated using block kriging. The impedance ratios computed for the resulting model are shown in Figure 14c. Although still different from the reference data, they were decreased around well I1 as requested. This induced a decrease in the seismic data mismatch term by about 3% (Fig. 16, step 2). However, this manual adjustment contributed to deteriorate the production data match: the production term in the objective function increased by 40% (Fig. 16, step 2).
Step 3. The final step consisted in alleviating this unwanted effect. The whole matching workflow (Fig. 1) was run again starting from the model determined at step 2. Ten new gradual deformation parameters were introduced to vary the spatial distribution of facies IE and HS in units MT2_3B and MT2_1&2 as well as the spatial distribution of porosity within facies IE. At this stage, the data to be matched encompassed both the available production and 4-D inverted seismic data. The minimization process was repeated until the two terms of the objective function achieved the values given in Figure 16, step 3. The seismic match kept essentially unchanged but the main outcome was the 33% decrease in the production data mismatch term. Examples of improvement are shown in Figure 15c for the water flow rate at well P2, in Figure 17a for the RFT pressures at well P1 and in Figure 17b for the bottom hole pressures at well I1.
These results illustrate the ability of the proposed matching workflow to handle both production and 4-D inverted seismic data. They also point out the issues related to the number of parameters and the significant uncertainty in data, mainly 4-D seismic data. Improvements could have been probably obtained by adjusting the petrophysical properties within the other MT2 units and by refining the parameterization techniques (for instance, by varying the mean properties over specific regions). A typical aspect of the case studied is the poor quality of the 4-D inverted seismic data. This obviously Step 1 -Model constrained to production data only Step 2 -Adjusted model Step 3  RFT pressure at well P1 after 527 days of production and bottom hole pressure at well I1. The reference data are plotted in red, the corresponding responses simulated for the model constrained to production data only in yellow and the ones for the model constrained to both production and 4-D inverted seismic data in blue.

Figure 16
Evolution of the objective function for the three successive steps considered in this study. Purple bar: production data mismatch term. Blue bar: seismic data mismatch term.
made the matching of these data more difficult. For now, the seismic data mismatch was expressed following the usual least-square formulation. Alternative metrics could help to reduce further this mismatch.

CONCLUSION
The simulation workflow introduced in Section 1 is based upon two scales: a fine geological scale and a coarse reservoir scale. At the fine scale, the dimensions of grid blocks are typically 0.5-1 m in the vertical direction and up 10 to 50 m in the horizontal direction. At the coarse scale, they are a few meters in the vertical direction and up to 100 to 200 m in the horizontal direction. Such a design was favored to make data integration at the proper scale. In addition, considering data of different types like production and 4-D seismic data led us to revise the usual history-matching practices. The formulation of the objective function and the development of customized parameterization techniques turned out to be two key factors controlling the efficiency of the matching process. The potential of the proposed matching methodology was addressed through two examples. The first one, which refers to the synthetic PUNQ-S3 case, highlighted the added value of 4-D seismic data for improving the reliability of reservoir models and consequently of the forecasts estimated from these models. The second example was a field case operated by Statoil. Again, we applied the workflow developed in the first section to determine a model respecting as well as possible production data and acoustic impedances. This second study stressed the difficulties to be faced when dealing with a real field case. These ones can be related to the uncertainties in seismic-related data as well as the uncertainties in modeling.
Restricting the number of scales to two in the matching workflow proposed is clearly a simplification. For instance, one can argue that the seismic scale differs from the geological scale: the vertical seismic resolution is a dozen meters. Thus, the simulation workflow should be supplemented by an additional upscaling step. A more tricky issue is about the use of core and log data to constrain the geological model. A core plug is only a few centimeters long while logs investigate areas over 10 centimeters to a few meters. Rocks can be very heterogeneous at the sub-geological scale because of laminations or bed bounding surfaces and this can have a significant influence on flows. Therefore, some form of upscaling is also necessary to transfer core and log information towards the geological model. Some of the reasons why this step is usually avoided are listed below (Pickup and Hern, 2002). Flows in reservoirs can be mainly driven by large scale features like faults or channels, so that the core/log scale can be considered as negligible. There may be not enough data to build a reliable model at the core/log scale. Even if there were enough data, the resulting model would comprise a few billions of grid blocks. This number would have to be drastically reduced before performing flow simulations, calling again for some appropriate upscaling. This upscaling is time-consuming and is not guaranteed to properly capture the equivalent properties to be input into the coarse reservoir model. This is all the more true for multi-phase flow, which is probably more demanding in terms of describing heterogeneity.
Clearly, there is a demand for more sophisticated techniques capable of accounting for multiple scales when building geological models, performing flow simulations and historymatching. Multiscale methods may be among the most promising avenues for further improvement.