Advanced Integrated Workflows for Incorporating Both Production and 4D Seismic-Related Data into Reservoir Models.

Résumé — Boucle de calage avancée pour construire des modèles de réservoir contraints par les données de production et les attributs sismiques — Les modèles de réservoir sont utilisés pour prédire les quantités d’huile qui peuvent être produites ou pour évaluer l’intérêt de différents scénarios de production. Ils sont considérés comme fiables lorsqu’ils respectent l’ensemble des données collectées au préalable sur le champ : les données statiques comme les données de diagraphie ou les mesures faites sur des échantillons de roche extraits des puits et les données dynamiques comme les pressions ou les débits. Depuis la fin des années quatre-vingt-dix, les données dynamiques comprennent aussi les données dites de sismique 4D ou de sismique répétitive. Les besoins en modélisation de réservoir ont motivé le développement de chaînes de simulation spécifiques qui permettent d’élaborer des modèles réservoir respectant toutes les données collectées. Cet article présente une chaîne de modélisation dédiée à la construction de modèles réservoir contraints à la fois par les données de production et les données inversées de sismique 4D. On parle ici de sismique inversée car Abstract — Advanced Integrated Workflows for Incorporating Both Production and 4D Seismic-Related Data into Reservoir Models — Reservoir models are used for predicting future oil recovery or for evaluating alternative field management scenarios. They can be considered as reliable when they account for all available data collected on the field: data are split into static data such as logs or measurements carried out on cores extracted from wells and dynamic data such as pressures and flow rates. Since the late nineties, the latter also consist of 4D seismic data. This motivated the development of very specific workflows, which yield reservoir models respecting all collected data. In this paper, we focus on workflows for building reservoir models consistent with both production and inverted 4D seismic data. Seismic data are referred to as inverted since we do not consider the amplitudes of the seismic traces, but the acoustic impedances or velocities derived from amplitudes. Then, two application cases are presented. The first one is a synthetic case inspired by typical North Sea Brent fields. It aims to stress the potential of the proposed approach for determining reservoir models respecting production data. The second one is also rooted in a real field case, but focuses on the matching of impedances.


INTRODUCTION
In many fields of geosciences like reservoir engineering or hydrology, due to the ever-increasing improvements in the price-performance ratio of computers, there is a growing demand for numerical models to represent geological formations. These models contribute to improve our understanding of underground fluid displacements. When considered as reliable, they serve as a basis for simulating and predicting fluid flows. For instance, reservoir models are used to forecast potential oil-and gas-recovery efficiencies, for well optimization, and for economic studies. Besides, aquifer models make it possible to predict the migration of pollutants or CO 2 plumes.
In a few words, a geological model is a three-dimensional grid with petrophysical properties such as porosity, permeability and facies assigned to every grid blocks. It often comprises several millions of grid blocks. The spatial distributions of petrophysical properties are very difficult to characterize for at least two reasons. First, petrophysical properties may be very heterogeneous (Dutton et al., 2003;Eaton, 2006). Second, there is a deep lack of quantitative data to calibrate models with realistic values. Data are usually split into two groups: static and dynamic data. Static data do not depend on time and are directly related to the property to be modeled. They include direct observations on well-exposed outcrops, well-log data, measurements performed on cores extracted from wells and baseline 3D seismic. Integrating static data into geological models usually calls for geostatistical techniques (Goovaerts, 1997;Chilès and Delfiner, 1999). On the other hand, dynamic data vary with time because of fluid displacements and are nonlinearly related to the target property. They encompass bottom hole pressures, flow rates, tracers... In the specific case of reservoir engineering, these data are tagged as production data. Since the late nineties, dynamic data also consist of 4D seismic data, that is, 3D seismic repeated over time. 4D seismic is a potentially powerful source of data for monitoring pressure, saturation and temperature changes related to fluid movements over large areas. Its main value is the additional information about reservoir connectivity, flow barriers or bypassed areas (Kawar et al., 2003). In the context of CO 2 sequestration, 4D seismic can help to quantify the amount of CO 2 injected and to detect leakages. The integration of dynamic data into geological models is an inverse problem solved using optimization techniques (Yeh, 1986;Tarantola, 1987;Sun, 1995). In reservoir engineering, it is also known as history-matching since it originally focused on the integration of production history (Nelson, 1960;Jacquard and Jaïn, 1965;Oliver et al., 2008). Mathematically, history-matching relies on the minimization of an objective function quantifying the least-square mismatch between the actual data and the corresponding simulated responses: the reservoir model is adjusted step by step until a reasonable data match is achieved. The final constrained model respects the required data, hence is more consistent. It can thus be used to forecast fluid displacements.
Because of the lack of observations and measurements, the constrained model is still uncertain. To make it more reliable and decrease the uncertainty in the resulting predictions, all available data must be integrated. The matching of production data only calls for quite simple workflows involving geological modeling, upscaling and flow simulation. As soon as 4D seismic is also considered, advanced workflows must be designed to incorporate both production and seismic data in a consistent way. The improvements in the quantitative processing of 4D seismic motivated many works about this topic (Landa and Horne, 1997;Pianelo et al., 2000;Gosselin et al., 2001;Waggoner et al., 2002;Stephen and McBeth, 2006;Roggero et al., 2007;. However, there is still no definite matching methodology and several issues remain to be addressed. 4D seismic can be accounted for under various forms: amplitudes (Arenas et al., 2001) or acoustic impedances derived from amplitudes in a preliminary step , impedances at base and monitor times or impedance differences between base and monitor times, impedances given in depth (Gosselin et al., 2001) or time (Fornel et al., 2007; domain. All of these may significantly impact the number of data, but also their informative content. In addition, the definition of the 4D seismic error term into the objective function is not straightforward. First, the weights balancing datasets as different as production and 4D seismic data must be properly estimated. Second, the least-square formulation used to measure the seismic error may be inappropriate to capture the relevant information (Aanonsen et al., 2003). Difficulties to achieve significant decrease in the seismic data mismatch were reported in several studies . Another issue is related to scales. Simulation workflows must handle several scales. Flow simulations are usually performed at the coarse scale while the petroelastic model, which provides the seismic responses, applies to the fine scale. Last, the widespread spatial coverage of seismic gives information about petrophysical properties in the inter-well regions. Capturing information in these areas calls for more flexibility in the model parameterization.
In this paper, we first present a general workflow for building reservoir models respecting both production and inverted 4D seismic data. Seismic data are referred to as inverted since we focus on the acoustic impedances or velocities derived from the inversion of the amplitudes of the seismic traces, and not on the amplitudes themselves. Basically, this workflow calls for geological modeling, upscaling, flow simulation, downscaling and petroelastic modeling. Then, we recap the main features of a few parameterization techniques developed especially for adjusting the spatial distribution of petrophysical properties within geological models while preserving their spatial variabilities. The following sections present two application cases. The first one is the "Brugge field" prepared for benchmarking purposes (Peters et al., 2009): the dataset to be matched encompasses 10 years of production history. The second case is also a synthetic case built from a real field. The reference dynamic data include both production and 4D seismic data.

METHODOLOGY
The joint integration of production and 4D seismic data into reservoir models involves an iterative optimization process since flow simulation is nonlinear. In practice, a forward simulation workflow is built from a sequence of activities: it eventually yields numerical responses that are compared to the data collected on the field. Typically, these numerical responses encompass the production responses provided at once by the flow simulator and 4D seismic responses derived from a petroelastic model.

Workflow
We develop a general workflow for incorporating both production history and inverted 4D seismic data (acoustic impedances, velocities or time shifts) into reservoir models. Its main steps are recapped hereafter (Fig. 1): -Geological modeling at the fine scale. As there is not enough data to accurately describe the distributions of porosity and permeability within the reservoir, we refer to a stochastic framework. Thus, petrophysical properties are considered as random variables whose spatial variabilities are characterized from the collected static data. A usual geological modeling process is recapped hereafter (Journel et al., 1998;Doligez et al., 2007). First, we generate a facies realization, which respects the observations at wells and the facies proportions estimated from logs and seismic. Algorithms often used to produce facies realizations are the Sequential Indicator Simulation method (Goovaerts, 1997), the Truncated Gaussian method or its extended version, the Pluri-Gaussian method (Chilès and Delfiner, 1999). Other newer techniques as yet in an immature stage can also be used: multi-point statistics (Strebelle, 2002;Salles et al., 2007), process modeling with cellular automata, etc. Porosity realizations are then drawn conditionally to the available porosity data to populate the facies realization applying for instance the Sequential Gaussian Simulation method (Goovaerts, 1997), the Turning Band method (Chilès and Delfiner, 1999) or the FFT-MA method (Le Ravalec-Dupin, 2005). At this point, permeability is either (1) computed from empirical relationships between porosity and permeability derived from data measured on cores or (2) simulated conditionally to the permeability data and the previously simulated porosity realization. The so-obtained model is called a geological model: it corresponds to a fine scale and usually encompasses a few millions grid blocks; -Performing a flow simulation at the fine scale is usually too CPU-time demanding. Therefore, the fine geological model is upscaled (Renard and Marsily, 1997). This process results in a coarse model, called reservoir model, which has to reproduce as closely as possible the flow behavior simulated for the fine geological model. This reservoir model is associated to a coarse scale and consists of about 100 000 grid blocks; -Flow simulation at the coarse scale. Flow equations are solved to simulate fluid displacements over time within the reservoir model (Aziz and Settari, 1979). The simulation outputs are production responses such as bottom hole pressures, flow rates or water cuts at wells. Besides, we obtain pressure, saturation and temperature grids describing the variations in these variables over the whole reservoir model at base and monitor times; -Downscaling of the simulated pressure, saturation and temperature grids. Seismic responses (velocities, impedances) are computed at the fine scale and depend on pressure, saturation and temperature grids at base and monitor times. Therefore, these grids must be first transferred from the coarse scale to the fine scale. This is the purpose of the downscaling step. It can involve a simple mapping, meaning that all fine grid blocks in a given coarse grid block are attributed the same pressure, saturation and temperature as the coarse grid block. Although this approach ensures mass conservation, it does not succeed in capturing the influence of fine-scale heterogeneities. To circumvent this drawback, a flow-based downscaling method was Flow chart describing the workflow followed to match both production and 4D seismic data.
proposed by Castro et al. (2006). However, its extension to any kind of flow patterns was not straightforward. Hereby, we employ a more flexible downscaling algorithm, which also calls for flow simulation (Enchery et al., 2007); -Petroelastic modeling (PEM). The PEM links the reservoir domain to the elastic one. It converts 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 seismic responses (P-and S-wave velocities, impedances). For simplicity, we neglect temperature variations and assume that the porous medium is isotropic and linearly elastic. The PEM considered involves a set of theoretical or empirical relationships. First, Hertz relations are applied to account for pressure effects: t is time and superscript ini stands for initial, that is time 0. K M and µ M are the drained bulk and shear moduli of the dry rock. P eff is the effective stress: it corresponds to the difference between confining pressure P c and pore pressure p. As fluids flow and pore pressure changes, the effective stress varies. h K and h µ are the Hertz coefficients with respect to the drained bulk and shear moduli. Assuming that the vertical stress is constant and that there is no lateral deformation (oedometric conditions), the above effective stress ratio is expressed as a function of pore pressure difference Δp: ν is the Poisson coefficient: it is a dimensionless coefficient varying between 0 and 0.5. Second, we refer to Gassmann equations (1951) to describe fluid substitution. They apply to homogeneous and isotropic rocks with pore space fully connected. These equations yield the low-frequency undrained bulk and shear moduli, denoted K and µ, respectively: K Gr and µ Gr are the bulk and shear moduli of the solid matrix. ϕ is the porosity and K F is the compressibility of the fluid saturating the pore space. Again, as fluids move into reservoirs, saturations, hence fluid compressibility, vary with time. Last, Pand S-wave velocities are derived from: where ρ is the density of the saturated rock. Impedance is the product of velocity and density. The PEM results in fine grids populated with velocities or impedances at the required base and monitor times. A low/high/band pass filtering step is added to ensure that the ultimate acoustic responses are in the same bandwidth as seismic data. Last, the acoustic velocity and impedance grids can be transferred to the time domain (Fornel et al., 2007;.

Inverse Problem and Parameterization
The difference between the actual dynamic data d obs and the corresponding simulated responses g(v) is measured in a least-square sense by the objective function: v is the set of unknown parameters and g is the operator mapping the parameter space to the dynamic data space. Covariance matrix C D quantifies the experimental and theoretical uncertainties. Various parameters of the workflow presented above can then be adjusted so as to minimize the objective function. A key one when dealing with fluid displacements in underground geological formations is the spatial distribution of petrophysical properties.
An "optimal" geological model is identified from the minimization of the objective function. However, this problem is known to be "ill-posed", meaning that there is no unique solution, that there may be no solution at all and that the identified model is highly sensitive to variations in the dynamic data. Several regularization approaches (Neuman, 1973) have been investigated to find solutions of the historymatching problems. One possibility consists of penalizing deviations from a prior model v prior , which is assumed to be known. This can be achieved by adding a term describing prior information into the objective function (Tarantola, 1987): Covariance matrix C V characterizes the uncertainties in the prior model. Another widely used approach is the Tikhonov regularization (Tychonoff and Arsenin, 1977). Then, the objective function is written: Regularization can also be performed by narrowing the parameter space. This alternative also calls for some knowledge about the properties of the searched solution. Examples of restrictions on the parameter space consist of setting upper and lower bounds to limit parameter variations or of assuming some kind of geological continuity. This motivated the development of geostatistics-based parameterization techniques such as the pilot point method (RamaRao et al., 1995), the gradual deformation method (Hu, 2000), the probability perturbation method (Caers, 2003) or the co-simulation perturbation method (Le Ravalec-Dupin and Da Veiga, 2011). These techniques rely on various geostatistical notions. The pilot point method refers to kriging to constrain realizations: the pilot point values are adjusted to minimize the objective function. The gradual deformation method is based upon the linear combination of Gaussian random functions with identical covariances. The simplest gradual deformation scheme consists of combining two independent sets of normal deviates z 1 and z 2 as: This relation ensures that z is also a set of normal deviates, whatever the value of the θ deformation parameter. Then, the value of the deformation parameter can be adjusted to reduce the data misfit. The normal deviates are used as a seed to draw a realization from a given random function, for instance to populate the reservoir model with porosity values. The continuous variations in the deformation parameter induces a continuous change in the resulting realization (Fig. 2). Therefore, the optimization process boils down to pick up the deformation parameter minimizing the objective function. Further details about the practice of the gradual deformation method in inverse modeling can be found in Hu (2000), Hu et al. (2001), Hu and Le Ravalec-Dupin (2004) and Hu and Jenni (2005). The probability perturbation method is very similar since it depends on the linear combinations of probabilities. Again, the combination weights are parameters that can be changed to reduce the data misfit. Besides, the co-simulation perturbation method involves correlation coefficients that can be adjusted during the optimization process. Therefore, these methods are used to modify the petrophysical properties populating the geological model from a reduced number of parameters while preserving spatial variability. They can be incorporated into matching workflows to minimize the objective function. Note that the objective function is restricted to the data mismatch term when applying geostatistical parameterization techniques. However, the convergence rates of such matching processes is often slow. Although the first iterations induce a significant decrease in the objective function, the following ones are less efficient. This is even worse when dealing with 4D seismic data. The matching process can be accelerated by applying the parameterization methods listed above to given sub-domains of the geological model. Nonetheless, this additional degree of freedom is not enough. The matching process can be further improved by refining the selection of the sub-domains to be modified. Early approaches focused on simple Voronoi polygons around wells (Hoffman and Caers, 2005;Gervais et al., 2007) or drainage areas between injectors and producers (Le Ravalec-Dupin and Fenwick, 2002;Hoffman and Caers, 2005). Then, focusing on the integration of 4D seismic data,  proposed to identify areas with strong seismic data mismatch and to define the target subdomains from the streamlines arriving to these areas. Another prospect to make the production and 4D seismic data matching process more efficient is the parameterization technique itself. Two options were investigated to significantly perturb petrophysical properties over the selected sub-domains. The first one relies on the local perturbation of the mean of continuous petrophysical properties such as porosity and permeability, while the second one focuses on the local perturbation of facies proportions. Variations in mean properties can be performed by referring to the pilot point method. In this special case, a secondary variable, that is the mean property, is also accounted for (Le Ravalec-Dupin, 2010; . Then, the mean values are adjusted so as to reduce the data mismatch and the resulting changes are mapped to the whole model applying co-kriging Gradual deformation process applied to a realization with a deformation parameter increasing from -1 to 1. instead of kriging. In addition, perturbing facies proportions can provide a lot of flexibility since facies strongly impact fluid flows. Ponsot-Jacquin et al. (2009) suggested to vary the ratio of the proportions of a group of facies, named "association", to the proportions of a larger group of facies, named "selection". The "association" is a sub-group of the "selection". The way to regroup facies depends on the case studied. For instance, the "selection" can correspond to the group of the most influential facies, that is the group of facies significantly impacting the objective function, and the "association" can include the facies with similar petrophysical properties. The ratio method makes it possible to vary the proportions of the facies included in the "association". These variations are then easily propagated to the other facies of the "selection" while the proportions of the facies not included in the "selection" are unchanged. The ratio method is very attractive because of its simplicity. However, when applied locally, it generates discontinuities at the boundaries of the modified sub-domains, which may be undesired. Then, a variant of it, referring to indicator cokriging, was introduced by Tillier et al. (2010) to preserve the continuity of facies proportions over the whole reservoir.

APPLICATION CASES
The complete matching workflow or only part of it was successfully applied to real field cases (Le Roggero et al., 2007;Pourpak et al., 2009; as well as laboratory cases (Soltani et al., 2010). Two new application cases are described below. Although synthetic, they remain realistic since they are inspired by real fields.
The first case illustrates the use of the workflow described in Figure 1 to match production data only. It outlines the needs for properly selecting the weights, which balance the contributions of the various data types accounted for: pressures, water and oil rates. In addition, it stresses the key role of parameterization. The second example focuses on the match of 4D seismic related data without explicitly accounting for production data. 4D seismic being mainly used for monitoring fluid front movements, this case aims to show that matching seismic data solely makes it possible to better capture the water breakthrough times observed at wells.

Brugge Case: Matching of Production Data
The Brugge test case is a synthetic oil field built for benchmarking purposes (Peters et al., 2009). We shortly describe the available dataset and explain the procedure developed to match the target production data.
The Brugge field is characterized by an E-W elongated half-dome structure. It is bounded by a large fault at its northern edge and has one modest throw internal fault. Its dimensions are approximately 10 × 3 km 2 . This field was inspired by typical North Sea Brent fields, i.e., reservoirs with delta plains and barriers. It consists of four main formations whose petrophysical properties are reported in Table 1. The Schelde and Maas formations have good reservoir properties unlike the two others. The Brugge reservoir is produced from twenty producers drilled close to the top of the dome and ten injectors all around (Fig. 3). A fine geological model with 20 million grid cells was initially built, with average cell dimensions of 50 × 50 × 0.25 m 3 . Then, this model was upscaled to a grid with 450 000 cells. The fine geological model was considered as the true model. Numerical well logs were extracted from it and were used to estimate the distributions of sedimentary facies, porosities, permeabilities, Net to Gross ratios, water saturations, gamma rays, bulk densities, etc. Once the well logs were generated, the true model was discarded. Then, the upscaled model was considered to be sufficiently reliable to compute production data consistent with the fine geological model. A flow simulation was run for the upscaled model to generate 10 years of production history. Noise was added to the resulting Bottom Hole Pressures (BHP) and oil rates to be more representative of a real case. Finally, the available data were the BHPs, oil rates plus all simulation input parameters (geometry, relative permeability curves, PVT data, well constraints, etc.). The only unknown parameters were the petrophysical properties within the whole model. Besides, several random geological realizations were generated from well logs from various stochastic simulation schemes. At last, Peters et al. (2009) delivered 104 realizations defined over a grid with 60 000 cells and consisting of nine layers (two in Schelde, three in Waal, three in Maas, and one in Schie). To sum up, the available data were the well logs, the first 10 years of production history and the 104 random realizations generated on a coarser grid.
Thus, we developed a workflow to determine porosity and permeability models that yield the best possible match of the 10-year production history. For simplicity, we assumed that there was a single facies. We started from an initial porosity model and applied the gradual deformation method (Hu, 2000) to modify it during the optimization process. This approach is appropriate within a two-order statistical framework, meaning that the mean and variogram of the target property must be specified. In the case studied, instead of focusing on well logs, we took advantage of the realizations with a single facies among the 104 given to compute the mean and variogram of porosity. We selected a porosity/permeability regression model whose coefficients were attributed values calibrated from well logs. All flow simulation parameters were set to the ones initially provided in the Brugge benchmark. Thus, a total number of 18 gradual deformation parameters (two for each of the 9 layers) was adjusted during the optimization process.
During a first optimization process, we defined the objective function as the sum of the least-square mismatch on oil rates, water rates and BHPs over all wells with equal weights. After nine loops, the optimization stopped. The evolution of the objective function during this first process is shown in Figure 4.
After nine loops, the objective function was reduced by more than 50%. Careful analysis of the results shows that most of this reduction was due to a better match of BHPs, while the contribution of rates was almost unchanged. To improve the rate match, we defined a new objective function by increasing the weights attributed to rates. Starting from the model obtained after the first nine loops, 11 more loops were carried out with the new objective function. Figure 5 shows that this second optimization process yielded an additional 25% reduction of the objective function. As the weights were modified, the value of the objective function changed: its value for loop 9 in Figure 4 is different from that of loop 0 in Figure 5.
The matched model obtained after the 20 loops reproduced the observations much better than the initial one, as illustrated in Figure 6 for the BHP at producer BR-P-5, in Figure 7 for the oil rate at producer BR-P-14, in Figure 8 for the oil rate at producer BR-P-18 and in Figure 9 for the water rate at producer BR-P-18. Focusing more specifically on well BR-P-5, we check that the overall trend is well captured in the long term: the pressure support due to water injection is properly reproduced. However, the match is not as good in the early times (< 1 000 days). This discrepancy can be explained by a poor modeling either of permeabilities locally around the  well or of interactions with neighboring wells. Improvements could be achieved by refining parameterization, for instance by performing local deformation around the target well.
Therefore, generally speaking, the final match is reasonable. However, when wells are considered independently, it is not as good as it could be and there is a need for varying the parameterization technique to add degrees of freedom to the optimization problem. As previously explained, gradual deformation parameters were used to modify entire layers. This deformation process should be localized around wells. In other words, it would be well-suited to split the layers into sub-domains and to assign deformation of parameters to each of them. This would induce an increase in the number of parameters to be accounted for, but would help to improve the data fit.

Synthetic Case: Matching of 4D Seismic Related Data
The second reservoir example studied is also inspired by a real field case. It encompasses 7 facies, all consisting of shales and sands. Their distributions are constrained by nonstationary proportions derived from the baseline seismic data acquired before production started. A specific feature is the very strong contrast between the petrophysical properties of shales and sands: permeability ranges from 0 to 5 000 mD Objective function versus chain number. Figure 6 Observations (red), initial model (dotted black) and matched model (plain blue) for the BHP at producer BR-P-5.

Figure 5
Objective function versus chain number. and porosity from 5% to 29%. The spatial distribution of facies is characterized by an anisotropic spherical variogram. The main correlation lengths are 1 320 m and 660 m along the horizontal diagonal axes and 5 m along the vertical axis. For simplicity, the geological model and the reservoir model are the same meaning that there is a single model. The grid includes 47 × 50 × 45 grid blocks of dimensions 33 × 33 × 1 m 3 . Production is driven by a producer and a water injector opened 6 months after the beginning of production (Fig. 10).
We first generated the reference reservoir model, considered as the true one. A facies realization was simulated on the basis of the pluri-Gaussian method (Doligez et al., 2007).
Then, the so-obtained facies model was populated with porosity realizations. Last, relationships between porosity and permeability were applied to assign permeability values to all grid blocks. At this point, we applied the workflow depicted in Figure 1 (minus the upscaling step which is useless here since the CPU-time for a flow simulation is low) to get reference production (flow rates, water cuts) and seismic responses (acoustic P-impedances at base and monitor times, monitor time being set to one year after beginning of production). We focused on a vertical cross-section including the injector and the producer (Fig. 10). P-impedance variations indicate fluid and pressure changes in the reservoir (Fig. 11).  Observations (red), initial model (dotted black) and matched model (plain blue) for the oil rate at producer BR-P-14.

Figure 8
Observations (red), initial model (dotted black) and matched model (plain blue) for the oil rate at producer BR-P-18.
A negative variation corresponds to a gas saturation and/or pressure increase while a positive variation suggests a water saturation increase and/or a pressure decrease. For the P-impedance sections shown in Figure 12, we observe that the positive (blue) patch on the left top side (close to the injector) is due to the water saturation increase. There is also a pressure increase, but its influence is overbalanced by the water saturation increase since impedances increase. Just below, there is a second blue patch to be related to the pressure decrease in this area. A red patch appears on the right, close to the producer: it corresponds to a gas saturation increase. In the middle, the alternating blue and red stripes result from competitive effects between pressure and gas saturation.
For our tests, we assumed that the spatial distribution of facies was unknown while all other parameters were known. Then, we defined an objective function with the reference P-impedance data at times T = 0 and T = 366 days. The purpose of this study was to investigate whether we could capture the reference water breakthrough at 900 days by matching only the impedances associated to early times. As above, we referred to the gradual deformation method to adjust the facies model and minimize the objective function. The 45 layers of the reservoir models were stacked into 5 units with equivalent thicknesses. Then, we considered two gradual deformation parameters per unit to modify the whole facies model from 10 parameters. Gradual deformation based matching processes rely on a sequence of optimizations. In the case studied, we ran successively 10 optimization processes, which resulted in an objective function decrease of 66% (Fig. 13). Observations (red), initial model (dotted black) and matched model (plain blue) for the water rate at producer BR-P-18.

Figure 10
Vertical cross-section between injector and producer.

Figure 11
Interpretation of the variations in V p /V s and impedance ratio against saturation and pressure changes.   Objective function against the number of flow simulations.
The P-impedance differences obtained for the reference, the initial and the final "optimal" facies models are compared in Figure 14. The reference impedance mismatch was clearly reduced close to the injector (left), but there is no improvement close to the producer (right). As stated above, we wanted to investigate how the match of impedances impacts production responses, more specifically, water breakthrough. The water breakthrough time for the initial facies model was 1 300 days. For the final "optimal" model, it was reduced to 1 080 days, which is close to the reference breakthrough time of 900 days. Therefore, although the matching of seismic-related data was not perfect, it contributed to make the model more reliable, hence more predictive.

CONCLUSION
In this paper, we developed an integrated workflow for reconciling reservoir models with both production and inverted 4D seismic data. It involves a sequence of simulation activities among which geological modeling, upscaling, flow simulation, downscaling and petroelastic modeling. Such simulation workflows yield numerical responses, which are then compared to the data to be matched. The data misfit can be reduced by varying any parameter involved in the workflow. The advantage of workflows is that any parameter change can be consistently propagated through the whole modeling process. Then, we recapped the main outlines of a few geostatisticsbased parameterization techniques, which make it possible to vary fine geological models from a limited number of parameters while preserving their spatial variabilities. Driving the required changes from these techniques at the first geological modeling level of integrated workflows makes it possible to get reservoir models consistent not only with production and inverted 4D seismic data, but also with static data. This integrated modeling process was applied to two test cases. The first one focused on the matching of production data only and the second one to the matching of 4D seismic inverted data. The integration of seismic data into reservoir models was shown to significantly enhance their reliability.