Monitoring of SAGD Process: Seismic Interpretation of Ray+Born Synthetic 4D Data

Résumé — Monitoring de procédé SAGD : interprétation sismique de données 4D synthétiques ray+Born — L’objectif de cette étude est d’évaluer quelle information de production peut être déduite d’une campagne sismique 4D durant le procédé de récupération par injection de vapeur SAGD ( Steam-Assisted Gravity Drainage ). En plus des hétérogénéités réservoir d’origine géologique, de nombreux facteurs interagissent pendant la production thermique d’huile lourde et de bitume, ce qui complique l’interprétation des données sismiques 4D : variation de la viscosité de l’huile, des saturations en fluide, de la pression de pore, etc. Abstract – Monitoring of SAGD Process: Seismic Interpretation of Ray+Born Synthetic 4D Data – The objective of this study is to evaluate which production information can be deduced from a 4D seismic survey during the Steam-Assisted Gravity Drainage (SAGD) recovery process. Superimposed on


INTRODUCTION
The performance of heavy oil production by Steam-Assisted Gravity Drainage (SAGD) recovery process is affected by reservoir heterogeneities. Pressure and temperature variations during thermal production operations induce stress changes in the reservoir and in the surrounding media. These modifications of the stress state may imply deformations that can, in turn, have an impact on reservoir production. These changes have therefore an influence on both densities and seismic velocities, consequently on the wave propagation into rocks and fluids along production and 4D seismic data. On such heavy oil reservoirs produced by SAGD, the literature shows that associated 4D seismic data are analyzed either in a qualitative way [1,2], or in a more precise/quantitative way by using PetroElastic or Rock Physics Modelling [3][4][5]. Moreover, additional information of higher frequency content than traditional seismic has proven to be of great help for understanding the reservoir heterogeneities and production problems (e.g. [6], combining standard and crosswell seismic over 500 Hz).
For this study, a model representative of an Athabasca heavy oil unconsolidated sands reservoir was constructed [7]. First, the initial static model was constructed with a geostatistical approach, considering five lithofacies for the reservoir (three sandstones from coarse to fine-grained, and two shaly facies) and three lithofacies (from sands to shales) for the overburden up to the surface. Second, the thermal production of heavy oil was then simulated along calendar time with two coupled fluid-flow and geomechanical models. The use of a negative thermal expansion coefficient for the shales allowed the possible collapse of the shale materials and thus a potential pass through for the steam. Elastic models were then generated along calendar time: the elastic parameters were carefully estimated according to the oil viscosity (temperature and pressure dependent), the oil and water saturations and the geomechanical effect through a stress Hertz coefficient.
Four production stages are considered here: the initial stage before steam injection, and the stages of 1st month of production, 6th month of production, and 3rd year of production.
The 4D seismic modelling mimics a permanent surface acquisition survey. To focus on the reflections generated within the reservoir zone, a 3D target-oriented modelling tool is chosen. This tool is based on the ray+Born formalism [8][9][10] which allows us to compute the P-wave elastic response, by correctly handling the seismic amplitudes as a function of source-receiver offset. In particular, the effects of the overburden and of the survey are taken into account. Some real noise is then added to the computed seismic data.
Based on these 4D stacked synthetic seismic data, an interpretation is carried out to answer these points: -first, better understand the evolution of the 4D seismic amplitudes as a function of calendar time in terms of production (pressure, temperature and fluid saturations) and mechanical effects (stress and strain); -second, provide guidelines for the interpretation of real 4D seismic data; reservoir heterogeneities of geological origin, many factors interact during thermal production of heavy oil and bitumen reservoirs, which complicate the interpretation of 4D seismic data: changes in oil viscosity, in fluid saturations, in pore pressure and so on. This study is based on the real Hangingstone field case of the McMurray formation in the Athabasca region (Canada). In previous works, an initial static model (geology, petroacoustic and geomechanical) has been constructed and a thermal production of heavy oil with two coupled fluid-flow and geomechanical models has been simulated. Seismic parameters (density, compression velocity and shear velocity) of the saturated rocks have then been computed from mechanical and reservoir parameters at several stages of the production. A repeated acquisition survey is modelled at different stages of SAGD production. This is performed using a 3D seismic modelling approach. To focus on the reflections generated within the reservoir zone, a target-oriented modelling is chosen. It is based on the ray+Born approach which permits to compute the P-wave elastic response by correctly handling the seismic amplitudes as a function of source-receiver offset. Real incoherent noise is added to the zero-phase synthetic data to produce a more realistic result. The noise-free and the noisy synthetic data are processed to get stacked and time migrated images. A simple processing workflow leads to image the steam chamber development, in particular its V-shape in radial section, and to observe time-lapse in the reservoir zone. An interpretation work is then carried out. Some seismic attributes like RMS values of amplitude changes between stages, energy, time differences of reservoir bottom between stages, etc. are computed from the synthetic (noise-free and noisy) seismic data. Some of these attributes prove to be robust to the noise and to show some production effect. Possible trends between these attributes and the modelled reservoir/geomechanical properties (lithofacies, pressure, temperature, steam saturation, etc.) are also evaluated. Finally, geobodies are extracted from the seismic attributes.
-third, assess the interest of permanent acquisition technology in heavy oil reservoirs produced by SAGD.
In this paper, we first briefly summarize these early works on model building and steam injection/production, then we present our results on 4D seismic modelling and interpretation related to some seismic attributes and their link with the reservoir/geomechanical properties (qualitative or quantitative), as well as geobody computations (see Fig. 1 for the general scope).

A REALISTIC MODEL OF A SAGD WELL PAIR [7]
The model represents the Hangingstone field located in Northern Alberta (Canada) in the Athabasca oil sands province. The reservoirs are made of unconsolidated sands from the fluvio-estuarine McMurray Formation (Mannville Group, Lower Cretaceous). Reservoirs from this area exhibit a complex internal architecture which is associated with heterogeneities (shales) specific to fluvial and estuarine environments. The heavy oil reservoir is at a relatively shallow depth of 260 m in average. Its thickness is close to 50 m.
The simulation of thermal production of heavy oil in the reservoir is fully described in [7]. The approach is based on a fluid-flow model explicitly coupled with a geomechanical one.
The main steps of the simulation workflow are summarized below.

Geological Model
A 3D detailed geological model of the reservoir has been constructed on a very fine grid by a geostatistical approach, constrained by both horizontal and vertical wells (Fig. 2): -first, the geological model was defined at a very fine scale near the well bores in order to preserve the description of heterogeneities (shales); -second, it was populated with lithofacies (see Tab. 1) and two initial petrophysical properties (porosity and permeability). The overburden up to the surface was modelled on a much coarser grid in the vertical direction (with sandstones, shales and sands, shales lithofacies).
Finally, a SAGD well pair was extracted from the whole model. General scope of the seismic modelling and interpretation study.

Reservoir/Geomechanical Model
This fine-scale geological model was then used to build a mechanical model in the reservoir part. Poro-elastic properties were assigned to the cells according to the lithofacies and associated porosity. This assignment was based on a consistent interpretation and integration of available log data and core data. At a temperature of 10°C, the oil viscosity is about 2 000 000 cP and its density is about 8°API. Production data available for this study consisted of steam injection and oil production rates. A simulation was performed for a selection of time periods of steam injection on the SAGD well pair, with a fluid-flow model coupled with a geomechanical one. The modelling focused on the first six months of production and on a later period comprised between one and three years. An explicit coupling scheme was performed with an update of the reservoir permeabilities at chosen steps of simulation. The reservoir simulation results are in good agreement with field data. Figure 3 presents the simulated extension of the steam chamber after six months and three years of production. The 3D envelope corresponds to a minimum temperature of 100°C, which means that every cell inside the envelope has a temperature ranging between 100°C-280°C (maximum temperature observed in the reservoir during production). This figure shows a high degree of heterogeneity of temperature distribution along the well pair, especially for short durations. Indeed, after six months the steam chamber development seems to be confined vertically from sections 1 to 11. These sections coincide with the presence of heterogeneities (shales) at the heel of the well (see Fig. 2).

4D Elastic Model
Because there are three main sources of stress dependency of wave velocities (changes in porosity with stress, existence of Distribution of heterogeneities in the reservoir (top) and on selected sections and along SAGD well pair (bottom). See Table 1 for lithofacies description. Reservoir cube of 100 m wide (X), 820 m long (Y), 50 m high (Z).
grain contacts and presence of cracks) a stress-sensitive rockphysics model was required. The Hertz-Mindlin's contact theory was used as a first approach. This model is based on the evolution of the contact surface between two spherical particles and accounts for the impact of mean effective stresses on both P-and S-wave velocities. As the impact of temperature on these velocities is not well established, a simple model issued from literature was used. The effects of fluid saturations on effective bulk modulus, then on velocities, were inferred from Biot-Gassmann theory. The initial seismic parameters (density, compression velocity and shear velocity) of the saturated rocks were estimated from initial mechanical and reservoir parameters in coherence with the well log data. Before the warm-up phase, only water and oil components are present in our fluid model, but heavy oil has been considered almost a solid due to its high viscosity. Therefore the generalized Gassmann's equations [11] were handled under the simplifications of a homogeneous grain rock and of a uniform mascroscopic fluids repartition.
Another challenge was that the determination of seismic parameters requires dynamic mechanical parameters (Young's modulus, Poisson's ratio), whereas mechanical parameters used in mechanical modelling are static.
Each material was assigned the elastic properties corresponding to the initial water saturation set to 0.15 in the reservoir sands materials, and set to 1 in all other materials, as observed from well logs. The 3D full model was then populated with density, P-and S-velocities (Fig. 4). In the overburden, the P-to S-velocities ratio is in a 2.5 to 3.1 range, whereas in the reservoir it falls between 2.0 and 3.0.
Concerning the 4D model, the basic idea was that the seismic properties of the oil saturated rocks are related to seismic properties of heavy oils which depend on density, composition, temperature and gas/oil ratio.
The evolution of physical parameters during steam injection was computed thanks to the coupling of the reservoir and geomechanical models. These changes were used  Steam chamber extent corresponding to temperatures over 100°C after six months (left) and three years (right) of production. to dynamically update the seismic velocities both in the overburden and in the reservoir: -temperature, saturation and pore pressure evolutions, computed by the fluid flow simulation, were taken into account in the reservoir only; -stress and strain evolutions, computed by the geomechanical simulation, were taken into account in both the reservoir and the overburden regions. The evolution of these physical parameters drives the modification of mechanical parameters (dry modulus, fluid modulus and density) which pilot saturated modulus. The evolution of both saturated modulus and density pilots wave velocity changes. P-and S-velocities decrease inside the slim core of the steam chamber associated with the presence of gaseous water and in heated-up shales (elsewhere it depends). Our objective is to compute the wave reflections occurring in an elastic target zone. We consider a 3D seismic modelling based on the ray+Born formalism [12] which corresponds to a linearization of the wave equation around a reference medium. On the one hand, the wave propagation between the surface and the target is simulated with the help of a 3D ray-tracing in the reference medium. On the other hand, the modelling of the reflections/diffractions in the target is performed with the help of amplitude computing using the Born approximation [13]. The ray+Born technique allows us to compute the elastic response, by correctly handling the seismic amplitudes as a function of source-receiver offset. In particular, the effect of the overburden and survey is included. The multiple reflections and refractions are not modelled. The main advantage of this method is that the simulation is fast computing; and the main disadvantage is that the modelling is not suitable for models with strong velocity variations.

4D SYNTHETIC PERMANENT
The Born approximation decomposes the model m into the sum of a reference medium m 0 and a perturbation target δm. The reference medium m 0 corresponds to the large spatial wavelengths of the model and governs the wave propagation in terms of traveltime, geometrical spreading and wave front directivity. We call Green's function the seismic response of a unitary point computed in the reference medium for an impulse source. The perturbation target δm corresponds to the small spatial wavelengths of the model and influences the wave reflection and diffraction amplitudes. The Born approximation is valid for relative perturbations of low amplitude and of small size compared with the wavelength of the wavefield. More precisely, the validity condition may be written as [14]: where λ is the wavelength of the wavefield, R is the significant size of perturbation elements and δp/p 0 is the average value of perturbation parameters. In the case of a very small significant size of perturbation elements, the average value of perturbation parameters can be not very weak such as a value of 0.1. This statement has been confirmed by many studies on different types of elastic model [15,16].
We used an in-house software for computing the P-wave elastic response. The ray-tracing is performed in the smooth reference medium m 0 : traveltimes, ray amplitudes and slowness vectors are delivered. These quantities of incident and diffracted rays make up the Green's functions (one for each couple source-target point and each couple receiver-target point).
For one source s, a seismic trace U is associated with each receiver r and corresponds to the displacement vector in function of time t. We note T the traveltime, Ampli 3D the ray amplitude and S(t) the emitted wavelet. The perturbation δm/m 0 is represented by relative perturbations in P-impedance δIp/Ip 0 , in S-impedance δIs/Is 0 and in density δρ/ρ 0 where the subscript 0 refers to the reference medium. Each trace U associated with a couple (s,r) is composed of the contributions of each point x Ω of the target Ω: where * denotes the convolution operator and δ denotes the Dirac function. The term Diffraction is expressed as: where θ is the diffraction angle derived from slowness vectors and γ 0 is the ratio of S-velocity and P-velocity at x Ω in the reference medium. We can note that the contribution of the point x Ω is only taken into account if the dip perpendicular to the sum of the slowness vectors has been selected in the modelling parameters [17].
The perturbations [δIp/Ip 0 , δIs/Is 0 , δρ/ρ 0 ] are discretized on fine parallelepiped meshes [δx, δy, δz]. On the other hand, the Green quantities are calculated on coarse parallelepiped meshes [Δx, Δy, Δz] of the target. In order to obtain the contributions of the perturbations in the trace, the Green   Surface acquisition pattern (sources in green, receivers in black) and target location (shadow zone, 820 m long by 100 m wide).
quantities are interpolated for the fine meshes using the slowness vector derivatives. The reflection events are correctly built for a dip thin plate discretized in parallelepiped meshes if the mesh is fine enough. More precisely, the (x, y, z) dimensions of the mesh must respect some criteria which depend on the significant wavelength λ of the P-wave propagated field: δx and δy are required to be less than λ/8, δz is required to be less than λ/16. In addition, several works have proved that the Born modelling is not accurate for large incidence angles [18,19,8].

Modelling Parameters
The features of the 4D elastic model are that the reservoir zone is very heterogeneous and that the overburden zone is close to a stratified medium. Each of the four models is defined in P-velocity, S-velocity and density from the surface to the bottom of the reservoir (Fig. 5). The extensions of the model are Y = 1 970 to 2 790 m and X = 712 to 812 m in our local coordinate system. The models are initially discretized with a sampling of 0.484 m in depth, 20 m in the Y-direction (parallel to the well paths), and 1 m in the X-direction (perpendicular to the well paths). The well paths are located around the depth of 310 m.
Acquisition parameters are selected to image a 100 m wide and 820 m long target, the target depth thickness remaining that of the reservoir model (close to 50 m). In doing so, the image will contain the entire well pair. We used an example of a dense land acquisition survey design, described in Figure 6. The source area located right above the target is included in the receiver area, which is larger. Sources and receivers are at the ground surface. The receivers record the vertical displacement of the computed wavefield. The source type considered is a vertical point force whose signature is a 5-to 320-Hertz band-pass. For each model, we handled a 3D seismic profile of 336 shots with 1 550 fixed receivers.
In the ray+Born modelling technique, the reflection traveltimes, as well as the geometrical spreading, are governed by the reference model, whereas the reflection amplitudes are mainly linked to the perturbations of the target. Special attention is paid to the smoothing of the full model because the surface-target traveltime differences between successive stages of production must be preserved. The four 3D SAGD well pair velocity models going from the surface to below the reservoir are smoothed one by one. Dealing with P-wave reflections, the smooth reference model must only be defined in P-velocity for the full domain.
At this stage of the work, we controlled that the traveltime differences between the four models remain after smoothing. We observed time differences in the range of 1 ms, which corresponds to the mean range that could be expected according to the P-velocity models converted in time.
Once the reference model is built, the P-impedance, S-impedance and density perturbations of the target are computed. For each model, the perturbations [δIp/Ip 0 , δIs/Is 0 , δρ/ρ 0 ] must be defined on a fine mesh in the target zone. The smoothed S-velocity and P-velocity ratio γ 0 is also needed in the target zone.
As we were dealing with a maximum frequency of 320 Hz and a mean velocity of 2 500 m/s in the target, we needed to take a fine mesh of dimensions δx = δy < 0.97 m and δz < 0.48 m. We chose δx = δy = 1 m and δz = 0.484 m in order to have to resample only in the Y direction.
The relative perturbations were then computed and we obtained average values less than 0.1 for P-impedance and S-impedance and less than 0.04 for density. But for the third year of production model, we obtained some individual S-impedance perturbation values greater than 0.1 (Fig. 8) because of some strong S-velocity contrasts between sands and shales at this stage of production. That is why we performed an additional test of the Born modelling validity in this case. We calculated one shot using the ray+Born modelling for a 1D column and for a 2D section of the model corresponding to the high S-impedance relative perturbations. These P-wave shot gathers were compared to the full-wave shot gathers modelled with the help of finite-difference modelling (Fig. 9,  10). Though some high values of S-impedance perturbations are present, we concluded that the Born modelling is still valid for near and medium offsets in our case.
Different sizes of coarse meshes were tested: 20 × 20 × 9.68 m 3 and 15 × 15 × 4.84 m 3 . We selected the lower size because 1D target results obtained by ray+Born modelling fit well with the ones obtained by 3D full-wave modelling in a stratified medium.
Another point that we paid attention to is the dip of reflectors to be simulated. The target model at the 3rd year of production presents the strongest dip on the edge of the V-shaped steam chamber. The dip could rise there to 30 degrees. So we took  Smoothed models (section X = 762 m parallel to the well paths) at initial stage (left), 3rd year of production minus initial stage (right): P-velocity (top), velocity ratio γ 0 (bottom).
into account, during the modelling process, a dip aperture limited to 40 degrees.

Seismic Results
For each of the four productions stages, shot gathers were then modelled according to the chosen acquisition geometry.
The generated seismic traces are thus free of noise and contain only primary reflections originating from the reservoir zone (no surface, multiple or converted waves).
The synthetic data obtained with the 3D ray+Born modelling are multi-offset and multi-azimuth seismic datasets. Each 3D shot gather contains 1 550 seismic traces. Figure 11 shows part of a shot gather located at the acquisition center. As the Born modelling breaks down for large incidence angles, it is necessary to deal with a limited offset range (0 m to around 600 m).
The models of the four stages of production used for the seismic modelling exhibited differences in the reservoir zone. As observed for models, there were more significant differences in the seismic data between the 3rd year of production and the initial stage. We looked at the differences induced in the synthetic shot gathers in order to verify that the 4D effect is present in time-lapse and in amplitude variation. It can be noted that the maximum time difference is very small (close to 1 ms). On the other hand, seismic amplitude variation is stronger (Fig. 12). Model with high S-impedance perturbation values: ray+Born (left) versus Finite differences (right). 1D model (column X = 752 m Y = 2 780 m, 3rd year of production). 130 Hz Ricker wavelet.

Adding Noise
Our aim was to add realistic noise on prestack seismic data. Two kinds of noise can be encountered for real data: random and coherent. The random noise is generally present with a signal to noise ratio lower for large offsets than for near offsets. The coherent noise often comes from a physical or mechanical event.
For the first work on these synthetic data, we chose to take into account only random noise. We used random noise data recorded by receivers during a real acquisition survey. Four different noise gathers were required for the synthetic prestack data (one gather per stage of production). For each stage, we added random noise on synthetic prestack seismic data with a frequency-dependent signal-to-noise ratio: around 0 dB for frequency 50 Hz and 10 dB for frequencies 120-200 Hz (Fig. 13).

Process Workflow
The noise-free data are processed to get stacked and time migrated images in order to be able to run a seismic interpretation study in relation with the SAGD production as described    Time (s) Figure 12 Shot gather located at Y = 2 300 m, X = 767 m: initial stage (black color) and 3rd year of production (grey color). Comparison of amplitude spectrum of modelled shot (blue color) and noise data (grey color) at initial stage: difference of 0 dB at 50 Hz and of 10 dB at 200 Hz.
in the second part of this paper. We used a commercial software for the processing of the 4D seismic dataset. The synthetic data obtained with the 3D ray+Born modelling contain only the P-wave reflections that occur in the target reservoir zone. That is why we could perform a simple but robust processing workflow for each stage of production: -first, we binned the traces with a bin size equal to 5 × 5 square meters. The seismic dataset was then sorted by Common Mid Point (CMP) position defined by X and Y locations; -secondly, we applied a normal moveout correction on each CMP with a RMS velocity model derived from the smoothed P-velocity model used in the modelling. The smoothed velocity differences of around -50 m/s between the third year of production and the initial stage led to RMS velocity differences of around -10 m/s. The corrected traces were stacked for the near offset range from 0 to 300 m; -finally, we needed to perform a 2D time migration on the stacked sections, as many diffraction events are visible in the Y direction. This fact was strongly visible for the 3rd year of production where the dip of the steam chamber edges is thirty degrees. We applied the migration on all the Y sections of the stack that are perpendicular to the well paths. In order to be able to study the effect of noise on interpretation results, the same processing workflow as for the noise-free data was performed on noisy data. But several "smiles" of migration were obtained due to the high energy of local noise events. That is why we needed to filter the low frequencies of the noisy synthetic data. In addition, it has been decided to filter the high frequencies in order to be closer to a real acquisition. A band-pass filter of 5-20-200-220 Hz has been applied to the noisy data before migration and also to the noise-free data in order to obtain comparable data.

Post-Stack Seismic Data
The processing workflow allowed us to obtain well focused reflectors for the noise-free data. The migrated stacks of the four stages are shown in Figure 14. There are no significant differences on the migrated stack between the initial stage and the first month of production for most of the CMP (Fig. 15). But at Y = 2 100 m, differences are visible just above the bottom of the reservoir around the location of the wells and could be related to the abnormal growth of the steam chamber (overpressured zone). Concerning the larger differences between the initial stage and the third year of production, they are more visible, on the lower right part (T ≈ 0.265 s) of the X = 760 m section, around the heel of the wells, and on the upper middle part, above the well pair around the top of the steam chamber (T ≈ 0.25 s). As expected, the migrated stack of the noisy data was more biased than for the noise-free data. Particularly, some noise residuals appear like reflectors at the initial stage and the first month of production stage (Fig. 16). The seismic interpretation presented in the second part of this paper will give more explanations of these differences. But, we must keep in mind that the large seismic coverage improves the resolution of both the seismic response and the seismic differences.

INTERPRETATION OF 4D SYNTHETIC SEISMIC DATA
Both noise-free and noisy synthetic migrated stacked seismic data were considered for interpretation, at the four available production stages. Let us remind that all seismic data are filtered with a band-pass of 5-20-200-220 Hz. In the following, we first briefly present the 4D synthetic seismic data as well as the reservoir/geomechanical properties issued from the model. Second, we present some of the interpretation results: -computation of some seismic attributes from the synthetic seismic data; -evaluation of the quantitative link (linear or not) of these attributes with the reservoir/geomechanical properties issued from the model; -geobodies.

Available Reservoir/Geomechanical and Seismic Data
All figures in this part and Section 3.2 present inline section X = 760 m of the seismic cube, from Y = 1 970 to 2 790 m (820 m long) and from time 225 ms to 275 ms: this inline section is along the well pair path. Figure 17 shows (at initial stage and 3rd year of production) some of the reservoir/geomechanical properties that will be used for the interpretation of the seismic attributes: lithofacies, temperature, pore pressure, oil and steam saturations. Other reservoir/geomechanical properties like mean effective stress, volumetric strain and water saturation were also available for interpretation. All these properties are issued from the 4D geological model in depth: they were transferred in time according to the associated velocity model, and then rescaled on the seismic grid. Figures 18 and 19 show the seismic amplitudes at the four production stages, respectively for noise-free data and noisy data. And Figure 20 present the seismic amplitude changes after one month and three years of production, computed on noise-free and noisy seismic data. As can be seen from Figure 18 and Figure 20, there are already some amplitude changes after one month of production at the heel of the well pair: these changes are linked with the presence of a shale level just above the injector well  Noise-free migrated stack Y-sections: initial stage, 1st month, 6th month and 3rd year of production (from top to bottom).

Seismic Data
(see Fig. 17a, b and Fig. 2). This shale level acts as a barrier and prevents the steam chamber to develop vertically; therefore, very early in the production, high pressure values appear above this level, due to the fact that fluids in the corresponding materials are heated by conduction but cannot be produced; as well below the shale level, high mean effective stress values and high temperature values are observed. This impacts both densities and velocities, thus seismic amplitudes. These seismic amplitudes and delays changes at the heel of the well pair increase with calendar time.
In the middle of the well pair, the reservoir is of very good quality, showing no shale levels above the well injector. In these areas, changes in the seismic character (amplitudes and delays) are mainly visible after the 6th month of production.
From Figures 19 and 20, it appears that the impact of the noise is very strong at all stages: for example the noise completely overrides the seismic changes at the beginning of production. Of course, on real seismic data, one would apply noise removal techniques before interpretation of the seismic data. However, the geometry and volume of the steam zone cannot be clearly delineated even on noise-free seismic data. Therefore, this call for the use of seismic derived attributes instead of the sole seismic amplitudes, when trying to interpret the seismic surveys in terms of geomechanical and reservoir properties along production.

Computation of Seismic Attributes
Several seismic attributes [20] were computed, e.g.: -3D attributes (energy, frequency filters, RMS values of amplitude changes between stages, attribute changes between stages, similarity); -2D attributes (maps of time differences between stages, reservoir traveltime changes between stages). Some of them are briefly discussed in the following in terms of production content, robustness to noise etc.

Energy
This attribute enhances, among others, lateral variations within seismic events. The response energy also characterizes acoustic rock properties and bed thickness.

Figure 16
Noisy migrated stack Y-sections: initial stage, 1st month, 6th month and 3rd year of production (from top to bottom).  Figure 21 shows the energy on the noise-free and on the noisy data at initial stage and at 3rd year of production. First, we see that this attribute is sensitive to the production (change of energy between initial stage and third year of production) and robust to the noise: even if the noise is visible, it is not strong enough to override the change linked with production. Second, by comparing this figure to Figure 17 (lithofacies and for example steam saturation), we see that energy combines both reservoir information on lithology (high values of energy corresponding to shaly content at initial stage) and production effects.
To get rid of the lithology effect on energy, the changes of energy between the three production stages and the initial stage (see Fig. 22) were computed. The high values  of energy changes after three years of production can be visually linked to the steam saturation envelope (Fig. 17i, j). Thus, a cutoff on energy change could give an idea of the steam area.
According to the way the noise was added to the seismic data (i.e. independently from one stage to another), the energy change appears slightly sensitive to the noise.

Frequency Filters
Nine band-pass filters were applied on the noise-free and on the noisy data at the four production stages. To sample the whole seismic frequency range (10 Hz to 220 Hz) and preserve the same octave number, these band-pass filters were chosen as triangular as possible with frequency peaks at 10, 14, 20, 28, 40, 57, 81, 115, 163 Hz.   Figure 23 shows the result for band-pass filters 28, 57 and 115 Hz on the noise-free data, respectively at initial stage and 3rd year of production. Results on noisy data are not given here. But from all results, one could see that: -there is very little information up to the 20 Hz peak; -the frequency filtered data are difficult to visually link with reservoir or production information; -there is no visible production effect below the 28 Hz peak; -above the 28 Hz peak, the production effect is increasing with frequency;  Seismic amplitude changes after one month of production (top), 3 years of production (bottom), computed on noise-free seismic data (left) and noisy data (right). Inline X = 760 m along well paths of length 820 m, between 225 ms and 275 ms.

Figure 21
Energy at initial stage (top) and 3rd year of production (bottom), computed on noise-free seismic data (left) and noisy data (right). Inline X = 760 m along well paths of length 820 m, between 225 ms and 275 ms.
-the noise strongly affects the seismic filtered data up to 20 Hz peak frequency and is clearly visible up to 57 Hz peak frequency (81 Hz maximum frequency), which is coherent with the noise frequency content; -the noise overrides partly the production effect for low frequencies.
So, from these 4D synthetic seismic data, the frequency filter seismic attribute appears sensitive to the production effect but visually unlinked with reservoir/production properties. It is also sensitive to the noise, mostly at low frequencies. And the noise may override the production effect on the seismic, thus biasing further interpretations.

RMS Values of Amplitude Changes Between Two Production Stages
Root Mean Squared (RMS) values were computed on seismic trace amplitude changes between production stages and initial stage, over a local neighbourhood. Figure 24 presents the RMS values of amplitude changes after one month of production and three years of production, on noise-free and noisy data.
We see that this attribute accounts for production. The higher values after one month of production (in green) correspond to the shale level at the heel of the well pair (see Fig. 17a, b) which prevents the steam to expand vertically. After three years of production, the higher values (in blue and purple) can be visually related to the highly steam saturated zone (Fig. 17i, j).
But unfortunately, this attribute appears very sensitive to the noise: this is critical at the beginning of production since the tiny production effects of the first month of production are completely overridden by the noise.

Computation of Time Differences Between Stages
In order to carefully estimate the time differences between any stage of production and the initial stage, the migrated stack  Changes of energy between 3rd year of production and initial stage, computed on noise-free seismic data (left) and noisy data (right). Inline X = 760 m along well paths of length 820 m, between 225 ms and 275 ms.  cross-correlations of the three stages of production with the initial stage were first computed. This was performed for a time window of 20 ms around the reservoir top and around the reservoir bottom; the time of maximum auto-correlation of the initial stage was then subtracted to the three times of maximum cross-correlation. Figures 25 and 26 show the corresponding time difference maps, respectively at top and bottom of the reservoir, between initial stage and three production stages, both for noise-free and noisy data.
These time differences between stages computed from the seismic data appear interesting to monitor the production, as well as are the reservoir traveltime changes not shown here. The general trend is that time differences are negative at the heel of the well pair (where a shale level prevents the steam chamber to expand vertically toward the surface), and positive in the middle of the well pair.
At reservoir top (Fig. 25), the production effect is only visible after six months of production with time differences less than 0.13 ms in absolute values. At reservoir bottom (Fig. 26), local differences are already visible along the well paths for the first month of production (comprised between -0.6 and 0.28 ms), and more accentuated for the sixth month of production (between -1.9 and 0.33 ms). After the 3rd year of production, time differences reach -2 ms at the heel of the well pair and 1 ms in the good reservoir area.
Noise has a very strong impact, mostly at the beginning of production and on the seismic cube borders (due to insufficient stack fold). Therefore, using these attributes on real seismic data calls for a careful removal of the seismic noise. This is even more crucial at the early stages of production with very tiny production effects. However, the reservoir bottom time difference map can obviously help to delineate, very early in the production, the areas where the steam chamber cannot grow properly like the shaly area at the heel of the well pair (with strong negative time differences).

Seismic Attributes Summary
This brief description of a few seismic attributes already shows that some seismic attributes are visually/qualitatively related to the lithofacies, but more importantly to the geomechanical/production properties. Working on real seismic data, one would better use 3D attributes like energy or energy changes between stages, or 2D attributes like maps of time differences between stages at reservoir top/bottom, since these seismic attributes are more robust to seismic noise. Other seismic attributes like 3D cubes of time differences between stages, amplitude envelop, instantaneous frequencies/phases or spectral components should be evaluated for comparison.

Quantitative Relationships Between Seismic
Attributes and Reservoir/Geomechanical Properties The next step was to check out for potential quantitative links between the geomechanical/production properties and the stages from noise-free data: the amplitude changes (named D_Ampl), the three frequency filters changes (named D_f28, D_f57, D_f115), the energy changes (named D_Energy) and the RMS values of amplitude changes (named RMS_D_Ampl). As can be seen, the relationships are of poor quality and cannot be used for estimating the geomechanical/production properties from the seismic attributes: the conclusion is the same for all others crossplots, things being even worst when considering noisy seismic data. Accounting for the lithofacies or the production stage does not improve the relationships. Thus, even if some seismic attributes appeared to be visually/qualitatively linked with the geomechanical/ production properties, they are not quantitatively related to them. This confirms the need for PetroElastic or Rock Physics Modelling when trying to quantitatively interpret 4D seismic data of heavy oil reservoirs produced with SAGD.

Geobodies
Two geobodies based on cutoffs (defined from the crossplots) on specific seismic attributes are proposed here. The first geobody is based on a cutoff on the energy changes between stages, keeping only values higher than 4.2e22. Figure 28 shows the corresponding geobody from noise-free seismic data for three years of production. As shown, this geobody corresponds to high values of temperature (greater than 157°C) and provides partly the top of the steam chamber. But it is slightly sensitive to noise (Fig. 29); and it is not efficient at the early stages of production, since no geobody can be extracted up to 6 months of production.
A second geobody has been established from a principal component analysis. A principal component [21] is a specific linear combination of attributes. In our case, the first principal component Prin1 is based on the seismic attributes changes between stages (mainly on seismic amplitudes and higher frequency filters at peaks 57 and 115 Hz) and can be expressed as:    stages, frequency filters between stages at peaks 115 and 57 Hz, RMS values of amplitude changes between stages, frequency filter between stages at peak 28 Hz, energy changes between stages. Keeping only values greater than 9.1 of this principal component provides the second geobody shown in Figure 30. This geobody corresponds to the very high steam saturation areas, and also to the higher values of temperature (greater than 179°C), higher values of mean effective stress, already partly produced areas (oil saturation 25% lower than over the whole reservoir). It is not efficient for detecting areas with steam saturations less than 50% (e.g. for the 1st month of production). Nevertheless, as observed from Figure 31, this geobody is stable to the noise which is of great interest if working on real seismic data. And moreover, this geobody proves the interest of the high seismic frequency content (in our case between 40 and 163 Hz, but mostly above 81 Hz) for delineating the very high steam saturation areas, and thus the potential of technologies such as permanent seismic acquisition as compared to traditional seismic acquisition.

CONCLUSION
The purposes of the work presented here were, first, to better understand the links between heavy oil reservoir SAGD production and associated seismic amplitude changes along calendar time, second, to investigate permanent seismic acquisition interest for monitoring the steam chamber growth in heavy oil reservoirs. This work was done on a reservoir model of such a Canadian Athabasca reservoir (Hangingstone field). It involved seismic modelling followed by interpretation.
The ray+Born seismic modelling allowed us to create 4D seismic data corresponding to different stages of production of the selected Canadian reservoir model. The pre-stack seismic data represent the P-wave response of the four elastic models at initial stage, at first month, sixth month and third year of production. The seismic amplitudes and the traveltimes are correctly handled for the near and medium offsets.
The synthetic seismic data that we modelled were clean data as they did not contain any multiple, converted or refracted events. We were able to obtain a large seismic fold order for a binning of five by five square meters because the acquisition pattern was dense. So, we were in an ideal case for the post-processing. We performed a time migration on stacked data which allowed us to image the V-shaped steam chamber development and to observe time-shifts in the reservoir zone. In this way, we could follow the evolution of the reservoir over time. In order to obtain more realistic seismic data, random noise was added to the pre-stack data. The same post-processing workflow as for noise-free data was applied on noisy data. The migrated sections were then of lower quality, some reflectors being biased. Nevertheless, the time and amplitude differences between the four stages of production partially appeared.
Other processing tests could be performed on these synthetic data. For example, we could take into account a sparser acquisition and study its impact on migrated sections. What could also be done is to estimate the RMS velocity in a blind context. Some other migration techniques such as 3D time migration or depth migration could be tested, as well as azimuth trace sorting. The analysis of amplitude versus offset or incidence angle could allow us to estimate the ratio of P-velocity and S-velocity which is a useful parameter for interpretation.
Moreover, full-wave seismic modelling should be useful to study S-wave parameters (steam effect?), as well as visco-elastic seismic modelling to study attenuation (heavy oil effect? weathered zone?). Based on the 4D modelled seismic data (noise-free or noisy), some interpretation aspects were studied in order to provide guidelines for real cases. A first look at the seismic data showed that the geometry and volume of the steam zone cannot be clearly delineated even on noise-free seismic data, thus calling for enhanced interpretation methods. Moreover, the impact of the noise appeared very strong at all stages: this real noise completely overrides the seismic changes at the beginning of production.
A few seismic attributes like energy, frequency filters etc. were then studied. Some of these seismic attributes appear visually/qualitatively related to the lithofacies, but more importantly to the geomechanical/production properties. When working on real seismic data, one would better use 3D attributes like energy or energy changes between stages, or 2D attributes like maps of time differences between stages at reservoir top/bottom, since these seismic attributes are more robust to seismic noise. In future works, others seismic attributes should be evaluated, like amplitude envelop, instantaneous frequencies/phases, spectral components, or 3D cubes of time differences between stages that could maybe help for delineating heterogeneities in the reservoir.
Even if visual links were found out between these seismic attributes and geomechanical/production properties, no quantitative relationship could be defined: this confirms the need for PetroElastic or Rock Physics Modelling in 4D seismic data interpretation of heavy oil reservoirs produced with SAGD.
However, based on the seismic attributes, two geobodies were proposed. The second one, which is stable to the noise, corresponds to the very high steam saturation areas. And as this geobody accounts for high seismic frequencies in its construction (in our case between 40 and 163 Hz, but mostly above 81 Hz), it proves the potential of technologies providing high frequency contents for monitoring the steam chamber growth in heavy oil reservoirs, as compared to a traditional seismic acquisition of much lower frequency content. The next step would be now to confirm all these results on real 4D seismic data.