Reservoir Parameters Quantification from Seismic Monitoring Integrating Geomechanics

Résumé — Quantification des paramètres réservoir à partir des mesures du monitoring sismique, intégrant l’aspect géomécanique — Les mesures sismiques répétées dans le temps peuvent aider au suivi des mouvements de fluide au cours de la production. Le monitoring sismique peut potentiellement améliorer la récupération et optimiser les schémas d’exploitation des champs existants et des nouveaux champs. La production du réservoir induit des changements en saturation, en pression et en contraintes qui peuvent influencer le processus de propagation des ondes dans la roche. L’influence des contraintes effectives moyennes, dues au changements en contraintes totales moyennes et/ou en pression, sur les propriétés élastiques du réservoir n’est pas toujours prise en compte de manière appropriée. La méthodologie proposée évalue la contribution d’une approche géomécanique sur le calcul de paramètres élastiques. L’implémentation de cette méthode est composée de plusieurs étapes. Tout d’abord, pour prendre Abstract — Reservoir — Seismic measurements acquired at different stages in the life of a reservoir can monitor the fluid distribution over production time. Seismic monitoring has the potential to significantly enhance recovery and optimize exploitation schemes in existing and new fields. The production from the reservoir induces changes in saturation, pore pressure and stresses, which may influence the process of wave propagation in rocks. The influence of mean effective stresses, due to changes in mean total stresses or/and pore pressure on reservoir elastic properties are not always taken into account in a proper manner. The proposed methodology evaluates what is the contribution of a geomechanical approach on the computation of elastic parameters. The implementation of this method is composed of several steps. First, in order to take into account multiphase fluid flow, pressure and saturation are computed through reservoir simulation. Then, the computed pressure is used as a load in the geomechanical modeling. Using poroelasticity theory introduced by Biot, the contribution of rock deformation to fluid flow is evaluated. This allows the simulation of the stresses and pressure distributions throughout the whole geological structure over production time. The following step consists in generating elastic parameters as function of reservoir effects using rock physics. In order to take into account the mean effective stresses on seismic velocities, the Hertz-Mindlin’s contact theory is used. The saturation effect on seismic velocities is then considered by Biot-Gassmann’s relation. This methodology has been validated on real repeated data for monitoring an underground gas storage. This integrated procedure is then applied to several scenarii of production. Thus, the sensitivity of elastic parameters has been analyzed in order to differentiate the different reservoir effects. Then, using elastic parameters, seismic modeling allows the generation of seismic responses at well location that reveal the patterns of expected seismic monitoring results. Some seismic attributes, such as time-shift delay have been measured. They show that careful processing of seismic data is required for seismic monitoring. This integrated procedure applied on real data for monitoring an underground gas storage leads to demonstrate the play of differentiated attributes involving P and S -waves to separate mean effective stresses effects from saturation ones. Using elastic modeling, the impact of offset changes was demonstrated to measure reliable time-lapse parameters like time-shift delay attribute and


INTRODUCTION
Seismic monitoring (time-lapse or 4D) has emerged as one of the most important technical developments in the oil and gas industry for this decade.This technique has the potential to significantly enhance recovery and optimize exploitation schemes in existing and new fields.It aims at monitoringby repeated well seismic, 2D or 3D seismic surveysseismic changes, velocity and density, related to fluid, stresses and temperature changes during the production of a field.
Changes in density and velocity result in impedances changes which, under favorable conditions, can be detected in seismic data.The probability of success of the technique 556 Abstract -Reservoir Parameters Quantification from Seismic Monitoring Integrating Geomechanics -Seismic measurements acquired at different stages in the life of a reservoir can monitor the fluid distribution over production time.Seismic monitoring has the potential to significantly enhance recovery and optimize exploitation schemes in existing and new fields.The production from the reservoir induces changes in saturation, pore pressure and stresses, which may influence the process of wave propagation in rocks.The influence of mean effective stresses, due to changes in mean total stresses or/and pore pressure on reservoir elastic properties are not always taken into account in a proper manner.The proposed methodology evaluates what is the contribution of a geomechanical approach on the computation of elastic parameters.The implementation of this method is composed of several steps.First, in order to take into account multiphase fluid flow, pressure and saturation are computed through reservoir simulation.Then, the computed pressure is used as a load in the geomechanical modeling.Using poroelasticity theory introduced by Biot, the contribution of rock deformation to fluid flow is evaluated.This allows the simulation of the stresses and pressure distributions throughout the whole geological structure over production time.The following step consists in generating elastic parameters as function of reservoir effects using rock physics.In order to take into account the mean effective stresses on seismic velocities, the Hertz-Mindlin's contact theory is used.The saturation effect on seismic velocities is then considered by Biot-Gassmann's relation.This methodology has been validated on real repeated data for monitoring an underground gas storage.This integrated procedure is then applied to several scenarii of production.Thus, the sensitivity of elastic parameters has been analyzed in order to differentiate the different reservoir effects.Then, using elastic parameters, seismic modeling allows the generation of seismic responses at well location that reveal the patterns of expected seismic monitoring results.Some seismic attributes, such as time-shift delay have been measured.They show that careful processing of seismic data is required for seismic monitoring.This integrated procedure applied on real data for monitoring an underground gas storage leads to demonstrate the play of differentiated attributes involving P and S-waves to separate mean effective stresses effects from saturation ones.Using elastic modeling, the impact of offset changes was demonstrated to measure reliable time-lapse parameters like time-shift delay attribute and amplitude variations.
heavily depends on many factors, like reservoir parameters (depth, rock and fluid properties, pressure, etc.), nature of the recovery processes and the repeatability of the different seismic surveys (Wang, 1997).A careful feasibility study on the field of interest, coupled with a clear reservoir objective, is required to give a realistic estimate of what seismic monitoring can provide for reservoir management.
The Céré-la-Ronde underground gas storage reservoir, in the Paris Basin, a test site to study and enhance reservoir seismic monitoring method, is a water-bearing sandstone reservoir in a faulted anticline structure (Fig. 1).The depth of reservoir is about 900 m and its thickness is about 25 m, at the structure top (Fig. 2).
Regarding seismic monitoring, the site presents favorable characteristics-shallow reservoir and gas injection process (Wang, 1997)-, and unfavorable ones-low repeatability due to unpredictable statics.The different repetitive seismic data (sonic logs, well seismic, walk-away and 2D seismic) acquired so far have generated a qualitative seismic interpretation of the location of the gas bubble, by studying fluid saturation (Meunier and Huguet, 1998).However, between 1994 and 1997, two sonic logs showed subtle differences in P velocity (Vp) not explained solely by saturation variations.Changes in pore pressure and stresses also influence the reservoir elastic properties.Hence, we used geomechanical modeling, to evaluate quantitatively, how exploiting the gas reservoir impacts the seismic measurements.
Several previous studies have shown the interest of using both fluid saturation and pore pressure to interpret time-lapse seismic data (Lumley et al., 1997;Johnston et al., 1999).
Often, when the influence of mean effective stresses on seismic velocities is studied, only pore pressure is considered and mean total stresses are seen as constant.As stresses evolve during reservoir exploitation, we will consider variations in the mean effective stresses, induced by changes in both the mean total stresses and pore pressure distribution, using geomechanical modeling.
This paper presents a methodology developed in order to integrate geomechanical modeling in the computation of seismic velocities.It implies the combination of geomechanics with geophysics.The implementation of this method is composed of several steps.Our method begins by computing, in a reservoir simulation, pore pressure and saturations.Pore pressure is a key input in the geomechanical modeling that produces mean effective stresses.These and the saturations are used to update seismic velocities in accordance with rock physics theory.In the final step, the introduction of a wavelet allows seismic modeling and the study of seismic attributes.The workflow of the methodology is summarized in Figure 3.

GEOMECHANICAL MODELING
The methodology aims at evaluating what is the contribution of a geomechanical modeling, assuming linear poroelasticity, on the computation of seismic velocities.This approach implies a forward modeling first involving reservoir simulation.To improve fluid flow description, pressure computed by reservoir simulation, is used as a load in the geomechanical modeling.

Reservoir Simulation
Reservoir simulation describes fluid mass balance (changes in pressure, temperature and multiphase pore fluid saturation) due to reservoir exploitation.Darcy's law can be substituted into the fluid flow mass balance, to give a diffusivity equation, which in monophasic case, expresses as: (1) where C t = C fl + C p is the total compressibility.

Geomechanical Modeling Principles
The contribution of rock deformation to fluid flow is studied using the poromechanical theory introduced by Biot (Biot, 1941).In his description, the porous medium is considered as the superposition of two phases: a solid one (the rock) and a fluid one (Boutéca, 1992).Hence, any variation of fluid mass is due to either rock deformation or fluid flow.Both the deformations and the variation of the fluid mass (m) characterize porous medium transformations.The state φµ Ġeology

Linear Poroelasticity
The poroelasticity model assumes infinitesimal strains and fluid mass variations.In these conditions, the existing relation between: is assumed to be linear.Then, the constitutive law is written: (2) C is the 4 th order stiffness elastic tensor for a dry anisotropic solid.M is the Biot modulus given by: (3) where K fl and K mat are respectively the bulk modulus of fluid and of rock matrix.In the isotropic case, B is reduced to α, the Biot's coefficient, and Equation (2) becomes: (4) (5) and the Biot's coefficient (α) is given by: (6) Combining Equations ( 4) and ( 5) gives: (7) We note that the material constitutive law (7) can be expressed using Biot effective stress, defined as: (8)

Fluid Flow Equation
Using Darcy's law, fluid flow mass balance and the second state Equations of poroelasticity given by Equation ( 5), we obtain the coupled diffusivity equation: (9) Equations ( 7) and ( 9) indicate coupling between strains and pressure.The constitutive laws given by Equation ( 4) and diffusivity equation allow using numerical simulation to evaluate mean effective stresses in a geological structure.

Geomechanical Behavior Modeling of Céré-la-Ronde
For the geomechanical modeling, a finite element code has been used.The reservoir simulation model has been used and embedded with overburden, underburden, for the geomechanical modeling.To build the geomechanical model, we need to define the mesh, the geomechanical properties per layer, the initial state, the load and the boundary conditions (Vidal et al., 2000).Well A

Figure 4a
Reservoir simulation modeling grid.

Figure 4b
Geomechanical modeling grid with corner cells.

Mesh
For the geomechanical modeling, the model has to cover the whole geological description, from the surface to a depth of 1500 m (Fig. 4b).To build the reservoir part of the geomechanical model, the reservoir model in Figure 4a has been used and the shape of cell has been slightly adapted in order to get corner cells and continuous structure, as shown in Figure 4b.To complete the definition of the geomechanical model, overburden and underburden layers were added, according to the geological layering.
The whole model is 7000 m-long, 1500 m-depth, has 37 cells in the X-direction and 27 cells in the Z-direction.The size of the 2D-model is 37*27 = 999 active elements, with 3126 nodes (high order elements).

Geomechanical Properties
The geomechanical properties needed for the modeling are described in Table 1.In poromechanical models, properties are defined statically.For reservoir and caprock layers, geomechanical propertiers are extracted from core measurements (when available) but for the other layers, poromechanical parameters could only be derived from sonic logs.We must note that this approach can only give a trend for mechanical properties because from Vp and Vs, we access only undrained dynamic parameters.The change from dynamic to static parameters is still empirical and we have used also, characteristics of analogs.

Initial State
Initially, vertical stresses (σ v ) are determined by rock densities and horizontal stresses (σ h ) by an estimated stress ratio (σ h /σ v ).The initial pore pressure has been fixed in the structure: for reservoir levels, pore pressure is given by reservoir simulation and for other layers the gradient is hydrostatic.

Loading
The reservoir exploitation, gas in and gas out, modifies the geomechanical equilibrium of the structure.The load for the geomechanical modeling is determined from cycling variations of the pore pressure in the reservoir (computed by the reservoir simulation).The Figure 5 shows over time, the cycling variations of pressure in the reservoir at well location.

Boundary Conditions
The lateral extension of the 2D-plane geomechanical model corresponds to one of the reservoir.For geomechanical considerations, it can be important to model a larger extension.In this study, in order to keep a simple geomechanical model, we do not explicitly consider the sideburden.Several sensitivity studies allow us to quantify the influence of this assumption.Some available pressure measurements in the upper aquifer have been used in order to validate which kind of boundaries allows such scaling.The variation of pressure measured in the upper aquifer is best explained by the boundary conditions given in Figure 6.For this set of boundary conditions, there is no lateral displacements at external boundaries.Then, the geomechanical modeling has been performed by a finite element numerical code to compute pressure and mean stresses in any layers.

Results
The results of the geomechanical modeling (i.e., pore pressure and mean total stresses variations) have been plotted over depth at the well A, for different production times (Figs.7a and 7b).In Figure 7a, pressure in the reservoir level (R1-R2-R3) is taken from reservoir simulation and used as a input for the geomechanical modeling.Over the reservoir R1, pore pressure is computed by the hydromechanical code in the overburden layers.Pressure in R1 varies from 7.8 MPa to 11.4 MPa.In Figure 7b, mean total stresses have been computed by the hydromechanical code.We can see that mean total stresses are not constant and vary from 16.8 MPa to 18.2 MPa.The mean effective stresses (σ eff ) is defined as difference between total mean stresses and pressure.
It is useful to note that this geomechanical modeling performed here is site dependent.The obtained mean effective stresses depend on the nature of the reservoir rocks, the reservoir depth, the nature of the overburden rocks, the pressure distribution, the regional stresses state, etc.

Figure 7a
Computed pore pressure at well A for different production times.

Figure 7b
Computed mean total stresses at well A for different production times.
contact models based on Hertz-Mindlin contact theory (Mindlin, 1949).In this technique, two identical spherical grains of radius R are deformed by normal and tangential forces.The radius of the contact area is a function of the mean effective stress (Fig. 8).Therefore, effective shear and bulk moduli are also linked to the mean effective stresses.Using both moduli, Vp and Vs velocities may be computed.Hertz-Mindlin theory assumes that velocity varies with σ eff raised to the 1/6 th power: (10) Some laboratory measurements on samples gave a smaller exponent: 0.09 for Vp and 0.13 for Vs (Fig. 9).
If the initial velocity (V 0 ) is known, the new velocity (V 1 ) can be computed as function of initial mean effective stresses (σ eff0 ) and updated mean effective stresses (σ eff1 ): (11) The velocities P and S are linked to K u and µ by the following equations: (12)

Impact on Seismic Velocities
To model the impact of geomechanical approach, we have considered a reservoir layer with a constant velocity, V 0 = 3100m/s and with a thickness of about 30 m.We have compared the difference for time delay attribute (DT) between two approaches: the geomechanical one and the nongeomechanical (Fig. 10).

Geomechanical Approach
In this approach, both pressure and mean total stresses are varying: they are modeled by the geomechanical modeling.The computed DT is 235 µs (one way travel time) for a 3.6 MPa pressure variation.

Nongeomechanical Approach
In this approach, only pore pressure is varying and the mean total stresses are considered as constant.The computed DT is 392 µs (one way travel time) for a 3.6 MPa pressure variation.
Not considering a geomechanical approach (mean total stresses considered constant) as commonly done, might lead to a time delay attribute error of about 75% in comparison with the geomechanical approach where the mean total

Figure 9a
P velocities variations versus effective pressure, for three rock samples, from wells G, A and B (Zinszner, 1997).

Figure 9b
S velocities variations versus effective pressure, for three rock samples, from wells G, A and B (Zinszner, 1997).
stresses are varying.Considering only pore pressure effect (without mean total stresses effect) may lead to an error on saturation influence evaluation.
For the geomechanical approach, the computed DT is 235 µs (one way travel time), and 470 µs (two way travel time), for the greatest pore pressure variation (3.6 MPa pressure variation).The average 2D time-shift resolution is about 250 µs.This preliminary result from 1D modeling allows us to conclude that the pore pressure variations are too weak to be detected on 2D seismic.However, a well seismic with a better seismic resolution may allow us to detect pressure effect.

Impact of Gas Saturation on Seismic Velocities
Rock physics models also help evaluate the effect of gas saturation on effective bulk modulus.If the saturation state is known, the undrained bulk modulus can be used in the Biot-Gassmann's Equation (Gassmann, 1951) to calculate the drained modulus: (13) The most common in situ problem is to predict the change from one fluid to another.One procedure is simply to apply the Biot-Gassmann equation twice.So, after applying mean effective stress effect on velocities, we transform the resulting undrained bulk modulus from the initial fluid saturation to the new fluid-saturated state.In Equation ( 13), K mat may be estimated by effective media model, knowing mineral composition (Mavko et al., 1998).The effective bulk modulus of pore fluid, K fl , is evaluated by Reuss average knowing water and gas saturation (Mavko et al., 1998).Gas density and bulk modulus are determined under reservoir conditions (Batzle and Wang, 1992).Knowing saturation and porosity, dynamic elastic moduli (Vp, Vs) can be updated using Equations ( 12) (Fig. 11).

Rock Physics Modeling Procedure
The workflow of the rock physics modeling is illustrated in Figure 12.
Figure 12 The rock physics modeling scheme of the two reservoir effects : mean effective stresses and gas saturation on seismic velocities.

Validation on DSI (Dipole Sonic Lmaging) Data.
Figure 2 shows the 1993 Vp log on which the methodology was applied (Fig. 2).This (1993) is the initial state-no gas in the reservoir.Figure 13 shows P velocity in 1994 and 1997, and the results of rock physics modeling.In

Gas saturation variations σ eff variations
Hertz-Mindlin's model Biot-Gassmann's Equation applied twice Comparison of DT evolution for a geomechanical approach and for a nongeomechanical approach.With V 0 = 3100m/s.Figure 13e, only the fluid substitution effect has been applied to the 1993 data.The 1994 and 1997 logs (Fig. 13d) show, in the sandstone interval (923 m-930 m), a subtle difference not explained by the saturation effect only.In Figure 13f, where σ eff and the saturation effects have been both applied to the 1993 data, we observe a difference between the curves similar to that in Figure 13d.Thus, real data validate the model obtained by combining both effects on seismic velocities.

Elastic Parameters Modeling
These results prompted us to use the methodolgy on three values of mean effective stresses (minimum, mean and maximum) and three saturation states (Sw=100%, gas height = 8 m and 16 m) (Table 2).These nine production possibilities will help evaluate the impact of both pore pressure and saturation in a quantitative way (Vidal et al., 2001).
The computed elastic parameters are P and S impedances (Ip and Is), Poisson's Ratio (PR), µρ, λρ, λ/µ (Lame's coefficent).These attributes were crossplotted to identify which (if any) define reservoir properties.Figure 14, some combinations of time-lapse attributes involving shear wave or shear modulus, shows that discriminating different reservoir effects is possible (i.e., between the gas saturation, the σ eff , and the two combined effects).

Convolutional Seismic Modeling
We next used 1D convolutional synthetic modeling to obtain synthetic traces for each case.This step convolves the reflectivity series with a wavelet to generate seismic amplitudes.A Ricker wavelet of 80 Hz central frequency was used.Figure 15 shows the traces for the nine possibilities.The effect of σ eff on seismic amplitude is obviously weak but the addition of gas seems to change amplitude.Nevertheless, we estimate a time-shift attribute (DT) induced by change in σ eff to be larger than that induced by gas saturation (Table 3).
Thus, even when the change in mean effective stresses effect has only minor impact on amplitude, it is apparent on the time-shift attribute.(1994 and 1997).

Elastic Seismic Modeling
We mentioned previously that involving S-wave is really useful to discriminate the different reservoir effects, so we carried out elastic modeling.Figure 16 demonstrates by showing difference sections between monitor data and base data.As expected, σ eff (Fig. 16a) has little impact on amplitude and the addition of gas (Fig. 16b) gives larger changes.A time-shift attribute (DT) between the different surveys is often computed on stacked data in seismic monitoring studies.DT measurement is based on the correlation function between the reference seismic trace (before gas injection) and a repeated seismic trace (after gas injection or after pressure increase).We measured it for stacked data and also prestack data to analyze DT versus offset (Fig. 17).In our study, DT measured with near traces is better than with all traces because the former equals to the theoretical one computed from 1D convolutional synthetic seismic trace.With acoustic modeling, a constant value of DT was observed.That shows that variation of DT with offset on elastic modeling is due to converted waves.Our reservoir is thin, and DT is below sampling rate (1ms) so we need to be as accurate as possible.Thus we recommend high-precision processing and that DT be measured on partial stack with near traces.

DISCUSSION AND CONCLUSION
A methodology has been developed to combine geomechanics and geophysics for reservoir seismic monitoring.The results from the geomechanical modeling were then used for rock physics modeling.The results from this second modeling have been interpreted in terms of P velocity and time-shift attribute (DT).We concluded that failure to consider a geomechanical approach may cause important errors in DT computation.Modeled Vp logs have been compared to sonic data acquired at different time during production.The final result-modeled Vp log with Hertz-Mindlin's contact model and Gassmann's fluid substitutionagrees with real data and consequently validates the method.Many elastic parameters were estimated using rock physics modeling.This modeling using synthetic logs computation based on real data, demonstrated the play of shear waves or shear modulus data to differentiate σ eff effects from saturation ones.We also work with synthetic seismic data to study direct seismic attributes especially time-shift delay.Because differential attributes, which involve S-waves, are the most efficient for discriminating reservoir properties, elastic seismic modeling is necessary to study both saturation and σ eff effects.Our case study showed that the time-shift attribute induced by P-wave propagation needs careful processing and proper method to be measurable.It should be computed on partially stacked data with near offset traces.1D convolutionnal modeling for each production case with an 80-Hz Ricker wavelet.
Figure 17 Time-shift attribute (DT) versus offset for σ eff effect and for gas saturation effect.For comparison, the value of DT computed for full-stack trace is 530 µs for σ eff effect and 790 µs for gas effect.
Figure 1Location of the Céré-la-Ronde gas storage site.
The reservoir simulation gives pressure and pore fluid saturation.The description of pressure is used as input for the geomechanical modeling.The reservoir simulation uses 3D finite volume code.The model covers 336*228 km 2 and uses about 12 000 active cells.The most refined level of the reservoir model (zone of interest) is split into three reservoir layers: R1, R2 and R3, separated by shaly layers: C2 and C3 (Fig.2).So far, only the reservoir R1 is used for gas storage.The size of this refinement level of the model (zone of interest) is 36*30*7 active cells.Each cell is 200*200 m in the X-Y-plane and varies between 5 m for R1 and 80 m for R3 in the Zdirection.The model matches transmissibility through faults, rock compressibilities, relative permeability and local/global barriers.The main observations used in the history matching are pressures, injected and withdrawn gas volumes, watercuts in the wells and case hole saturation logs.The Figure4apresents one section, issued from the 3D model, with 36 cells laterally and 7 cells vertically, which represents 7000 meterlong and 120 meter-height.
Figure 2Composite log display of the Céré-la-Ronde reservoir showing lithology log, porosity, Vp, Vs and density logs from 1993 data, before gas injection.

Figure 5 Cycling
Figure 5Cycling variations of the pore pressure at injector well, used as a key input for the geomechanical modeling.

Figure 11 P
Figure 11 P velocities variations as function of water saturation and porosity for Céré-la-Ronde case.
Figure 14Time-lapse µρ, time-lapse λρ, time-lapse Ip for three cases: σ eff effect, gas saturation effect and their combined effect.