4D Joint Stratigraphic Inversion of Prestack Seismic Data: Application to the CO2 Storage Reservoir (Utsira Sand Formation) at Sleipner Site

Résumé – Inversion stratigraphique jointe 4D de données sismiques avant sommation : application au réservoir de stockage de CO 2 (Formation Utsira) du site de Sleipner – Le monitoring sismique est couramment utilisé dans l’industrie pétrolière pour suivre l’évolution des propriétés des réservoirs au cours de leur production. Nous présentons ici une méthodologie d’inversion stratigraphique 4D qui fournit une estimation des variations d’impédances des ondes P et S dans le réservoir, par inversion de données de sismiques répétées avant-sommation. L’inversion 4D est implémentée dans le domaine temps et nécessite une loi de mise à l’échelle des temps de trajets pour chaque jeu de données de sismique répétée (aussi appelé “millésime”) afin de mettre en correspondance les temps d’arrivée d’événements homologues observés sur les jeux de données appelés référence et monitors . Cette opération est souvent appelée “ warping ”. L’inversion 4D est une méthodologie qui comporte trois étapes : la première étape consiste à inverser chaque millésime sismique séparément, pour produire autant de distributions en impédances P et S que de jeux de données considérés. La seconde étape utilise l’information des impédances P disponibles pour résoudre le problème du warping , Abstract – 4D Joint Stratigraphic Inversion of Prestack Seismic Data: Application to the CO 2 Storage Reservoir (Utsira Sand Formation) at Sleipner Site – Seismic monitoring is commonly used in the oil industry to follow the evolution of reservoir properties during production. We present a methodology of time-lapse (or 4D) stratigraphic inversion, which is able to provide an estimation of P- and S-wave impedance variations in the reservoir by inverting prestack time-lapse seismic data. The 4D inversion is implemented in the time domain and requires a time scaling law for each repeated seismic dataset in order to adjust the arrival times of homologous events observed in the so-called reference and monitor datasets. This operation is often referred to as the warping problem. The 4D inversion is a 3-step methodology.


INTRODUCTION
One of the main advantages of time-lapse (or 4D) seismic is the ability to follow the evolution of the reservoir properties over time while continuing hydrocarbon production. For that reason, it is now a well-established and mature technology. Hydrocarbon production induces fluid substitution and possible pressure variations which are likely to modify the amplitudes and traveltimes of seismic signals recorded during monitor seismic surveys. Fine analyses of these changes can give valuable information about the time evolution of the reservoir and guide the optimal positioning of future wells to improve hydrocarbon production.
The geological storage of CO 2 in reservoir layers induces the same kind of changes in the seismic properties of the host formation. In this case, 4D seismic can greatly help to follow the CO 2 plume expansion and migration.
4D seismic technology requires a careful pre-processing of the seismic datasets to compensate for differences in acquisition pattern and recording parameters (source signature, frequency content) between successive campaigns. In our model-based stratigraphic inversion approach, seismic data are time-migrated prior to being used to build limited angle stacks as inputs to the inversion procedure.
A second crucial element of 4D technology is the warping correction, an operation designed to adopt a common time scale to match the reflections in one dataset to their counterparts in a former or subsequent dataset. To solve this problem, Martinson and Hopper (1992) developed a nonlinear correlation algorithm to align two adjacent seismic traces by determining a scaling factor, allowing, within some limits, to interpolate seismic traces while keeping the dip and amplitude of the seismic reflectors. This "cross-correlation" method was adapted for time-lapse seismic data: it consists in searching the time delay allowing two signals to be superimposed with a minimal difference. It is based on a similarity criterion between two signals and it gives good results when distortion between signal amplitudes is weak. This method belongs to a family of techniques that exploit the waveform similarity of two traces. Wolberg (1990) and Hall et al. (2002) show two applications of this type of time-alignment. To align two seismic traces, some other methods are based on a global nonlinear optimisation, such as the Needleman-Wunsch algorithm (Needleman and Wunsch, 1970).
The alignment methods presented above take into account the cinematic effects only, but are free of any physical constraints: the dynamic effects (amplitude variations) are not considered. To face this kind of problem, Williamson et al. (2007) implemented an approach that links velocity changes and time shifts to align amplitudes of time-lapse seismic data.
Over the past few years, IFP Energies nouvelles has developed another family of warping algorithms based on elastic impedances resulting from single stratigraphic inversions. The latter integrate geological information such as mode of sedimentary deposition and well log data (Tonellot et al., 1999). The technique has been successfully applied to S-wave impedances obtained by independent stratigraphic inversions of multi-component data, namely PP and PS synthetic seismic data (Agullo et al., 2004). The alignment methods rely on the computation of a scaling factor by minimizing the traveltime dissimilarity of the impedances expressed in different traveltime bases. Nevertheless in this case, the method does not take into account the fact that the time-scaling factor depends on the P-to S-wave velocity ratio and should have a realistic physical value.
A further evolution of the IFP Energies nouvelles methodology was inspired by scanner technology used in medical imaging to follow 3D displacements in the human body. The warping problem consists in finding a scaling factor using 3D spline surfaces (Kybic and Unser, 2003). This technique was adapted in a trace-by-trace algorithm applied to P-wave impedances obtained through single stratigraphic inversions of base and monitor vintages (Tonellot et al., 2010). In this approach, the time scaling factor is a function of time-shift and impedance values which depend directly on P-wave velocity (assuming an unchanged density). The estimated time-scaling law is then used to jointly invert all time-lapse seismic data. many P-and S-wave impedance distributions as datasets considered. The second step uses the available P-wave impedance information to solve the warping problem which is crucial to the third and final step: the joint inversion of all available seismic vintages. This 4D inversion sequence was applied to seismic datasets recorded on the Norwegian CO 2 storage reservoir of Sleipner field located in the North Sea. The latter is becoming a reference industrial site for the long-term storage of carbon dioxide (CO 2 ) in a saline aquifer, the Utsira sand formation. We focused our 4D inversion study on the 1994 and 2006 vintages acquired two years before and ten years after the beginning of CO 2 injection, respectively. The warping correction resulted in a time-scaling law with a maximum pushdown effect of about 45 ms at the base of the Utsira aquifer. The joint 4D inversion results show more consistency than the single 3D inversion results: the 4D inversion notably provides P-wave impedances for the CO 2 -saturated sandstones which are close to the values derived from rock physics studies.
Concerning inversion itself, some authors have noticed an improvement by introducing a coupling between the vintages in the seismic inversions, for example through the inversion of the amplitude differences between the baseline survey and the monitor survey (Buland and El Ouair, 2006;Sarkar et al., 2003;Lafet et al., 2009), or at least through the use of the baseline inversion result as an initial model for the monitor survey inversion (Sarkar et al., 2003). Notice that this coupling requires the preliminary time-alignment of the base and monitor seismic data before proceeding with the rest of the procedure. IFP Energies nouvelles methodology benefits also from joint inversion of seismic vintages in order to use redundancy of information in geological units where fluid injection or production has no impact (sometimes below the reservoir unit for example, or in the overburden). Compared with other approaches, IFP Energies nouvelles methodology does not require time-alignment before inversion: scaling law is taken into account during the inversion process.
The above 4D inversion was applied to the storage reservoir at Sleipner site (located offshore Norway) where carbon dioxide (CO 2 ) is injected in a saline aquifer in supercritical state for a long-term storage (Chadwick et al., 2009). Seven seismic surveys were carried out between 1994 (before CO 2 injection) and 2008 to monitor the time evolution of the CO 2 plume. Rabben and Ursin (2011) insist on the processing of seismic amplitude in order to obtain reliable seismic amplitudes for seismic inversion process. 2D full wave form inversions were performed at Sleipner by Gosselet and Singh (2008) but this approach is computer intensive. 3D single prestack stratigraphic inversions were recently performed using the 1994 and 2006 seismic data . These inversions yielded P-wave and S-wave impedance cubes which turned out to be very useful to delineate the CO 2 plume in a much better way than the seismic amplitudes. The P-wave impedances were used to solve the warping problem.
The massive CO 2 injection in the CO 2 storage reservoir at Sleipner site (more than 8.4 Mt up to 2006) is responsible for a drastic decrease of P-wave velocity producing misalignments of the seismic events located inside and below the CO 2 plume compared to the reference 1994 seismic data, before CO 2 injection (Arts et al., 2008). Another difficulty of the CO 2 storage reservoir at Sleipner site is the accumulation of CO 2 beneath shale layers generating strong reflections quite different from those of the 1994 vintage. To overcome these issues, we applied the patented warping method of Tonellot et al. (2010) which can handle strong pushdown effects and amplitude variations. The resulting time-scaling law was subsequently used in the joint stratigraphic inversion of the 1994 and 2006 vintages. We present here the joint inversion results of the base and monitor surveys as well as their interpretation in terms of impedance variations between vintages. The 4D inversion results are compared to the single inversion result of 1994 and 2006 vintages showing that the joint inversion methodology benefits from redundancy information below and above the reservoir, where CO 2 injection has no impact.

THREE-STEP TIME-LAPSE STRATIGRAPHIC INVERSION WORKFLOW
The joint model-based stratigraphic inversion methodology of prestack seismic data used in this work (Tonellot et al., 2010) consists of three steps summarised in Figure 1 and explained below.
Step 1: Sequential Stratigraphic Inversions The stratigraphic inversions run at step 1 follow the methodology described in Tonellot et al. (2001). Each dataset (baseline and monitors) is inverted in its own traveltime basis. The results of these inversions are elastic earth models parameterized in terms of P-wave and S-wave impedance, one for the baseline dataset and the others for the monitor datasets. Therefore, at this step, no direct comparison (such as impedance variation estimations) can be performed between the elastic earth models.
Step 1: Sequential prestack stratigraphic inversions of the reference baseline and the n monitor seismic vintages Step 2: Estimation of the time scaling laws t 1 (t 0 ), ..., t n (t 0 ) related to velocity variations → PP-traveltimes of base and PP-traveltimes of monitor are related → Time shift estimations Δt 1 = t 1 − t 0 , ..., Δt n = t n − t 0 , from impedance results Step 3: Joint inversion of all 4D seismic data represented on various traveltime basis → Joint re-estimation of optimal impedance and density models using time shifts and applying common constraints → Optimal impedances and density models estimation (l P , l S , ρ) in their own time basis (t 0 , t 1 , ..., t n )

Figure 1
Time-lapse stratigraphic inversion workflow in three steps. For each step, the result is indicated italicized. t 0 stands for two-way traveltime for the base survey, t i (i ≠ n) stands for the two-way traveltime for the monitor survey.
Step 2: Estimation of the Scaling Law Estimation of the scaling law, also called "warping", is performed using the P-wave impedances obtained at step 1 from the base and monitor surveys by computing a time scaling law between the two surveys. This step does not only take into account the time shifts between the impedance distributions but also accounts for impedance variations related to velocity changes. The estimated time scaling law then provides time shifts between base and monitor traveltimes that are constrained with physical variations of the Pwave velocity. In our methodology, warping is formulated as a differential equation that links time shifts and velocity changes. The problem is solved by an optimization of a non-linear quadratic problem. For each trace of the seismic data cube, we seek to determine the function t 1 (t 0 ) (where t 0 is the time of a given seismic event in the baseline survey, and t 1 is the time of the same seismic event in the monitor survey) that minimizes the objective function F (1): (1) where v 0 and v 1 are the P-wave velocities in the base and monitor datasets, respectively, and α is such that: i.e., t 1 (t 0 ) is constructed from an initial profile t 0 perturbed by a quantity d depending on parameter α.
As input of the warping process, the initial time profile t 0 can be chosen from a linear interpolation between the time pickings of two horizons, or trivially as the identity function.
The implementation of the warping method presented here using the P-wave impedances obtained in step 1 assumes that density changes are negligible compared to P-velocity changes. This assumption is reasonable in two cases: either P-wave velocity decrease is predominant over density decrease (for example, when gas is injected into a saline aquifer), or density variations are indeed very weak compared to velocity variations.
Compared to other methods directly using the seismic amplitudes (Williamson et al., 2007 for example), our approach is more robust because inverted P-wave impedances benefit from the integration of a priori geological knowledge in the sequential inversion process. Therefore, the estimated impedances are less noisy and exhibit a higher frequency content; they also show more lateral continuity than seismic amplitudes. Besides, variations of incidence angle due to 4D effects do not need to be taken into account since information coming from partial angle stacks is integrated into inverted P-wave impedances.
Step 3: Joint Stratigraphic Inversion The joint stratigraphic inversion consists in inverting two or more time-lapse prestack seismic datasets simultaneously.
The inversion methodology requires a priori information which integrates various kinds of knowledge recorded at different scales: well log information, Vertical Seismic Profiles (VSPs) if available, picked seismic horizons, stacking velocities which are combined to build a priori 3D elastic models for each vintage. We describe the methodology for one baseline and several monitor vintages. The final results are 3D optimal models parameterized in P-wave impedance, S-wave impedance and density, consistent with all the input data. Similar to the poststack (Brac et al., 1988) and prestack (Tonellot et al., 2001) stratigraphic inversion techniques, the joint stratigraphic inversion methodology is based on a Bayesian formalism, where uncertainties are described by Gaussian probabilities with covariance operators C d in the data space and C m in the model space.
The stratigraphic inversion is formulated as a non-linear least-squares local optimization problem which is iteratively solved using a conjugate gradient method. The global objective function to minimize for the joint inversion of the baseline and n monitors is defined by Equation (2): (2) with: • , computed in the t i two-way traveltime basis; • , computed in the t 0 twoway traveltime basis; • t 0 and t i linked through the time scaling law estimated at step 2. i subscript refers to the considered seismic vintage, either the baseline (0 subscript) or monitors (i subscript for i ≠ 0). All m models are described and updated in the base traveltime basis (t 0 ). However, in order to compute J i Seis seismic terms, monitors are mapped in their own traveltime basis (t i ) using the time scaling law (estimated at step 2). • For each seismic vintage, the seismic term J i Seis measures the misfit between the real seismic data (d θ obs ) and the synthetic seismograms (d θ synth ) computed with the m i current model. Note that the seismic data are always expressed and computed in their own t i traveltime basis to avoid any transformation of the initial pre-processed seismic data. Synthetic seismic data are computed for a given incidence angle by a 1D-convolution model: where R(m, θ) is the Aki-Richards reflection coefficients series (Aki and Richards, 1980) corresponding to the current elastic model m at incidence angle θ; and W θ is the optimal wavelet for that angle.  (Utsira Sand Formation) at Sleipner Site C d is the matrix describing the uncertainties on seismic data. Because seismic noise is assumed to be uncorrelated from one trace to another, the C d matrix is diagonal, with a variance that is a function of the signal-to-noise ratio in the seismic data. • The "geological" term J i Geol measures the misfit between the a priori model m_prior i and the predicted impedance model m i according to the norm defined by the inverse of the model covariance operator C m . Note that all a priori and current models (base and monitors models) are described and updated in the t 0 baseline traveltime basis. The elements of the covariance matrix C m are exponential operators which express the confidence the user has: -in the a priori model geometry through a correlation length parameter, -in the a priori model values through the standard deviation which is the deviation of the optimal model relative to the a priori model. The developed methodology offers advantages compared with some other approaches: -seismic data always remain in their own traveltime basis and seismic amplitudes are preserved from pre-processing until limited angle class building; -the warping provides time-alignment which guarantees realistic physical results taking into account the pushdown effect (due to CO 2 injection) and P-wave impedance variations; -the joint inversion methodology is able to invert globally all vintages, assuming the warping is solved for each monitor vintage; -3D constraints can be introduced in some specific stratigraphic units where no changes are expected due to the field operations (for example below the reservoir unit).

CASE STUDY: THE CO STORAGE RESERVOIR AT SLEIPNER SITE
The Utsira Sand formation is a relatively shallow saline aquifer located in the North Sea, on the Sleipner field (Fig. 2). It is used for long-term underground CO 2 storage. Since 1996, more than 11 million tons of CO 2 have been injected in this formation. The Utsira Sand reservoir is made of poorly consolidated sandstones (Zweigel et al., 2004) of high porosity and high permeability with thin shale intra-layers inside (Arts et al., 2008). It is topped by a thick shaly layer.
Seven seismic data acquisitions were performed over the injection area between 1994 and 2008 in order to monitor the evolution of the CO 2 plume. The seismic images obtained highlight the strong effect of CO 2 injection on the seismic response, both on seismic amplitudes and traveltimes inside the CO 2 plume (Eiken et al., 2000).
The application of IFP Energies nouvelles methodology presented in the last paragraph allows a quantitative comparison of impedance variations to be performed.

Input Data
The required input data for our time-lapse stratigraphic inversion workflow are seismic data, well log data, picked seismic horizons and stratigraphic information. We considered seismic datasets of 1994 (before injection) and 2006 (after ten years of injection). We chose 2006 dataset because it was the last available dataset when our study started and it showed the CO 2 -related maximum pushdown effect.
The used seismic data consist of three limited angle stacks for each seismic vintage as well as the P-wave velocity cubes obtained by velocity analysis. More precisely, the seismic data are time-migrated, NMO-corrected data cubes which were partially stacked over the following incidence angle ranges: [6°-16°], [17°-27°] and [28°-38°]. The data considered have 249 × 468 inline and crossline dimensions. Seismic bin size is 12.5 m × 12.5 m and time sampling rate is 2 ms. The P-wave velocities are used to constrain the very low frequencies (up to several Hertz) of the a priori model, as described in Nivlet (2004) and Delépine et al. (2009).
Two sets of four horizons were picked on seismic data over the area (Fig. 2): one set on the 1994 vintage and another set on the 2006 vintage. The Utsira Sand formation is delineated by its top (Top Utsira) and base (Base Utsira) horizons. The Base Utsira horizon was difficult to pick on the time-migrated seismic data of 2006 because of the defocusing of the seismic image due to the presence of CO 2 in supercritical state. Two additional horizons were picked below and above the saline aquifer at a time of about 720 and 1330 ms, respectively. They are used to bound the inverted time window.
The chosen horizons are used to delineate the geological units of the a priori models.
This structural information is completed by information on the sedimentary deposition mode inside each structural unit, which eventually allows us to define seismic surface correlations that will be used in the construction of the a priori geometrical model. The seismic wavelets required by the inversion process for each angle stack are derived from well log information. On the CO 2 storage reservoir at Sleipner, only two wells are available within the area of interest for our study (Fig. 2): -well 15/9-13, a vertical exploration well; -well 15/9-A16, a horizontal injection well which is only partly located inside the 4D seismic area. P-wave sonic and density logs are available for these two wells. This information is completed with S-wave sonic data coming from a third well (15/9-A23) located outside the area of interest in the vicinity of the area encompassed by the 3D seismic survey. At the top Utsira, the interwell distance between well 15/9-A23 and well 15/9-13 is about 1 600 m. However, since it is the only well on the CO 2 storage reservoir at Sleipner with S-wave information, we have used it.

Results of the Time-Lapse Stratigraphic Inversion Workflow on the CO 2 Storage Reservoir at Sleipner
The results of the joint inversion of the Sleipner seismic data are presented according to the three steps of the methodology introduced earlier.
Step 1. Sequential Stratigraphic Inversion Results Two 3D stratigraphic inversions were separately performed with the 1994 and 2006 seismic data to estimate the P-and S-wave impedances. The complete inversion parameters and results of these two inversions are presented in . Figure 3 shows the P-wave impedance results for a cross-section through the CO 2 plume. Clochard et al. (2010) stress that prestack stratigraphic inversion provides optimal elastic impedance models which are very useful to characterize the CO 2 plume. For instance, the S-wave impedance distribution is used to assess the pushdown effect associated with the gas injection. As a first approach, it can be used to delineate the extension of the CO 2 plume.

Step 2. Warping Results
Warping is performed with the P-wave impedance results derived from step 1 (Fig. 3). Application of our warping methodology is legitimate because in the case of CO 2 storage reservoir at Sleipner, CO 2 is injected at supercritical state in a saline aquifer, and therefore density changes are negligible compared to P-wave velocity changes.
It results in the time scaling law presented in Figure 4, which represents a 2D slice of a 3D distribution of time shifts between the 1994 and 2006 seismic data (in the 1994 traveltime basis). It is to note that inline 1 832 (Fig. 4) intercepts the openhole section of the 15/9-A16 injection well. Figure 4 displays the pickings of the Top Utsira and Base Utsira horizons used during the warping process. As expected, the Top Utsira pickings are similar in 1994 and 2006, however, the Base Utsira pickings are different in 1994 and in 2006 because of the presence of the CO 2 plume. Figure 4 shows positive time shifts that can reach up to + 44 ms in the CO 2 plume. A positive time shift means that the 2006 seismic velocity is lower than the 1994 seismic velocity, an observation consistent with the injection of supercritical CO 2 in the reservoir.
Some negative values are nevertheless observed in Figure 4. This can be explained by artefacts of the seismic processing possibly related to the migration step, because a different velocity model was used for the two vintages. Optimal P-wave impedance models (g/cm 3 .m/s) obtained for inline section 1 832 from the 1994 (top) and 2006 (bottom) seismic data as a result of step 1 of the workflow (sequential stratigraphic inversions). The two inversion results are presented in their own traveltime basis.  Step 3. Joint Inversion Results In order to compare the joint stratigraphic inversion results and the sequential stratigraphic inversions, it is necessary to use same inputs (same optimal wavelets for each angle stack) and same inversion parameters (such as uncertainties on the a priori model). Joint inversion results are presented here with a number of fifty iterations, that is exactly equal to the number of iterations used for the sequential inversions. For more details about the sequential inversions on the CO 2 storage reservoir at Sleipner, readers should refers to Clochard et al. (2009). Because no change is expected below Base Utsira horizon during CO 2 injection, joint inversion is performed with the constraint of finding a common model below Base Utsira horizon. Figure 5 presents the joint inversion results along inline section 1 832 for the P-and S-wave impedances expressed in the same traveltime basis. Below the saline aquifer, the same results are obtained for the two vintages; this result was imposed by the inversion parameters. In the saline aquifer, the CO 2 plume can be seen very clearly on the 2006 results, particularly on the P-wave impedances with low P-wave impedance values (in green). Concerning estimations of the velocity and density in a sandstone saturated with CO 2 , different works have been previously done. For example, Carcione et al. (2006) consider a sandstone with a porosity of 35% with a water and a CO 2 saturation of respectively 40% and 60%; in this case, the impedance obtained is 2 200 g/cm 3 .m/s. This value is very near from those we obtained using inversion in the CO 2 plume.
Outside the plume, results of the two vintages are very similar in and above the saline aquifer. Such observation is very clear in Figure 6, where the variations in P-and S-wave impedances between 1994 and 2006 are displayed for the same inline section. The negative variation in P-wave impedance observed between the two horizons and the crosslines 1 040-1 200 is consistent with the injection of CO 2 at a supercritical state in a saline aquifer, since the associated decrease of P-wave velocity has a stronger effect than the accompanying decrease in density. This observation is also consistent with the push-down effect observed on the seismic data and with the result of the warping operation (Fig. 4).
In addition, Figure 6 shows relatively small changes in S-wave impedance in the CO 2 plume. This can be explained by variations of density before and after the CO 2 injection since no significant pressure increase was observed in injection well head on the CO 2 storage reservoir of Sleipner field (Alnes et al., 2008).
We note a decrease of P-wave impedance values of about -1000 g/cm 3 .m/s in the area invaded by the CO 2 plume, but also an increase of 500 g/cm 3 .m/s seemingly associated with the shale layers present inside the Utsira formation (Fig. 6).
The observed increase of P-wave impedance in the shale layers is abnormal because no pressure increase has been noticed on the CO 2 storage reservoir of Sleipner field. This artificial increase of P-wave impedance in the shale layers is explained by a lack of low frequencies in the seismic bandwidth and by the zero-mean property of seismic data: during the inversion process, a decrease in impedance is partially compensated by a corresponding increase in impedance. Such effect is attenuated by using a low frequency trend introduced in the P-wave impedance a priori model with the use of the stacking velocities. Comparing with previous work in 3D post-stack inversion (Delépine et al., 2010), such an effect is much more attenuated in our 4D pre-stack inversion results. To suppress it totally, additional constrains could be used in the inversion process.
As seen in Figure 7, the joint inversion shows weak seismic waveform residuals at the end of the inversion procedure compare to the residuals obtained after the first iteration of the inversion process.
The comparison is done between the 4D inversion results in P-wave impedance with the 3D impedance results warped in the 1994 time basis (Fig. 8). 4D inversion results show less P-wave impedance variations than 3D inversion results. More precisely, the lowest values are around 1 750 g/cm 3 .m/s in the 2006 sequential inversion results and 2 200 g/cm 3 .m/s in the 2006 joint inversion results. The highest values are about the same in the sequential and joint inversions (around 4 500 g/cm 3 .m/s), but are less frequent in the joint inversion results. We can explain this by the zero-mean effect mentioned previously. This effect is attenuated in the joint inversion process because more seismic data cubes are taken into account during the inversion process. Events in the 4D case seems to be better aligned especially below the saline aquifer. Another advantage of the 4D inversion is to use several seismic vintages to characterize the same part of the impedance model: a common problem such as the defocusing effect below the saline aquifer -due to the presence of CO 2 -is solved by using simultaneously the monitor and reference seismic vintages.

CONCLUSIONS
We developed a new time-lapse stratigraphic inversion methodology that jointly inverts prestack seismic data of different vintages for P-and S-wave impedances.
Compared with independent stratigraphic inversions of each vintage, the new methodology offers some advantages and improvements, namely: -the possibility of constraining the geological units which are supposed to remain identical in all vintages, in order to get the same optimal earth model in the areas where field operations have no effect; -a robust determination of a time-scaling law from P-wave velocities. Starting from a priori models of different vintages in the reference time basis, the warping allows us  Figure 6 Results from the joint stratigraphic inversion for inline section 1 832: display of the impedance variations (g/cm 3 .m/s) between 1994 and 2006 for the P-waves (bottom) and S-waves (top). to transform such models in their own traveltime basis to compute synthetic seismic data. Real seismic data always remain in their own traveltime basis and therefore have never been distorted by any warping; -working in the same traveltime basis allows us to obtain a quantitative estimation of impedances variations; it is a step further to reservoir quantification. We illustrated this time-lapse stratigraphic inversion methodology with seismic data acquired over the CO 2 storage reservoir at Sleipner site. The joint inversion procedure was shown to provide better results than independent inversions, notably by giving more realistic P-wave impedance values. Reservoir characterization and the gas injection scenario on the Sleipner field could therefore benefit from this new methodology. P-and S-wave impedance variations between 1994 and 2006 are crucial information to perform a quantitative interpretation of CO 2 injection and to estimate the in situ volume and/or mass of CO 2 (Dubos-Sallée and Rasolofosaon, 2010).
The results can also be used as a constraint of the history matching in reservoir simulation.