History Matching of Production and 4D Seismic Data: Application to the Girassol Field, Oﬀshore Angola.

Résumé — Calage simultané des données de production et de sismique 4D : application au champ de Girassol, Offshore Angola — La sismique répétée (4D) constitue une source d’information précieuse sur l’évolution dans l’espace et le temps de la distribution des hydrocarbures dans les réservoirs pétroliers. Le monitoring sismique améliore notre compréhension des mécanismes de production et permet d’optimiser la récupération des hydrocarbures. Bien que les données de sismique 4D soient de plus en plus utilisées par les compagnies pétrolières, leur exploitation est souvent qualitative, en raison du manque de techniques d’interprétation appropriées. Des travaux de modélisation récents ont démontré la faisabilité de l’intégration des données de sismique 4D pour mettre à jour les modèles de simulation de réservoir. Cependant, les méthodologies basées sur l’interprétation séquentielle des données 4D, par essais et erreurs et tests successifs de simulations de réservoir, exigent un effort conséquent et la mobilisation d’équipes pluridisciplinaires. Le développement actuel des méthodes de calage Abstract — History Matching of Production and 4D Seismic Data: Application to the Girassol Field, Offshore Angola — Time-lapse seismic provides a source of valuable information about the evolution in space and time of the distribution of hydrocarbons inside reservoirs. Seismic monitoring improves our understanding of production mechanisms and makes it possible to optimize the recovery of hydrocarbons. Although 4D seismic data are increasingly used by oil companies, they are often qualitative, due to the lack of suitable interpretation techniques. Recent modeling experiments have shown that the integration of 4D seismic data for updating reservoir flow models is feasible. However, methodologies based on sequential interpretation of 4D seismic data, trial and error processes and fluid flow simulation tests require a great effort from integrated teams. The development of assisted history matching techniques is a significant improvement towards a quantitative use of 4D seismic data in reservoir modeling. This paper proposes an innovative methodology based on advanced history matching solutions to constrain 3D stochastic reservoir models to both production history and 4D seismic attributes. In this approach, geostatistical modeling, upscaling, fluid flow simulation, downscaling and petro-elastic modeling are integrated into the same history matching workflow. Simulated production history and 4D seismic attributes are compared to real data using an objective function, which is minimized with a new optimization algorithm based on response surface fitting. The gradual deformation method is used to constrain the facies realization, globally or locally, which populates the geological model at the fine scale. Moreover, a new method is proposed to update facies proportions during the optimization process according to 4D monitoring information.


INTRODUCTION
Time-lapse seismic monitoring is a powerful technique to optimize field development. This additional source of information is particularly useful in deep offshore environments, when fields are produced from a limited number of wells with a complex architecture. As the distribution of heterogeneities and the production mechanism are very uncertain, 4D seismic data help to better capture the fluid flow behavior during production, and improve our knowledge of the geological framework.
Combined with reservoir modeling, time-lapse seismic monitoring enables reservoir engineers to improve reservoir characterization and reduce uncertainty in production forecasts.
The feasibility of using 4D seismic information to update the flow model has been already addressed in a number of studies. As an example, 4D seismic has been used successfully to update the Draugen field reservoir model (Guderian et al., 2003) and make predictions more reliable. For the Girassol field, 4D seismic has been used to rework the reservoir model and better locate the injected gas cap (Jourdan et al., 2006;Gonzalez-Carballo et al., 2006). For the Gullfaks reservoir management (Talukdar et al., 2008), 4D seismic has been tremendously helpful in getting a strong history match of the simulation model. The methodology involved the comparisons of saturation data with seismic responses using visualization tools. Another example stressing the added value of 4D seismic in reservoir management is given by Oliveira et al. (2007) for the Marlim field: incorporating de sismique 4D ont aussi contribué à mieux caractériser la distribution spatiale des hétérogénéités dans le réservoir. Le modèle géologique détaillé a ainsi été amélioré pour être cohérent avec la simulation des 4D seismic led to a better characterization of the geological model in a Turbidite system and helped history matching using secondary information. An integrated workflow for the processing of 4D seismic data is also presented by Lerat et al. (2010) for the monitoring of steam-assisted gravity drainage production of a heavy oil field.
However, the use of 4D seismic data often remains qualitative and subject to interpretation. Revisions of the reservoir model are usually made after processing and interpreting the monitor seismic, given some assumptions about the fluid flow mechanism. The multidisciplinary nature of 4D monitoring calls for combining the seismic, geologic and reservoir engineering know-how for a more reliable interpretation of 4D data. A forward interpretation of 4D seismic data without sufficient knowledge of the fluid flow behavior in the field may lead to indetermination: pressure effects may be difficult to discriminate from saturation changes for example. A better understanding of fluid changes requires iterating between reservoir characterization, fluid flow simulations and interpretation of 4D seismic data.
The development of assisted history matching techniques with 4D seismic data is a very active domain (see Gosselin et al., 2003;Lygren et al., 2005;Skjervheim et al., 2007;Stephen et al., 2008;Peters et al., 2010). It is a significant improvement towards a quantitative use of 4D seismic data in reservoir modeling. The approach described in this paper is based on a fully integrated assisted history matching process. The fine scale geological model and the fluid flow reservoir model are considered in the same inversion loop (see Mezghani et al., 2004;Langlais et al., 2005). Interpreted time-lapse seismic data are used to constrain the reservoir model and update the fine scale geological model in a quantitative way. As a result, a better reservoir characterization, and thus more reliable production forecasts, can be expected. This paper presents a successful application of these advanced methods for history matching production and 4D seismic data to the Girassol field. The main objective was to test the feasibility of the integration of 4D seismic based on assisted history matching techniques on a real field. The long-term perspective is the development of an operational monitoring methodology to facilitate the exploitation of new reserves located, in particular, deep offshore or in complex geological environments.

HISTORY MATCHING WITH 4D SEISMIC DATA
Currently, history matching is often performed after upscaling, at the fluid flow simulation scale only. The fine-grid geological model is hardly ever updated after history matching due to the lack of appropriate tools and methodologies. In the approach described in this paper, the geostatistical model and the fluid flow model are simultaneously updated during history matching. This approach entails the incorporation of the entire simulation workflow in the inversion loop, from the generation of the geological model to the numerical fluid flow simulation. Upscaling must be also included in this process to account for the different scales: the geological model is associated to a scale finer than that of the reservoir model.
History matching needs a high degree of flexibility in model parameterization, as history matching parameters can be at any level of the workflow. An essential aspect of the approach described hereafter is that it makes it possible to adjust various parameters of the geostatistical model like global or local modifiers of the facies realizations (on the basis of gradual deformations, see Hu, 2000), structural geostatistical parameters (facies proportions, correlation lengths, anisotropy directions) or petrophysical properties (permeability, porosity, etc.) populating a given facies realization. Usual history matching parameters such as aquifer strengths, relative permeability end-points or fault transmissivities can also be constrained simultaneously with the geostatistical model parameters.
The purpose of the Monitor approach (Mezghani et al., 2004) is to enhance the characterization of subsurface reservoirs by combining geological, geophysical and reservoir engineering data. It is based upon the integration of timelapse seismic attributes in the history matching workflow.
In this workflow, simulated production responses and variations with time of compressional and shear impedance responses are computed from a multi-phase fluid flow simulator coupled with a rock physics model based on Gassmann equations (Gassmann, 1951). The fluid flow simulator is used to simulate the pressure and saturation changes with times. These dynamic variables are then downscaled before being provided as inputs to the Petro-Elastic Model (PEM). They are used to compute synthetic compressional and shear impedances (Ip, Is) at the geological model scale. The following step consists of filtering the synthetic impedances in the bandwidth of the inverted seismic data. The mismatch between real data and numerical responses is estimated from a weighted, least squares objective function including both production and 4D seismic data. Finally, an optimization process is used to minimize the objective function and adjust the selected model parameters at the different scales.
In this iterative workflow, two complementary parameterization techniques are used for updating the reservoir model realization at the fine scale when geostatistical simulations have been used to describe the spatial distribution of heterogeneities: the gradual deformation method and the facies proportion calibration method.

Gradual Deformation Method (GDM)
Geostatistical modeling is now increasingly used for building detailed geological models which represent at best all the available static information. Well data, seismic data and geological information can be integrated to model the spatial distribution of heterogeneities consistently. Several realizations of facies or petrophysical property distributions may be generated to account for uncertainty.
However, in most cases, fluid flow simulations performed with these model realizations do not reflect the dynamic behavior observed on the real field. Reservoir engineers also need to modify these realizations to match dynamic information like well-test data, production history or 4D seismic data. The Gradual Deformation Method (GDM) makes it possible to adjust model realizations and improve history matching while preserving the overall geological structure (Hu et al., 2001;Roggero et al., 1998). The resulting deformations applied to the initial realization can be global or local, and are driven by a reduced number of parameters.
The main idea behind the Gradual Deformation Method is to control the stochastic simulation process from a continuous parameter, called the deformation parameter, instead of a random seed. Whatever the value of the deformation parameter, it corresponds to a new realization of the reservoir model. In the Gaussian framework, a convenient way of building a continuous chain of model realizations consists in combining two independent standard Gaussian random functions with identical covariances.
The Gradual Deformation Method, initially limited to Gaussian models, was generalized to the deformation of any kind of Gaussian-related model (e.g., lognormal model, truncated Gaussian model, etc.) using appropriate variable transformations. It was also extended to non-Gaussian stochastic simulations like the sequential indicator simulation and Boolean simulations (see Hu, 2000). Figure 1 shows an example of a continuous chain of realizations. In this example, a facies model realization of the Girassol field is gradually deformed, resulting in slight changes in the spatial distribution of facies in sandy channels.

Facies Proportion Calibration Method (FPCM)
The Gradual Deformation Method makes it possible to directly modify the geostatistical model realization throughout the matching process, while keeping the average geostatistical properties unchanged (variogram and facies proportions). However, in some cases, changing the geostatistical realization or the petrophysical properties is not enough for matching production and 4D seismic data. Uncertainty may be associated to facies proportions, and proportion trends may significantly impact the dynamical behavior of the field.
The methodology applied to the Girassol case  accounts for the parameterization of facies proportions, locally or globally, in order to constrain the model by production and 4D seismic data. The interested reader can refer to Ponsot-Jacquin et al. (2009) andTillier et al. (2010) for more details. This Facies Proportion Calibration Method (FPCM) is complementary with the Gradual Deformation Method: it makes it possible to drive the average spatial trends whereas gradual deformations affect the local distribution of heterogeneities with fixed facies proportions.
This approach consists in adjusting the average proportion of a given facies group with respect to the average proportion of a larger facies selection. Transformations of facies proportions can be applied globally or inside user-defined reservoir zones. The algorithm includes the following steps: -define reservoir zones in which facies proportions have to be constrained; -select the facies to be varied by the transformation, to form an entity called "selection"; -within a selection, associate a sub-group of facies to define a second entity called "association"; -define a global parameter as the proportion ratio between facies in the association and facies in the selection. Once the parameters listed above are defined, the algorithm makes it possible to fix or calculate new average proportions for the selected facies. Several zones with different transformations can be defined. The transformations are sequentially calculated following the order in which they were defined.

APPLICATION TO THE GIRASSOL FIELD
In this project, the main objective was to demonstrate the feasibility of integrating time-lapse seismic data to quantitatively improve reservoir characterization, using the proposed history-matching approach to the Girassol field data.
The project followed the workflow presented in Figure 2. Starting from an existing structural model, a detailed geostatistical model constrained to the base 3D seismic data was first built . The motivation was to obtain an accurate description of reservoir heterogeneity at various scales and good agreement with the High Resolution (HR) base seismic data before addressing the integration of timelapse seismic information. To reach this goal, a lithoseismic interpretation of the 3D HR base seismic data was conducted  Workflow used to build a detailed reservoir model of the Girassol field and perform the history matching of production and 4D seismic data.
to derive a non-stationary grid of facies proportions. In addition, the resulting facies proportions were adjusted so as to improve the match between the synthetic compressional impedances computed from the geological model and the real impedances obtained from seismic. Reservoir characterization uncertainty was accounted for on the basis of stochastic simulation: several facies model realizations were generated, all constrained to seismic and well data. These facies model realizations were populated with petrophysical properties (porosity, permeability) estimated from well data. Multiple realizations were used first to assess initial uncertainty on sand proportion and connectivity  and second during history matching when using the Gradual Deformation Method for constraining the model to dynamic information.
In the dynamic simulation workflow, a numerical uspcaling method was implemented to preserve the impact of heterogeneity at the coarser scale of fluid flow simulation . The integration of dynamic information was performed in two steps using an iterative optimization loop. In a first step, history matching was achieved with production data only. A limited period of the available historical data was used to keep the additional data for validation of simulated production forecasts. In a second step, 4D seismic data were used simultaneously with production data to update again the reservoir model. Finally, production forecasts simulated with and without integration of 4D information were compared to conclude about the added value of 4D data.

Girassol Field Presentation
The Girassol field is a complex and faulted turbidite field located in Block 17, deep offshore Angola. The reservoir is composed of good quality oil-bearing Oligocene channel-levee complexes and sand sheets extending over an 18 × 10 km area. Channel complexes, with general NE-SW and N-S orientations, encompass sinuous to meandering, stacked elementary channels and associated levees (see Fig. 3).

Channel margin and levee laminated sands
Channel fill sands Hemipelagites Debris-flows Raisson andTemple, 2004 Bouchet et al., 2004 La ter al mig rat ion : ag gra da tion rat io NW SE

Figure 3
Overall Girassol reservoir architecture (top view) and sedimentological model of lateral offset stacked channels (bottom view).
The base 3D High Resolution (HR) seismic survey was shot in 1999. Field development started in 2000 and the first oil was produced in December 2001. The field was initially close to the bubble point pressure with no gas cap. After three appraisal wells, the decision was made to launch a fast-track development.
Production gas was re-injected in the top of the reservoir, in a geological environment formed by lateral offset stacked channels (LOSC, see Fig. 3). A strong depletion was observed in the first months of development. Due to high connectivity in sand channels, released gas appeared in extended regions depleted by production wells. Water injection was initiated after 6 months of production to maintain reservoir pressure.
4D HR seismic was planned from the very beginning of the development, mainly to monitor the gas bubble extent, which is subject to uncertainty in this highly heterogeneous environment. The first repeated 3D HR survey was shot in December 2002, after only one year of production and six months after starting gas injection. A second monitor survey was shot in 2004.
Before the first monitor survey, water from aquifers was collected in two production wells, one in the north and one in the south of the field. A gas breakthrough was also observed during the first production year at one producer, south of the gas injection well. Then, after the first monitor survey, gas breakthroughs were observed in other wells in the north and injected water reached other production wells.
History matching was performed up to the first monitor survey (2002) with one year of production history (during 404 days, starting in December 2001). The production history included data for 12 production wells and 4 injection wells. The production data collected after 2002 as well as the second seismic monitoring survey (2004) were used at the end of our study to compare the simulated production forecasts with real data, hence to evaluate the reliability of our matched model in terms of predictions.

Geostatistical Geological Model
The aim was to build a detailed model of the Girassol field to describe the spatial distribution of heterogeneities at the fine scale. In the horizontal X and Y directions, a relatively short grid cell size was adopted (33 meters) as a compromise between an enhanced preservation of seismic information and the number of cells the software can reasonably handle. The vertical grid cell size in the channel complexes was set to approximately two meters to better represent the sub-seismic heterogeneity identified from well data. The resulting fine grid comprises about 28 millions of cells. The fluid flow simulation grid was defined by aggregating cells two by two vertically and three by three horizontally. Compelling this model to respect the base seismic data at the fine scale was considered as a prerequisite before going through the integration of time lapse data. A dedicated methodology was developed to allow the integration of seismic information within the geological stochastic modeling workflow. The modeling workflow was made of three steps: -Step 1: lithoseismic interpretation and definition of a 3D seismic constraint; -Step 2: reconciliation of seismic data with the geological model; -Step 3: geostatistical geological modeling with a 3D seismic constraint and well data. The result of the first step was a 3D volume of geological facies proportions. Seven lithofacies were considered after the analysis of a database of 40 wells including core data, sedimentological descriptions, laboratory petrophysical measurements and available logs . In the 3D volume, the inverted seismic data were used to perform a lithoseismic interpretation and predict six seismic facies associated with probability cubes (Nivlet et al., 2005a). This seismic information was then downscaled and combined with the lithofacies analysis to define the 3D seismic constraint made of seven cubes of facies proportions. As an example, Figure 4 shows two horizontal slices extracted from the volumes of average proportion of shale and sandstone, respectively. The good quality of reservoir volume delineation provides insights about the capability of seismic-derived information to constrain the geostatistical realizations.
The second step aimed to improve the quantitative match of 3D seismic attributes computed directly from the detailed description of the facies proportion model. This implied first ensuring optimal consistency between well and seismic data. Then, distributions of acoustic parameters were determined for each geological facies on the basis of well log data. This information was used to compute 3D synthetic cubes of seismic attributes, by combining average values of acoustic parameters with facies proportions. Finally, real and synthetic seismic attributes were compared and the resulting mismatch was minimized from an optimization algorithm by varying the 3D facies proportions. The comparisons between the predicted and real impedance cubes showed the capability of the geological model to reproduce High Resolution seismic information  as illustrated in Figure 5.
The third step consisted of populating the fine grid geological model with facies referring to the non-stationary truncated Gaussian method (Galli et al., 1993). This stochastic simulation process accounted for well data and the 3D volume of proportions derived from the HR 3D base seismic .
This methodology was successfully applied to the Girassol field. In the model realization shown in Figure 6, the seismic constraint defines the main spatial trends of facies proportions. In addition, the wells locally constrain the position of facies in the model. The modeling process captures Average shale a) and sandstone b) proportion models.

Figure 5
Cross-sections comparing the synthetic seismic impedances with the real impedances derived from seismic inversion.
the geometry of complex geological objects and reproduces the distributions of heterogeneities. As a result, the initial geological model preserves the geological organization in terms of channel complex and channel story.

History Matching of Girassol Production Data
In this first history matching phase of the Girassol project, production data only were considered while 4D seismic information was disregarded. Approximately one year of production history was considered, up to the first seismic monitoring survey. During this period, water or gas breakthroughs were observed at 3 wells only (see water cut and GOR data in (Fig. 7): water production was observed at PROD1 and PROD2 and gas production at PROD3. Note that water and gas breakthroughs were observed at other wells after the first 4D seismic survey, although this information was not used for this history matching study.
The production data to be matched included: -the static shut-in pressures recorded for 12 wells; -the pressure profiles recorded by Modular Formation Dynamics Tester (MDT) tools along well trajectory for 8 wells; -the dynamic flowing pressures measured in 12 wells; -the Gas Oil Ratio (GOR) data for the PROD3 well with gas breakthrough; -the water cut data for wells with a water breakthrough (PROD1 and PROD2). The assisted history matching was achieved by minimizing an objective function using an iterative optimization process. The weighted least-squares formulation was used to define the objective function as follows: (1) where: -nseries is the number of production data series to be matched. Each set corresponds to a type of data collected at a given well; -ntime j is the number of measurement times for the data series j; -P i,j obs is a production datum observed at time i for the data series j; -P i,j sim is the datum simulated for the same time; -σ i,j P is the standard deviation on production data errors; -w j P are additional weighting coefficients assigned to data series. Weighting coefficients were defined according to the objectives of the different history matching steps. The Powell's "dogleg" algorithm (Powell, 1970) was used to minimize the objective function, that is to minimize the production data mismatch. This gradient-based optimization algorithm is known to be efficient for solving inverse problems. In our approach, approximated derivatives of fluid flow simulation results with respect to the selected parameters were numerically computed using finite-differences. Additional simulations with perturbations applied to parameters were automatically managed by the optimization loop at each iteration. A polynomial model was then used to fit the simulation results and compute derivatives. The number of simulations to be performed was determined during the optimization process referring to an optimality criterion.

Workflow
An integrated simulation workflow was implemented to automatically chain geostatistical simulation, upscaling and fluid flow simulation in the same inversion loop. A facies realization is generated first using the truncated Gaussian approach. Then, facies are assigned average petrophysical properties (porosity, permeability, irreducible water saturation) Geostatistical realization describing the spatial facies distribution. It was simulated using the truncated Gaussian approach with non-stationary facies proportions derived from 3D HR seismic data.
computed from well data. Finally, permeability and porosity distributions are upscaled from the geostatistical stratigraphic grid to the reservoir grid and the fluid flow simulation is launched.
As the Girassol field is very heterogeneous, the upscaling method was carefully selected to preserve local flow barriers or sand connectivity at the coarse scale (see Ding et al., 2007).
This integrated assisted history matching approach enables the reservoir engineer to select inversion parameters in any component of the workflow. The geostatistical model can be updated at the fine scale using the GDM or the FPCM. In addition, regular history matching parameters can be defined in the fluid flow simulation data set. To maintain consistency between the various modeling scales, the entire workflow is simulated again as far as parameters are changed by the optimization loop.
For the Girassol study, the following history matching parameters were investigated: -gradual deformation parameters to control the facies realization in each reservoir unit; variogram ranges to account for uncertainty on geostatistical parameters; average petrophysical properties assigned to each facies in each simulation unit (porosity, horizontal and vertical permeability values); cut-off parameter based on shale proportion (Vclay) to remove non-reservoir grid-blocks at the fine scale; aquifer parameters like porosity and permeability; -fault transmissivity multipliers; -well skin factors. As there may be numerous inversion parameters, history matching was conducted in several steps with different sets of parameters. Sensitivity studies and reservoir engineering analysis were performed to select the most appropriate set of parameters at each optimization step.
In a first step, static pressures were matched to constrain various parameters either at the geological or reservoir model scale: volumetric-related parameters (cut-off values and facies porosities), facies permeabilities, facies realization and aquifer strengths. MDT pressure profiles recorded along well Water cut and GOR data (left) and well positions (right).
trajectories were introduced to constrain the hydraulic communication between complexes. In a second step, multiphase production data were used to match water and gas breakthroughs: water cut data for the PROD1 and PROD2 wells, and GOR data for the PROD3 well. Intensive gradual deformations were performed in order to constrain the model realization to the dynamic data. Other regular parameters like the facies-related petrophysical parameters were also optimized during history matching.

History Matching Results
A fluid flow simulation performed with the initial model showed a too fast pressure decrease compared to the observed data. It was mainly explained by the initial permeability values of sand facies assumed to be underestimated. These values were uncertain as they were measured on unconsolidated samples. The model realization also significantly contributed to hydraulic communication. Thus, the most influential parameters to match pressure data were sand facies permeabilities as well as gradual deformation parameters. Sand permeabilities were multiplied by about 2.5 after history matching to obtain a global improvement, while global gradual deformations made it possible to improve the pressure match well by well.
A good pressure match was obtained for most of the wells. An example is presented for a representative well in Figure 8 (bottom right), where initial and optimal simulation results are compared to real data. Static pressures are matched at the end of build up tests, which are simulated using refined time steps during well closures. Note that a strong pressure drop is observed in the reservoir during an initial period of about 200 days, and that pressure is maintained by water and gas injection afterwards.
The gradual deformation of facies realizations strongly impacted water and gas breakthroughs. As shown in Figure 8, a significant improvement is obtained for the water cut and GOR match, although a delay of about 60 days remains on water breakthrough in the PROD1 well. A satisfactory match of GOR data is obtained for the PROD3 well, even if further improvements should be expected after 4D seismic data matching.
The influence of the gradual deformation process on gas production is illustrated in Figure 9, which shows gas saturation maps in a cross section between the gas injector and production well PROD3. With the initial model realization, a shale barrier prevents the injected gas from reaching the producer, which is located in a lower sequence. After history matching, a vertical communication is created thus enabling the match of the observed gas breakthrough. Simulation results before and after history matching: water cut data (top views), GOR data (bottom left) and static pressure (bottom right).

History Matching of Girassol 4D Seismic Data
In the Girassol field, the first repeated 3D HR survey was shot after only one year of production and about six months after the start of gas injection. It was motivated by the need to monitor gas injection in an extremely heterogeneous environment. The excellent quality of the 4D response and of further processing provides useful information to update reservoir models (Jourdan et al., 2006;Gonzalez-Carballo et al., 2006;Gosselin et al., 2003).
In this section, we present the application of the Monitor methodology, which allows for history matching 4D seismic attributes together with production data. The objective is to validate the feasibility of a quantitative use of 4D seismic data and to evaluate the added value of 4D seismic in the reservoir modeling workflow. Starting from the model matched with production data only, it is expected that the integration of 4D seismic data will improve the reliability of predictions.
After a presentation of the 4D seismic data set, the history matching workflow and results are described in this paper.

Time-Lapse Seismic Data
Seismic data used for this Girassol field study consists of two time-migrated angle sub-stacks ([3°-22°] and [23°-37°]) acquired in 1999 (base survey), before starting production, and in late 2002 (monitor survey), after nearly one year of production. Both surveys cover the same area (about 100 km 2 ), and were collected and processed with a particular workflow to ensure the optimal repeatability of the 4D signal (Lefeuvre et al., 2003). The monitor pre-stack amplitude cubes were time-corrected (or warped) by computing a relative instantaneous velocity cube variation ΔVp/Vp (Williamson et al., 2007).
A joint pre-stack stratigraphic inversion (Tonellot et al., 2001) was applied to invert the base and monitor surveys (Nivlet et al., 2005b). As the low-frequency component of the impedance model is not updated during the inversion process, particular attention was paid to the construction of the a priori model. The a priori model used for the inversion of the base survey was updated to perform the inversion of the monitor survey (Nivlet et al., 2006). Relative changes in low-frequency P-impedance were determined using the lowfrequency component of the relative velocity cube variation ΔVp/Vp. The seismic inversion allowed us for estimating (P-and S-) impedance models explaining the observed seismic amplitudes (Nivlet et al., 2006). P-impedance differences between monitor and base surveys were computed to obtain P-impedance variations (ΔIP). An example is showed in Figure 10: it consists of a horizontal Gas saturation in a cross-section between injection and production wells, for the initial model (top view) and the history matched model (bottom). slice extracted from the channel complex area. Gas injection from well INJ1 (in purple), which induces a significant P-impedance decrease, can be seen at first glance in zone A. The extension of this negative 4D anomaly follows the path of good reservoir zones which correspond to high proportions of sands. Facies proportions are obtained from lithoseismic interpretation (Nivlet et al., 2005a). Finally, lateral extension of this 4D anomaly is also consistent with the posterior field history, since gas breakthrough was recorded only a few weeks after the monitor survey acquisition at production wells (in green) located in the North. In zone B, P-impedances increase around water injectors (in blue), which is consistent with the replacement of oil by water. Finally, in the South, the negative 4D anomaly close to a water injector (zone C) can be related to a pressure increase linked to water injection in the aquifer, combined with variations in salinity and temperature.
The 4D seismic attributes used for history matching are the P-impedance differences between the monitor and the base surveys , combined with the low frequency component of the relative instantaneous velocity cube variation ΔVp/Vp (Williamson et al., 2007). These attributes, obtained from the seismic inversion of time lapse data, were upscaled in the stratigraphic grid.

Modeling of 4D Seismic Attributes in the Workflow
In the Monitor approach, (Mezghani et al., 2004), the integration of 4D seismic data is performed by considering time lapse seismic attributes as part of the overall history matching process. Thus, this workflow has to handle the computation of synthetic seismic attributes that can be compared with the real data. To achieve this goal, a Petro-Elastic Model (PEM) is implemented to calculate time-lapse seismic attributes with respect to pressure and saturation changes obtained from the fluid flow simulation.
In the Monitor workflow (Fig. 11), several scale differences have to be taken into account. Fluid flow simulations are performed on a grid coarser than the fine stratigraphic geological grid. The PEM is applied at the fine scale to account for a detailed description of lithotypes and avoid upscaling of petro-elastic parameters. Finally, modeled seismic attributes include high-frequency information due to the heterogeneity of the geological model, and cannot be compared at once to real seismic data which are characterized by a limited frequency bandwidth.
The Monitor approach involves the steps listed below to address properly these scaling issues in the inversion loop (see also Fig. 11). These steps complement the ones already  -filtering of computed attributes in the seismic data bandwidth; -computing of an objective function based on a leastsquares formulation to quantify the mismatch between computed and real seismic attributes; -updating of reservoir model parameters at the fine-scale using an optimization algorithm to minimize the objective function.
The PEM used for the Girassol study is based on Gassmann equations (see Mavko et al., 1998;Gassmann, 1951) for modeling the fluid substitution effects. The pressure effect is modeled by proprietary correlations from Total, calibrated on laboratory and field data. Petro-elastic properties (bulk modulus or shear modulus) are associated with the different geological facies of Girassol, essentially for the pore volume variation with stress changes. The suitability of this PEM was checked with static well data and provided satisfactory results.
In the proposed workflow, the petro-elastic modeling is performed at the fine scale of the geological model. Therefore, pressure and saturation changes must be downscaled from the coarse fluid flow simulation grid to the fine grid. For the Girassol model, a coarse grid-block is an aggregation of several fine grid-blocks. The downscaling is performed by assigning a constant value extracted from each coarse gridblock onto the corresponding fine grid-blocks. The downscaling of saturation values takes into account irreducible water saturation end points defined at the fine scale as well as a pseudo-Vclay cut-off applied to the fine grid porosity to remove non-reservoir grid-blocks.
The filtering is performed on P-impedance differences ΔIP computed using the PEM. As a low frequency term was added to the measured ΔIP attribute to be history matched, only high frequency components must be filtered. According to the seismic data bandwidth, the band-pass filter is characterized by the frequency domain (0-0, 120-140) Hz.
One difficulty when comparing real and computed 4D data relies in the noise associated with measured data. This noise may have a negative impact on the objective function, which is based on differences between real and modeled data. Thus, a preliminary treatment of measured data noise is recommended to get a more reliable objective function evaluation. A median filtering in horizontal layers was used for this study to eliminate as much as possible the random component of 4D P-impedance measured data. Reservoir model

Figure 11
History matching of production and 4D seismic data based on the Monitor approach.

History Matching Steps
In this project phase, production data and 4D seismic attributes were history matched simultaneously, starting from the optimal model constrained with production data only. Additional difficulties arose when matching 4D seismic data. The objective function must be defined carefully, with an appropriate weighting of seismic attributes compared to production data. Moreover, as 4D attributes provide much more spatial information compared to well or production data, more flexibility is required in the model parameterization. As an example, time-lapse seismic provides information related to the geometry of geological sequences: the extent of reservoir bodies, geometry of preferential fluid paths, position of flow barriers, etc.
This specific difficulty motivated the investigation of new inversion parameters for matching 4D data. More particularly, the Facies Proportion Calibration Method (FPCM) was successful in introducing additional constraints based on geology.
As the spatial distribution of heterogeneity also plays a significant role in the characterization of sand connectivity, the facies model realization was also history-matched using the Gradual Deformation Method (GDM). Global, then local deformations were performed to refine history-matching in the vicinity of particular wells.
Finally, the pressure effect as modeled by the PEM was calibrated to better fit 4D seismic attributes.

Objective Function Definition
Compared to the previous objective function F P (θ) used for matching production data only (Eq. 1), the seismic data contribution F s (θ) is added to the objective function when matching simultaneously production and seismic data: ( 2) where: -nregions is the number of reservoir regions where seismic data ΔIP must be matched; -nvalues j is the number of observed seismic data values for region j; -ΔIP i,j obs is an observed seismic datum i for data series j; -ΔIP i,j sim is a seismic datum computed for the same time; -σ i,j S is the standard deviation on seismic data errors; -w j S are additional weighting coefficients assigned to seismic data series. Particular attention was paid to the definition of the objective function weights. A good way of weighting ΔIP data in the least-squares approach should be to normalize the residuals between observed and simulated data by uncertainties (i.e., However, as seismic inversion was performed sequentially for the base and monitor surveys, this uncertainty is difficult to obtain. In our study, residuals on 4D data match in the objective function (i.e. differences between observed and computed data) were normalized by an average uncertainty range, estimated from the noise level in non-reservoir parts. In the history matching workflow, the contributions of 4D data were computed in five different regions, to compare their influences in different parts of the reservoir. For each region, the objective function contribution was normalized by the number of observation points, then weighted according to the extent of the region.
Finally, multiplying coefficients were applied to balance the influence of 4D seismic data compared to production data. This weighting is semi-empirical, as sensitivity studies should be performed to test the impact of inversion parameter variations on the objective function: comparable variation ranges should be obtained for both production and 4D seismic terms to match all data with an equivalent efficiency (in terms of convergence rate). For this study, the contribution of 4D seismic data was required to be multiplied by 16 after a series of preliminary tests.
Although the gradient-based optimization algorithm was very efficient in matching production data, it was not appropriate to minimize the objective function including 4D seismic data. The large amount of data to be matched makes it impossible to compute gradients for each data point. In the 4D seismic history matching case, a new optimization algorithm based on global adaptive learning of the objective function (Feraille et al., 2004) was used. The basic principle is to build a proxy model of the objective function, based on existing simulations, to compute new parameter estimates. The proxy model is successively improved following an iterative process until a convergence criterion is reached. At each iteration, new simulations are added according to different criteria: part of the simulations is selected to minimize the objective function, and other simulations are designed to better sample the parameter domain and avoid local minima.

History Matching of Facies Proportions using the FPCM
The history matching of facies proportions was motivated by the need of new geologically-based constraints to better match the 4D seismic data. Comparisons between simulated and real ΔIP data after history matching of production data showed that the extent of the gas bubble and the repartition of reservoir volumes are not sufficiently honored. Selecting the inversion parameters used for matching production data only was not sufficient to explain 4D seismic data, and a better analysis of main reservoir heterogeneity was needed to improve the match of 4D data with a satisfactory convergence rate of the optimization process.
Geological features like last-channel plugs, abandonment phases associated with high shale and silt contents, or channel margin collapses made of shaly debris-flows are potential barriers to fluid-flow. In the geostatistical-based modeling workflow of Girassol, facies proportions defined in the 3D stratigraphic grid play a significant role in the description of these geological features. Higher proportions of shale or sand facies in a given reservoir zone may lead to the introduction of additional constraints. As an example, increasing the shale ratio in last-channel plug regions will reinforce the continuity of local flow barriers. A similar method can be used to increase the sand ratio in channels and improve reservoir connectivity.
In the Girassol field, facies distributions are non-stationary and are strongly linked to geological regions. As an example, high proportions of sands are usually found in channels, while proportions of shale are generally high in plugs. Optimization regions in Girassol field were defined by combining different methods, as illustrated in Figure 12.
Pseudo-Vclay thresholds were used to identify regions with high proportions of sands or dominated by laminated and shale facies. Some architectural elements like abandonment channels, interpreted from seismic data by Total, were used to delimit flow barriers. Finally, ΔIP thresholds were used to identify regions associated with strong 4D responses.
A total of five parameterization regions were selected to drive the average facies proportions . The selected inversion parameters were the average proportion of sands in regions with good reservoir properties and the average proportion of shale in regions with the highest pseudo-Vclay.
A total of 60 inversion parameters were defined for history matching after refining these regions by depositional sequences. Successive optimizations with different sets of proportion parameters were performed to better match 4D data, leading to a total of 600 simulations for this phase.
The general trend was to increase both sand proportions in channel regions and shale proportions in levees and barriers, hence resulting in more contrasted images. This result is clearly visible on Vertical Proportion Curves (VPC) extracted from the model before and after optimization (Fig. 13). Shale sequences were reinforced while sand volumes are more concentrated in channel sequences.
A strong impact on the geostatistical simulation can also be observed on the facies realization shown in Figure 14. Sand channels are better characterized with sharper contours and facies continuity is improved in reservoir parts. Volumes of laminated facies are reduced in levees. Last-channel plugs, which are associated to almost 100% of shale after optimization, are now well depicted in the final model. Significant improvements were obtained in matching both production and 4D seismic data. Figure 15 shows observed ΔIP data (from seismic inversion) compared with simulation results before and after optimization of facies proportions.
Negative ΔIP values (red) are generally due to gas appearance or pressure increase. Gas injection effects can be clearly identified in the middle. The matching of the gas bubble extent is significantly improved after the optimization of facies proportions. Negative ΔIP values in the North and South areas are essentially related to released gas, because the bubble point pressure is very close to the initial reservoir pressure. Water injection effects can also be observed from simulations, where water is injected in oil zones. Positive ΔIP values (blue) can be explained by the replacement of oil by water. Some improvements are also found around water injection wells (blue zones are reduced after optimization of facies proportions).

Constraining the Facies Realization using the GDM
The inter-sequence heterogeneity distribution plays a major role in the fluid flow behavior in the Girassol field. The facies model realization was shown to strongly influence water and gas breakthrough times during history matching. As preferential flow paths may be modified by the model realization, gradual deformations were used in addition to the calibration of facies proportions for history matching 4D seismic data.
Global gradual deformations were performed first, in parallel with the FPCM. About 250 simulations were run, with again an improvement of both production and 4D seismic data matching.
Then, local gradual deformations contributed to refine history matching in given zones. Local gradual deformations Vertical Proportion Curves before and after history matching -Central area, close to gas injection well. Comparison in one stratigraphic layer of initial realization (left) and final facies realization (right) obtained after calibration of facies proportions.
make it possible to update the model realization in a target zone only, without changing the other part of the model (Hu, 2000;Le Ravalec et al., 2000;Gervais et al., 2007). This provides more flexibility in history matching water or gas movements in selected regions. The objective of local deformations was to improve water arrival from aquifers to producers or between pairs of injector-producer, and also gas sweeping from the injection well towards producers. For this purpose, a total of 6 zones were constrained using local gradual deformations. Time-lapse seismic data provide spatial information on fluid movements in these zones, between wells for which water or gas breakthroughs are not yet observed in the considered history matching period. Therefore, the expectation was an improvement of model predictions based on 4D seismic information. A total of 390 new simulations were performed in this history matching step.

Calibration of PEM Parameters
Two parameters of the PEM, named "Overpressure" and "Underpressure", were introduced as history matching coefficients. A history matching step was performed from the previous model realization to obtain their optimal values. Starting from nil parameters values, the optimal results were obtained with parameters equal respectively to 0.844 and 0.219. With the initial parameter values, only oil and gas compressibility and pore volume were accounted for in the modeling of pressure effects, and pressure had no influence in aquifers. With the optimal parameter values, a strong impact of pressure on velocity is obtained, and thus on the compressional impedance, Ip. P-impedance decreases when pressure increases. In water injection zones, pressure and water substitution effects are in opposition. Therefore, the interpretation of P-impedance variations (4D data) may be wrong without a calibration of the PEM.
Using the new PEM parameters, 4D impedance variations can be better simulated in aquifers or close to water injectors. The impact on 4D seismic data matching was quite significant. In higher pressure zones such as gas injection areas, negative ΔIP values were reinforced, while negative ΔIP effects were attenuated in lower pressure zones where released gas appears. Around water injectors, positive ΔIP values (due to substitution of oil by water) were either attenuated or turned into negative values.
Local gradual deformations were again performed to account for the new PEM parameters and to try to improve the water breakthrough and seismic data matching in these areas. About 135 simulations were performed in this step.

History Matching Results
The objective function values obtained for the main optimization steps are shown in Figure 16 and seismic data contributions are presented. At the end of history matching with production data only (first step shown after the initial simulation), the objective function has been greatly improved (by about 80%). An improvement is also found in 4D seismic impedance matching, although these data are not accounted for in the optimization. During the matching of 4D seismic data, the objective function is again reduced by about 35% for the production term and by about 15% for the seismic term. Each optimization step improves both the production and seismic terms: the facies proportion calibration first, then the global and local gradual deformations and finally the optimization of PEM parameters. Note the strong impact of PEM parameters in the last step. Despite the significant improvement obtained, the final value of the seismic term still keeps a high residual value. It may be explained by uncertainty or noise associated with 4D seismic data and by modeling errors. It also highlights the difficulty to match 3D data pixel by pixel using a leastsquares objective function. A more suitable formulation and a better understanding of uncertainties should improve the history matching results. Figure 17 presents some examples of wells to compare simulation results before and after the matching of 4D data. GOR data for the PROD3 well are perfectly matched with the final run. The match of water cut data for the PROD1 well is also significantly improved after the integration of 4D seismic data. Water cut data for the PROD2 well and static pressure data for the PROD4 well are quite well matched, either before or after the integration of 4D data.
The main features detected from the simulated flow behavior, such as the gas injection path, the extent of the released gas area or of the water injection zone turn out to be in good agreement with 4D seismic data. Simulated and measured 4D seismic data were compared in 3D volumes by computing geobodies corresponding to regions of connected grid-blocks with ΔIP lower than a threshold of -50 g/cm 3 .m/s. Figure 18 presents these geobodies in a slice of 10 layers located in the top of the main reservoir complex, in the region of gas injection. Highly negative ΔIP values may be related to the appearance of gas (either injected or released gas) or by a strong pressure increase. Connected volumes are clearly identified: again, the simulation results are consistent with 4D seismic data. The geometry of injected gas bubble extent is particularly well reproduced. The main differences are indicated in Figure 18 by circles. In the North, a channel by-pass is found, although not observed on real data. The match could be improved locally by the introduction of local barriers. Faulting may play a significant role to explain these differences. The main differences shown in the South part may be explained by artifacts of the PEM, due to the underestimated pressure effect in regions of water injection in the aquifer, or to the poor quality of the seismic constraint below the water oil contact. Figure 19 presents the low ΔIP geobodies in a slice of 10 layers located in the bottom of the main reservoir complex. This region corresponds to the erosion of a lower channel complex by the main reservoir complex (see schematic view of reservoir architecture in Fig. 3). It plays a significant role ΔIP simulated data layer #10 to layer #20 Regions ΔIP < -50 ΔIP filtered seismic data layer #10 to layer #20

Figure 18
3D geobodies corresponding to highly negative values of 4D P-impedance differences, for a slice of 10 layers in the top of the main complex. The color code represents the layer number.

Figure 17
Simulation results before and after history matching of 4D seismic data: water cut data (top views), GOR data (bottom left) and static pressure (bottom right). in the vertical hydraulic communication between reservoir complexes in the Girassol field. Although the extent of the simulated volume is smaller than the observed one, the geometry of geological bodies strongly impacted by production is in agreement with the observed 4D seismic data. The final results are also presented on a stratigraphic grid layer which corresponds to a thickness of about 1 m to 2 m (Fig. 20). Calculated 4D impedances (centre view) are compared to real seismic data (left view). The final facies realization is also shown (right view) to better explain the impact of heterogeneity on 4D effects. A strong correlation between sand facies and 4D effects can be noticed, which indicates that the integration of 4D data is very informative and helps to constrain the geological model in this kind of highly heterogeneous environment.
In the top of the main complex, the geometry of the gas bubble, which is monitored by the most negative ΔIP values, was captured quite well (Fig. 20). The underlying fine scale geological model helps in describing the observed data. However, in some regions High Resolution details may be lost in the computed 4D attributes. Upscaling and downscaling of properties in the workflow may lead to some averaging in the modeling results.

Added Value of 4D Data on Simulation Forecasts
Simulations were run to forecast production until the second HR monitor survey (2004) to compare simulation results with real data. Predictions were obtained after history matching with and without 4D seismic data, to evaluate the impact of time-lapse seismic information.
In the South area of the main reservoir complex, water breakthrough was observed for three production wells (Fig. 21). At well PROD2, water production observed during the history matching period was very accurately matched. The water-cut prediction for this well is quite good. For the two other wells, PROD4 and PROD5, water breakthrough was observed just after the history matching period. Simulation also predicts the production of water, but with a large delay compared with observed data. Nevertheless, the integration of 4D seismic data improves significantly the production forecasts at well PROD4 by reducing the breakthrough time of about 10 months. Prediction of water-cut for well PROD5 is also improved but to a lesser extent. In the North part, both good matches and good predictions are obtained for well PROD1 (Fig. 21). 3D geobodies corresponding to highly negative values of 4D P-impedance differences, for a slice of 10 layers in the bottom of the main complex. The color code represents the layer number.

Figure 20
Match of 4D P-impedance differences for the first monitor survey (2002) and facies realization in layer #30 (region of gas injection). Predictions of water cuts obtained after history matching with and without 4D seismic data.  Comparison of observed and simulated 4D data (ΔIP) in a region between water injection and oil production, second monitor survey (2004).

Figure 23
Prediction of 4D P-impedance differences for the second monitor survey (2004) in layer #21 (region of gas injection). Real data (left) and simulated forecast (right).
The analysis of 4D seismic prediction in the South zone confirms the observations from well data. The second monitor survey (Fig. 22) clearly shows a channel between the water injection well W_INJ1 and the production well PROD4. Positive P-impedance variations (color coded in blue) enable us to monitor the water front. Negative variations are explained either by pressure increase or the appearance of gas saturation. Observed data are not fully reproduced by simulation: details are lost in the simulation results, perhaps because of upscaling and downscaling issues. Water arrival can be seen on simulation results, but the water front is not driven fast enough towards the production well. Thus, the water breakthrough occurs later compared to the observed data.
The prediction of the second monitor survey (2004) is also shown in one layer of the stratigraphical grid which corresponds to the region of gas injection in the main reservoir complex (Fig. 23). The overall prediction of the geometry of the gas bubble extent is satisfactory compared to the observed data.
To better evaluate the added value of 4D seismic data when incorporated into the history matching process, predictions of the second monitor survey (2004) were compared before and after matching the information provided by the first monitor survey (2002). As shown in Figure 24, a significant improvement of ΔIP prediction is obtained after matching the first monitor survey. More specifically, it can be noticed that the 4D seismic anomaly is more contained by the contour of channels, in better agreement with observed data. This result may be attributed to an improved characterization of the facies geological model obtained when constraining the facies proportions by 4D seismic data.

CONCLUSIONS
An assisted history matching technique has been proposed in this paper to update reservoir models with 4D seismic data. This approach represents a significant improvement in the reservoir modeling process, as it makes it possible to use time-lapse seismic data in a quantitative way. Prediction of 4D P-impedance differences for the second monitor survey (2004) in layer #30 (region of gas injection). Real data (left) are compared with predictions obtained after matching production data only (middle) and after matching the 2002 monitor survey (right). Dedicated parameterization techniques have been used to constrain geostatistical geological models at the fine scale. The Gradual Deformation Method makes it possible to update the facies realization globally or locally. Moreover, a new method, called the Facies Proportion Calibration Method, has been applied to constrain the spatial distributions of facies proportions. Both methods, which are complementary, provide great flexibility in the history matching of 4D information to account for information varying in space and time.
The feasibility of this approach has been demonstrated with a successful application to the Girassol field. Using the Facies Proportion Calibration Method, additional geologically-based constraints have been introduced to improve the matching of 4D data. Average facies proportions were updated according to 4D seismic data, over reservoir regions based on seismic interpretation surfaces, seismic attributes and fluid flow considerations.
Better production forecasts and improved predictions of the new seismic survey, shot two years after the history matching period, have been obtained using the model constrained with 4D seismic data. The 4D seismic data has also provided an improved characterization of the spatial distribution of heterogeneities in the field. The underlying fine scale geological model was updated in cohesion with the fluid flow simulation model and observed data.
Building an initial geological model constrained by the base 3D seismic information, using seismic facies analysis and calibration of facies proportions, was a key element of success in this project. Particular attention was also paid to the upscaling of petrophysical properties to maintain a reliable characterization of flow barriers and reservoir connectivity from the detailed geological model to the coarse reservoir model. For the Girassol case, the numerical upscaling approach and near well upscaling have shown significant improvements compared to other algorithms.
Finally, some critical issues have been identified from the experience of this project: -the uncertainty associated with the PEM may have a strong impact on history matching results, and calibrations of PEM parameters are strongly recommended; -the computation time and number of simulations required may be large. However, using parallel computing of fluid flow simulations on several processors can greatly improve the elapsed time of Assisted History Matching processes; -the downscaling of pressure and saturations from the coarse simulation grid to the geological grid leads to a loss of resolution when trying to match very detailed 4D monitoring images. Improvements can be expected with more sophisticated flow-based techniques, but additional research is still required; -the definition of the objective function has to carefully take into account uncertainty on observed data and the relative weights of different sources of data.