Seismic imaging of the carbon dioxide gas cloud at Sleipner using 3D elastic time-lapse full waveform inversion

A solution to the increasing greenhouse gas emissions is carbon capture and storage. At the Sleipner field, offshore Norway, carbon dioxide is separated from the produced gas and injected into the Utsira formation at approximately 1000 m depth. The Utsira formation is a high porosity saline aquifer that is 200–300 m thick in the injection area. In 1994, two years before the injection started, a seismic survey was acquired. The survey has been followed up by several seismic surveys acquired to monitor the migration of the injected gas. The full waveform inversion (FWI) method is a non-linear inverse method for estimating subsurface elastic parameters. Time-lapse FWI has recently been suggested as a monitoring tool for resolving time-lapse changes directly in the elastic parameter models. This paper presents an application of a three dimensional implementation of isotropic elastic time-lapse FWI applied on time-lapse field data from the Sleipner area. The baseline dataset is the 1994 dataset, whereas the monitor dataset is the survey acquired in 2006 after ten years of injection. Time-lapse FWI is used to estimate separate elastic parameter models for both the baseline and the monitor datasets. The results show that FWI is able to estimate detailed elastic parameter models that produce synthetic data that match the field data. The inverted models are used in a pre-stack depth migration method to yield detailed seismic images of the Sleipner before and after the injection of the gas. The gas cloud after ten years of injection is clearly visible on the monitor seismic images. The baseline and monitor seismic images show events and discontinuities that can explain the migration pathways of the injected gas through the Utsira formation. © 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license


Introduction
Carbon capture and storage (CCS) has been proposed as one of the solutions for lowering the greenhouse gas emissions to the atmosphere. The idea is to capture the greenhouse gases before emission and inject them into subsurface geological formations for safe storage into the future. Large-scale sites have been demonstrating the CCS technology for nearly two decades. At the Sleipner field, offshore Norway, carbon dioxide (CO 2 ) has been stripped off the produced natural gas and injected into the Utsira formation (Furre and Eiken, 2014). At the Snøhvit field in the Barents sea, CO 2 has been injected into the subsurface since 2008 . At In-Salah in Algeria, CO 2 has been injected onshore since 2004 .
Time-lapse seismic monitoring has emerged as an effective tool for monitoring hydrocarbon reservoirs during production (Landrø et al., 1999;Lumley, 2001). Recently, this method has been adopted for monitoring purposes at CO 2 storage sites (Furre and Eiken, 2014). Conventional methods for quantifying time-lapse seismic effects in seismic data involve pre-stack time-migrated data cubes in the migrated image domain (Greaves and Fulp, 1987;Landrø et al., 1999). In cases where the time-lapse changes induce strong changes in the rock properties, the linear assumption between the baseline and monitor model in these methods is violated. In addition, the migrated data cubes yield images in the time domain, which requires additional time-to-depth conversion. Therefore, more robust time-lapse analysis techniques are required to give reliable time-lapse images.
The full waveform inversion (FWI) method is a non-linear inverse method for estimating subsurface elastic parameters (Tarantola, 1984;Mora, 1987). FWI is formulated as an iterative optimization problem that seeks to minimize the misfit between field data and synthetic data (Tarantola, 1984;Mora, 1987;Fichtner et al., 2006). The method is, in general, computer intensive since synthetic data must be computed in each iteration of the optimization method. In addition, finding the model update is at least as costly as generating the synthetic data. Therefore, the first applications of the method were limited to acoustic wave theory and two dimensional numerical domains (e.g. Pratt, 1999;Dessa and Pascal, 2003;Chironi et al., 2006). As the available computer resources increased, the FWI framework was extended to involve elastic wave theory (Shipp and Singh, 2002;Sears et al., 2008Sears et al., , 2010. Elastic wave theory is favorable in FWI since it produces synthetic data that include more accurate physics that in the end yield synthetic data that is closer to the field data. The computational effort, however, is higher than for acoustic theory due to a more complicated equation to solve. To perform FWI on field data using two dimensional computational domains, the field data must be scaled to fit into the framework (Pratt, 1999;Ravaut et al., 2004). It turns out that this scaling may yield artifacts in the inverted models (Auer et al., 2013). Therefore, three-dimensional numerical domains are preferable. However, performing FWI in three dimensions has until recently been impossible in practice due to the required run times for the numerical methods used to generate the synthetic data. Sirgue et al. (2008) and Vigh and Starr (2008) demonstrated three dimensional implementations of FWI using both synthetic and field data. Their work has been followed up by several other studies (Sirgue et al., 2009;Plessix, 2009;Warner et al., 2013;Butzer et al., 2013;Vigh et al., 2014). In global seismology, three dimensional implementations have been used to create images of the crust of the Earth (Tape et al., 2010;Fichtner et al., 2010Fichtner et al., , 2013. Different approaches for performing time-lapse FWI (TLFWI) exist in the literature (Watanabe et al., 2004;Zheng et al., 2011;Raknes and Arntsen, 2014a). Compared to the conventional methods mentioned above, TLFWI is able to resolve time-lapse changes directly in the parameter models, in addition to yield results in cases where the linearity assumption is violated. Watanabe et al. (2004) used FWI on synthetic time-lapse cross well data to resolve time-lapse changes in a small area. Zheng et al. (2011) investigated different strategies for performing FWI on synthetic ocean-bottomcable (OBC) seismic time-lapse data. Queißer and Singh (2013) performed two-dimensional inversion of time-lapse data from the Sleipner area to investigate the migration of the injected CO 2 gas. Borisov and Singh (2013) applied 3D isotropic elastic TLFWI on synthetic data using the grid injection method. Raknes and Arntsen (2014a) used time-lapse data from a blow-out well in the North Sea to investigate the migration of the leaking gas into subsurface sand layers using TLFWI.
In this study, an implementation of 3D isotropic elastic TLFWI is applied on short-offset seismic time-lapse field data from the Sleipner area. The baseline dataset is inverted to obtain an elastic baseline model, which is used as initial model in a target-oriented FWI of the monitor dataset. The synthetic datasets generated using the inverted elastic models reproduce the time-lapse field datasets. The inverted baseline and monitor models are used as input in a pre-stack depth migration method to form novel and detailed seismic images of the gas cloud at the Sleipner area. The depth migrated seismic images show that the injected gas might have migrated upwards from the injection point in the bottom part of the Utsira formation through a fault. This paper is organized as follows: In Section 2 the theory behind TLFWI is briefly reviewed. In Section 3 the field datasets used in the inversion are described. The results for TLFWI is given in Section 4 whereas the results for the pre-stack depth migration method is given in Section 5. The results are discussed in Section 6 and the conclusion is given in Section 7.

Time-lapse full waveform inversion
The theory that forms the basis for FWI has been derived many times using different formulations. Here, a brief summary of the theory that is required to understand the inversion workflow is given. The reader is referred to Fichtner (2011) for a complete introduction to the method.
FWI is a classical non-linear inverse problem (Mora, 1987). The goal of the method is to estimate a synthetic model that can be used to compute synthetic data that match field data. FWI is founded on the assumption that elastic waves that propagate in a medium are approximated using a numerical solution of the elastic wave equation. The numerical solution of the wave equation accurately accounts for the complete physics of seismic waves propagating through the medium, and does not suffer from the restriction of, for instance, ray theory. The waveform misfit between the synthetic and the field data is therefore a result from undiscovered features in the synthetic model used in the numerical solution, and not from approximation errors, and can thus be used to estimate an accurate synthetic model of the medium.
In this work, a staggered-grid finite difference method is used to solve the elastic wave equation in an isotropic elastic 3D medium parameterized by the density ( ), the P-wave velocity (v p ) and the S-wave velocity (v s ) (Virieux, 1986;Holberg, 1987;Aki and Richards, 2002). The perfectly matched layer method is used to simulate a non-reflecting unbounded half-space (Berenger, 1994;Zhen et al., 2009), and the free surface is implemented by modifying the upper grid cells in the parameter models as described by Mittet (2002).
A normalized version of the least-squares norm is used to measure the misfit between the synthetic and the field data (Raknes and Arntsen, 2014a). The misfit is minimized using the iterative L-BFGS quasi-Newton optimization method (Nocedal and Wright, 2006), where the adjoint state method is used for computing the model gradients used to update the model in the next iteration (Tarantola, 1984;Mora, 1987;Pratt, 1999).
Inverting for three elastic parameters using FWI is challenging, particularly when working with streamer data (Raknes and Arntsen, 2014b). Therefore, FWI is used to invert for v p whereas and v s are updated using empirical relationships. To reduce the computational cost for performing 3D elastic FWI, the shot sub sampling method is used (Ha and Shin, 2013). The shots and the corresponding receiver traces are divided into four equally sized groups. Data from one group is used in each iteration of the minimization of the misfit. To reduce the computational cost for performing a single shot modeling, the full model is split into local models each which only includes a single shot and the corresponding receivers. The work flow for TLFWI can be divided into two steps: the baseline and the monitor model estimations. The first step consists of estimating the baseline model using the baseline dataset. The inverted baseline model is used as initial model for the inversion of the monitor dataset in the second step. To reduce the possibilities of time-lapse artifacts due to differences in the datasets not related to changes in the subsurface, a target-oriented inversion is performed for the monitor dataset, such that updates in a specific area containing the reservoir and its surroundings are allowed (Ayeni and Biondi, 2010;Zhang and Huang, 2013;Raknes and Arntsen, 2014a). At the end, the time-lapse anomalies are resolved by subtracting the two inverted models. A schematic overview of the time-lapse strategy is given in Fig. 1.

The Sleipner area
The Sleipner area in the North Sea west of Norway consists of the Sleipner Vest and Sleipner Øst field, in addition to some smaller fields (Fig. 2). The gas produced from the Sleipner Vest field contains about 9% CO 2 (Furre and Eiken, 2014). Due to environmental Fig. 1. A schematic overview of the inversion work flow for TLFWI. m0 is the initial model for the baseline inversion, d base is the baseline dataset, mn is the inverted baseline model, dmon is the monitor dataset, and m k is the inverted monitor model. and economical reasons the CO 2 is stripped away from the gas at the Sleipner Øst facility and injected into the Utsira formation at approximately 1000 m depth. The injection of the gas started in 1996 and up to 2010 about 12 million tons of CO 2 were injected into the formation (Furre and Eiken, 2014).
The Utsira formation consists of high porosity sandstones that in total are 200-300 m thick in the Sleipner area. The overburden is approximately 700-800 m and consists primarily of clay-rich sediments . The reservoir properties of the Utsira sand are good with a porosity of 35-40% and permeability of 1-8 Darcys (Eiken et al., 2000;Chadwick et al., 2004). With the assumed reservoir properties in the Utsira formation, the injected CO 2 is near the critical point (Eiken et al., 2000). The gas is injected at a higher pressure and temperature than the surrounding formation, and when the CO 2 gets in contact with the formation it gradually cools down and the density increases to approximately 700 kg/m 3 (Furre and Eiken, 2014). This is below the brine density of approximately 1030 kg/m 3 . The consequence is that when the brine is substituted Fig. 2. A map of the Sleipner area in the North Sea west of Norway. The Sleipner Vest field is a condensate field, and the Sleipner Øst field is a gas field. The other areas on the map are smaller fields. CO2 is stripped of the gas and injected at the Sleipner Øst facility. The gray shaded rectangle is the area where the dataset used in this study was acquired, and it covers the underground formations where the CO2 is injected.
Map courtesy of the Norwegian Petroleum Directorate. by the highly compressible CO 2 , the P-wave velocities in the rock are reduced. For moderate saturations the P-wave velocities may be reduced by as much as 30% . The S-wave velocities are expected to change only with a few percent (Carcione et al., 2006). Hence, the seismic signal from the injected gas should be clearly visible in the seismic data.

Acquisition and preprocessing
To monitor the migration of the gas from the injection point into the Utsira formation, several seismic surveys have been acquired over the area. A pre-injection survey was carried out in 1994, and post-injection surveys were performed in 1999, 2001, 2002and 2012(Furre and Eiken, 2014Furre et al., 2015). In this study the 1994 dataset is used as the baseline dataset, and the 2006 dataset as the monitor dataset. In this time period approximately 8.4 million tonnes of CO 2 were injected into the subsurface .
The baseline dataset was acquired using five receiver cables towed at a depth of 8 m. The cable length was 3000 m and the cross line separation between each cable was 100 m. The source was two 3400 in. 3 air gun arrays towed at 6 m depth. The shot interval was 18.75 m and performed with flip-flop shooting. The recording length for each shot was 5.5 s. During the acquisition of the monitor dataset, approximately the same sail lines were used, in addition to the same cable length and cable separation. The number of cables was eight, and the source was two 3660 in. 3 air gun arrays towed at 6 m depth and fired in flip-flop mode. The shot interval was equal to the interval in the 1994 dataset. The recording length for each shot was 6.0 s.
Before the datasets were released, the contractor applied a standard preprocessing sequence to both datasets. The sequence includes the following steps: 1 Restricted maximum offset to 1700 m, 2 Reduced the recording length to 2.3 s, 3 Applied a signature deconvolution and swell noise filter, 4 Applied a low-cut filter at 6.0 Hz, 5 Sampled the time-step to 2.0 ms, 6 Gained the data using a t 2 scaling factor, 7 Cross-equalization to match the two datasets.
To reduce the computational cost of performing elastic TLFWI, both datasets are decimated. An inversion area, which is smaller than the full acquisition area, is chosen such that it covers the area where the gas was injected. Inside this area, the baseline dataset consists of 852 shots and 570,840 data traces, whereas the monitor dataset consists of 1180 shots and 1,274,400 data traces. The shots in the 1994 dataset are spread over eight sail lines, and the shots in the 2006 dataset are spread over six sail lines (Fig. 3). Unfortunately, there is a gap in the sail lines for the eastern part of the monitor dataset (Fig. 3b). The consequence of this gap, in addition to the different number of receiver cables used in both surveys, is that the data fold is different between the two datasets (Fig. 4). The average data fold for the baseline dataset is 73.8 and for the monitor dataset 84.6. Since the data fold is an indication on how many traces are summed together to form the model gradient, the irregular distribution of the fold between the two datasets might thus contribute to add both acquisition imprints and acquisition differences on the FWI gradient.

Inversion of the field datasets
The inversion work flow that is used for the time-lapse field datasets from Sleipner consists of several steps that in the end yields high resolution elastic parameter models to be used in a depth migration method. The first step is a synthetic example where the implementation of TLFWI is tested with ideal assumptions. This is followed by steps that include regularization of the field datasets into the inversion framework, source wavelet estimation, and inversion of the time-lapse datasets. All steps are described in details in the following.

Synthetic example
Interpretation of the injected gas at Sleipner shows that the gas cloud is separated by small shale layers with a thickness of a few meters (Furre and Eiken, 2014). It is not obvious that the TLFWI method is capable of resolving thin layers in the elastic models. A synthetic model inspired by the geology in the Sleipner area is used to investigate the behavior of the TLFWI method. The synthetic model is built using parts from the SEG/EAGE Overthrust model and information from a well log in the area. The top of the model consists of smoothly varying structures (Fig. 5). A channel system is included at approximately 500 m depth. The bottom part of the model consists of horizontal layers with varying sizes. The target zone for the injection of the gas is the low velocity layer at approximately 1000 m depth.
Two monitor models are used to investigate the resolution of the inversion method. In the first model, the injected gas consists of a single layer that is 75 m high and has a horizontal extent of 3250 m × 1375 m (Fig. 6). In the second model, the injected gas is separated in three layers (Fig. 7). The vertical size of each layer is 12.5 m and the vertical distance between each layer is 12.5 m. As the grid sampling is 12.5 m on all axes in the model, the vertical size of the gas layers is as small as possible in the second model. The horizontal extent of the three layers are different. The injected gas is simulated by reducing v p with 0-300 m/s and with 0-300 kg/m 3 in the areas where the gas is positioned. Since the v s model at Sleipner is expected not to be significantly affected by the gas saturation, it is assumed constant between the baseline and monitor models.
To mimic the field datasets the following parameters for the synthetic datasets are used. For the baseline dataset the survey consists of 632 shots spread over eight sail lines. The shot interval in each sail line is 50 m, and the cross line distance between each sail line is 250 m. Five receiver cables, each with a length of 1700 m, are towed behind the source. Each cable consists of 152 receivers with a spacing of 12.5 m. The cross line distance between each cable is 100 m. The receiver cables and the source are towed at a depth of 12.5 m. The parameters for the monitor dataset are the same as for the baseline dataset, but eight receiver cables are used instead of   five. To simulate the difference in acquisition geometries between the field baseline and monitor datasets, the sail lines for the monitor survey are shifted by 75 m in the cross line direction. With this setup the shots are spread over an area of 4000 m × 2000 m, and covers the injected gas. The source wavelet is a Ricker wavelet with a center frequency of 7.0 Hz, and is equal in the baseline and monitor datasets.
The initial model for the baseline inversion ( Fig. 8) is made by applying a triangular smoothing operator on the true models. The source wavelet is assumed known in the inversion, so that the same wavelet as in the generation of the true data is used in the inversion. The same empirical relationships used to generate the true and v s model are used in the inversion to update the parameters.
The result for the baseline inversion ( Fig. 9) shows that the inversion is able to explain the horizontal layers in the model. The sharp interfaces are smeared out due to the frequencies used in inversion, which yield wavelengths of the same size as the layer thicknesses. The channel system is sharpened by the inversion. The target zone is well resolved, but the magnitude of the elastic parameters is incorrect in the layer due to the smearing effect. At depth the model is hardly updated, which is a consequence of the relative short offsets in the dataset. Footprints of the acquisition geometry are clearly visible in the inverted model. The reflections from the horizontal layers are clearly visible in the true shot gather (Fig. 9d). Since the initial model is smooth, the reflections are not in the shot gather from the initial model (Fig. 9e). The inverted shot gather (Fig. 9f) show that the inversion has introduced the reflections, and the data fit between the true and the inverted shot gathers is considered good.
A comparison of data traces from the two true monitor datasets is given in Fig. 10. The major differences between the traces are in the parts related to the reflection from the gas clouds. For the three layers model, the reflections have smaller amplitudes and oscillate faster than for the one layer model. The differences are increasing as the offset is increased (Fig. 10b). At large offsets the traces are phase shifted. The differences in the data are due to interbed multiples in the three-layer model, as well as interference of the waves that is reflected in the different gas layers.
The result for the inversion of the single layer model is given in Fig. 11, whereas the result for the three layer model is given in Fig. 12. In both examples the inverted time-lapse anomalies are correctly positioned, but are smeared out and affected by the footprints of the acquisition geometries. For the three layer model it is difficult to distinguish the separation of the three layers, and the middle layer is somewhat smeared into the two other layers. The magnitude of the inverted time-lapse anomalies is more or less  correct. Artifacts are visible above and below the time-lapse anomalies. The artifacts are positive anomalies and are more prominent for the three layer model compared to the one layer model. Some edge effects from the regularization are visible. Even though the position of the shots in the baseline and monitor datasets is shifted with respect to each other, the inversion resolves the anomalies without large artifacts.
The true (Fig. 12d) and inverted shot gathers (Fig. 12e) from the three layer model show that the inversion is able to fit the data. The residual shot gathers (Fig. 12f) show several events that are not well explained by the inversion. The "smiling events" following the direct wave are due to coarse shot sampling that results in small oscillations in the sea bottom in the baseline model. In addition, there are some remnants in the reflections from the gas layers. The shot gathers for the one layer model (not included) show the same behavior as the shot gathers for the three layer model.
The two synthetic examples show that elastic FWI is capable of resolving the time-lapse changes directly in the parameter models. The synthetic inversions are performed under ideal assumptions, such that the inversion results are somewhat the best possible outcome from the method. The most important observation is that the inverted time-lapse anomalies are smeared out due to the relative low frequencies used in the inversion. In addition, it is important to be aware of the time-lapse artifacts introduced by the inversion.

Data regularization
Data regularization of the field time-lapse datasets is a crucial step in TLFWI since FWI is a data driven method, and errors made at this point may have great impact on the final results. The regularization process consists of choosing the frequencies to be used in the inversion, moving the data traces from the continuous grid to the numerical grid, and lastly reverse the preprocessing sequence applied by the contractor to obtain field data as close to the raw data as possible.
The solution space for the FWI problem consists of local minima and a global minimum (Fichtner, 2011). The goal for FWI is convergence to the global minimum, or more likely to the neighborhood of the global minimum. Since the inverse problem is solved using a local optimization method, it may converge to a local minimum far away from the global minimum. To prevent this from happening it is important to use low frequencies in the inversion (Virieux and Operto, 2009), since at low scales the number of local minima is greatly reduced. Based on this fact, Bunks et al. (1995) suggested to perform sequential inversion runs starting at low scales and gradually increase the scales by introducing higher frequencies in the data in the following inversion runs. Sirgue and Pratt (2004) introduced a method to choose frequency bands to use with the approach introduced by Bunks et al. (1995). Using their method, the following three frequency bands are used in the inversion: 6.0-8.0 Hz, 6.0-11.0 Hz, and 6.0-15.0 Hz. The starting model for each frequency band is the final model from the previous frequency band.
The numerical grid used in the inversion is discrete and regular. To move the receiver traces in both datasets from the continuous and irregular real world grid, a simple nearest neighbor algorithm is applied. For each receiver the distances to all four grid corners in the horizontal plane are computed. The receiver is then moved, without phase or amplitude corrections, to the closest grid corner.
The only preprocessing step that is reversible is the time gain applied to both datasets. To reverse the gain, each trace is multiplied by t −2 . Interpolation is used to reduce the time sampling to 1.0 ms, which is the time sampling used in the modeling. The last step is to add time delay to the data, such that the first arrivals and the amplitude decay in the field datasets match the synthetic counterpart.

Source wavelet estimation
To estimate the source wavelet, the near-offset traces from several shots are stacked together and then everything but the direct arrival is muted away. The resulting dataset is used as input to FWI, where the source wavelet is inverted for. Since several frequency  bands are used in the inversion, a source wavelet is estimated for each frequency band.
The inverted source wavelets for both field datasets are given in Fig. 13. As the frequency band is increased the source wavelets are sharpened.

Inversion results
Before the inversion is started the empirical relationships between v p , and v s must be decided. The choice of empirical relationships to couple the v p model to the v s and models is not straightforward since it is difficult to find a perfect match for all rock types (Mavko et al., 2009). Based on the knowledge of the rock types in the area, the well-known Gardner's relationship (Gardner et al., 1974), which is valid for many rock types, is used to link to v p . This relationship is given as The relationship between v s and v p is governed by the so-called "mud-rock" line (Castagna et al., 1985), given as v s = 0.862v p − 1172. (2) In the above equations [ ] = m/kg 3 and [ The inverted baseline model after the inversion of the three frequency bands is shown in Figs. 14a and 15a. The initial model was made using a conventional tomography method, and thus it is a smooth model without clear structure details. Already using the first frequency band, structure details starts to appear in the inverted model. At this stage the high velocity layer at approximately 600 m depth is the most visible one, whereas shadows of the deeper layers are visible. As the data frequencies are increased in the inversion, more details show up, particular at the deepest parts of the model. In the final inverted baseline model several clear layers are visible. The most emphasized layers are at approximately 600 m, 900 m and 1200 m depths. The results from the baseline inversion are discussed in details by Raknes et al. (2015).
The size of the gas cloud at Sleipner is, with some uncertainties, known (Cavanagh and Haszeldine, 2014;Furre and Eiken, 2014). This information is used to create the target-oriented model constraint for the monitor inversion. The boundaries for the constraint are made wider than the gas cloud to allow the inversion to explain potential developments of the gas cloud that are not already known. In addition they are smoothed to avoid sharp edge effects during the model updates. The total size of the model constraint is 375 m × 1975 m × 3500 m.
In the inverted monitor model (Figs. 14b and 15b) both low and high velocity layers have been introduced in the model inside the constraint area. In the cross line direction (Fig. 15b) four low  velocity layers are visible. The low velocity layers are surrounded above and below by high velocity layers. In the in line direction (Fig. 14b) the same type of events are visible, but here the layers are not as clearly visible as in the cross line direction. The resolved time-lapse anomalies are clearly visible in time-lapse models (Figs. 14c and 15c). In the southern part there is a velocity increase in the upper part of the constraint area, and a velocity decrease in the lower part. In the northern part it is opposite. Strong positive anomalies are visible on the west and east sides of the constraint area. Fig. 16a and b show shot gathers from a single cable from the field baseline and monitor datasets, respectively. The figures show the data as they are used in the inversion, that is, after filtering and data regularization. It can be seen immediately that the baseline field data are affected by noise and artifacts from the data acquisition. Vertical stripes are visible which can be due to swell noise and  feathering of the streamers. Several reflected events are visible, but they are distorted by linear events that result in interference patterns in the reflections. The reflections from the gas cloud are clearly visible in the monitor data (marked with arrows in the figure).
Shot gathers from the initial models are given in Fig. 16c and d, whereas shot gathers from the final inverted models are given in Fig. 16e and f. The inversion is able to introduce reflections and position them correctly in depth. The refracted wave is modeled correctly, though with amplitude differences. For the monitor data, the reflections from the gas cloud are positioned correctly in depth. The major differences between the synthetic shot gathers and the field shot gathers are the presence of the linear events below the direct wave and the interference pattern. Despite these differences, the data fit is better for the synthetic data computed with the inverted models ( Fig. 16e and f), than for the data computed with the initial models ( Fig. 16c and d).

Pre-stack depth migration
The inverted baseline and monitor models are used as input in a conventional one-way shot-profile imaging method (Etgen et al., 2009) to create depth migrated seismic image cubes of the Sleipner area. Two seismic image cubes are generated. In the first migration the baseline dataset is migrated using the inverted baseline model (Figs. 14a and 15a). This baseline cube yields an image of the area before the injection started. In the second migration the monitor dataset is migrated using the inverted monitor model (Figs. 14b and 15b), to form the seismic monitor image cube after ten years of injection. Fig. 17 shows vertical slices through the baseline and monitor seismic image cubes. The two images are similar in the overburden, with small differences that can be explained by differences in the datasets. The interface at approximately 225 m depth includes some v-shaped events that can be related to a channel system at these depths (Raknes et al., 2015). A horizontal layer is visible at approximately 400 m depth. The velocities increase slightly in this layer in the velocity model compared to the overburden (Fig. 14a). The layer at 600 m depth includes unfocused events. Thus the inversion has not explained the velocity model well in this area. The top of the Utsira formation is clearly visible at approximately 800 m depth, followed by several layers at larger depths (Fig. 17a). In general the baseline seismic image is focused (Fig. 17a). The injected gas is clearly visible in the monitor seismic image (Fig. 17b). The gas is injected at approximately 1000 m depth, and from the seismic image it is clear that the gas has migrated into different layers as it has ascended from the injection point. Due to the seismic response the gas acts as a contrast fluid and thus emphasizes the structures it is injected into. The extent of the gas is different in each of the layers. The gas has evolved longer horizontally in the regions close to the injection point, than in the regions closer to the top of the Utsira formation.
Close-up seismic images of the Utsira formation and the area where the gas is injected are given in Fig. 18. In the middle of the images (approximately at y = 4100 m) a slightly slanted vertical event is visible in both the baseline (Fig. 18a) and the monitor (Fig. 18b) images. In the monitor image the event is enhanced compared with the baseline image. In addition, oscillating events not apparent in the baseline image are visible from approximately 1000 m depth and up to approximately 800 m depth in the monitor image.
In Fig. 19 horizontal slices through the top of the Utsira formation are shown. Several events (marked with arrows in the seismic images) are similar in the baseline and monitor images, and thus to get comparable seismic images, the magnitude of the black anomaly at approximately (x, y) = (900 m, 5100 m) is used as clip value in each seismic image. These clip values are used for the horizontal slices to be shown in the following. The seismic response of the gas cloud is clearly visible in the monitor seismic image. The major part of the gas is clustered together in the middle of the image, and several black circular dots are visible around the gas cloud. North and south of the gas cloud the same structure is visible in the two seismic images. The velocity is decreased compared to the surroundings in the areas where the gas is positioned.
In the horizontal slices at 881.25 m depth (Fig. 20) the gas cloud is bigger than at the top region. The cloud is consisting of two smaller clouds where the northernmost is the largest. In the southern part of both images a large magnitude anomaly is visible. The southern cloud follows the structure that is visible in the baseline model. The northern cloud stops at the east-west structure (at approximately y = 4000 m) visible in the baseline model. The v p model shows a decrease in velocities in the northeast region that correlates with the monitor seismic event. The same is true for the event in the southwest part of the slice. An interesting feature in the image is the circular black event positioned approximately in the middle of the baseline seismic image (marked with an arrow in the slices). The seismic response has increased considerably around the point in the monitor seismic image compared to the baseline seismic image. In general, the shape of the gas cloud is in agreement with the structure of the v p model.
The horizontal slices at 918.75 m depth (Fig. 21) show that the gas has migrated approximately 3.0 km in the north-south direction and 1.2 km in the west-east direction. By comparing the baseline and monitor seismic images, it is clear that the gas has migrated into a structure that is visible on the baseline seismic image. Fig. 22 shows horizontal slices in the bottom region of the Utsira formation. The seismic response from the gas cloud is clearly visible in the monitor image. By considering the overlay plot of the seismic monitor image and the monitor velocity model, several correlations are visible. The event in the northeast corresponds to a low velocity zone, whereas the middle region is a mixture of both high and low To demonstrate the quality of the monitor migration model, the monitor dataset is migrated using the inverted baseline model. The results of the migration is given in Fig. 23. By comparing these results with the already discussed results  it is clear that the latter results include better focused events and in general more details. This is clearly visible in the deepest horizontal slices, where the former horizontal slice (Fig. 23d) is severely smoothed compared to the latter counterpart (Fig. 22b). This can be explained by the simple fact the gas has been taken into account in the monitor migration model, and not in the baseline model.
The comparison of the common image point angle gathers (Fig. 24) between the migration of the monitor dataset using the baseline and the monitor model, show that for the monitor model the angle gathers are flatter than for the baseline model in the regions of the model constraint. This is an indication that the kinematics for the P-wave reflections are better explained in the monitor model than the baseline model. The small differences in the other areas are a consequence of differences in the datasets, and not the migration model since this is identical in these areas.

Discussion
The use of empirical relationships to update the parameters not inverted for and its impact on the inversion results are questionable. Raknes and Arntsen (2014a) performed a synthetic time-lapse sensitivity analysis using empirical relationships to update these parameters. They concluded that if wrong empirical relationships were used in the model updates, the resolution of the corresponding inverted time-lapse anomaly decreased. However, even though wrong empirical relationships were used to update both and v s , the inversion was able to resolve the time-lapse anomaly. In the preliminary work with the field datasets, and v s were estimated without success using FWI. The reasons for the failure are lack of Swave information in the data (since it is acquired in the water layer), in addition to the fact that is difficult to reconstruct using FWI (Virieux and Operto, 2009;Raknes and Arntsen, 2014b). The idea of using empirical relationships to couple the elastic parameters not inverted for is a way to circumvent the inversion problem in cases where estimates for the three elastic parameters are wanted.
One can question the validity of the empirical relationships used in this study (Eqs. (1) and (2)). The major issue with empirical relationships, in general, is that they are restricted to certain rock types. In this case, the choice of relationships was made on knowledge of the rock types in the area and the validity of the chosen relationships (Gardner et al., 1974;Castagna et al., 1985). An alternative to use pre-defined relationships is to estimate specific relationships using well log information from the area. The quality and the amount of well log data are, however, not enough for estimating reliable relationships. Another alternative is to use different relationships for different parts of the subsurface. This would yield a challenge in how to discriminate the relationships to the different regions of the model, particularly in the FWI framework where structures can be considerably shifted in space during the inversion. The FWI problem is a non-linear and ill-posed problem. Thus, a small change in the assumptions may have a severe influence on the results because the problem is solved using a local optimization method. For TLFWI applications this is a major issue since the inversion of the baseline and monitor dataset may lead to convergence to models that in the end give artifacts or false time-lapse anomalies in the parameter models. One could use a global optimization method in order to find the global solution. For the problem solved here it is impossible in practice, due to the extreme number of forward modelings required to arrive at the global solution. To circumvent this problem a target-oriented approach was used to update only certain parts of the model (Ayeni and Biondi, 2010;Zhang and Huang, 2013;Raknes and Arntsen, 2014a). It is important to have in mind that putting constraints on the monitor inversion can be problematic, since changes in other parts of the model are not updated during the monitor inversion. Furthermore, if the baseline model convergence is not good enough, then the constraints may yield time-lapse results that are a result of bad baseline convergence and not true timelapse changes. If, on the other hand, the position of the time-lapse area of interest is well-known, then using model constraints to focus the seismic energy and thus reduce the solution space for the inverse problem is a good option to obtain satisfactory inversion results, as shown by the synthetic examples.
The regularization of the field datasets to fit into the FWI work flow is essential to obtain reliable inversion results. Here, a nearest neighbor technique without amplitude or phase corrections was used for moving the receiver traces to a numerical grid point. This could potentially be problematic since FWI is sensitive to changes in both amplitude and phase. However, with the numerical grid used in the inversion, a receiver is moved a small fraction of the wavelengths in the data at the receiver locations. Due to coarse and irregular trace sampling it is not clear that more accurate interpolations schemes (Choi and Munson, 1998;Hindriks and Duijndam, 2000;Özdemir et al., 2010) could improve the data. The most questionable factor is the preprocessing sequence applied to the datasets before they were released. In an ideal setting the field data used in FWI should be raw and unprocessed data. The result of the preprocessing steps is that it can be problematic to match the field and synthetic data, and that important information is removed from the data. The consequence of the restriction of offsets and the reduction in recording length is that wide-angle data are effectively removed from the dataset. Furthermore, the swell-noise and the low-cut frequency filters remove important signals from the data.
It is a well-known fact that long-offset data as well as low frequent signals are important to resolve the subsurface, in addition to preventing the inversion from running into a local minimum (Virieux and Operto, 2009).
The inverted elastic models (Figs. 11,12,14 and 15) show that TLFWI is capable of resolving time-lapse anomalies in the elastic parameters. By a comparison of the synthetic inverted gas clouds (Figs. 11 and 12) and the field data gas cloud (Figs. 14c and 15c) it becomes clear that the inverted time-lapse models are influenced by artifacts. The artifacts in the synthetic and field time-lapse models are visible as positive anomalies above, in between, and below the time-lapse anomalies. They are more prominent in the field data example, which is expected due to the uncertainties in the data. The synthetic and field data artifacts are, however, to some extent similar. The high velocity layers in between the low velocity layers in the field models can be interpreted as side lobes to the inversion due to the frequency content of the field data used in the inversion. The same type of artifacts are confirmed by other studies (Zhang and Huang, 2013;Raknes and Arntsen, 2014a;Asnaashari et al., 2015). Since there are no physical explanations for a strong P-wave velocity increase due to CO 2 injection, the strong positive time-lapse anomalies on the east and west sides of the model constraint (Fig. 15c) can be explained by the difference in acquisition geometries between the baseline and monitor surveys. The monitor dataset covers a larger area than the baseline dataset (Fig. 3) and has a different data fold (Fig. 4). The data quality of the monitor data is, in addition, better due to better equipments used in the acquisition of the monitor data. Since the baseline dataset did not resolve the areas at the boundaries of the constraint satisfactorily, the monitor inversion is trying to explain the model in the edge areas. The consequence is that these events show up as false time-lapse anomalies in the models. Despite these shortcomings, the synthetic and the field example demonstrate the potential of using TLFWI to resolve time-lapse changes directly in the elastic parameters.
There is not a perfect match between the synthetic and the field dataset (Fig. 16). The inversion is able to reproduce the reflections from the sediments and the gas cloud, but the linear events and interference patterns below the direct arrival are not well explained. These events can be explained by converted wave energy due to layers close to the sea bottom. To be more precise, if there is a decrease with a sufficient contrast in v p and v s in a layer close to the sea bottom, then the downgoing P-wave is reflected back as converted P-wave when recorded by the receivers. Such a converted wave gives the same type of events visible in the field datasets ( Fig. 16a and b). For the field data, several layers are visible close to the sea bottom (Fig. 17). Moreover, the inversion has introduced layers close to the sea bottom with a decrease in v p (Figs. 14 and 15). The contrasts in v p and thus v s are, however, not high enough to generate a converted wave that yield the same type of events in the synthetic data as in the field data. The reasons are that the sea bottom, in general, is difficult to estimate correctly with FWI, and that the Castagna relation (Eq. (2)) has been used in the whole model. For the uppermost sediments, where the v p /v s ratio can be as high as 5.0, the relation is not valid. However, the inverted monitor model produces flatter common image point angle gathers than the baseline model (Fig. 24). In addition, there are correlations between the depth migrated images and the elastic monitor model (Figs. 19-22). The estimated velocities in the baseline and monitor models are within the interval that is expected (Eiken et al., 2000;Arts et al., 2004;Carcione et al., 2006). Furthermore, the decrease in the estimated v p velocities in the Utsira formation after CO 2 injection is of the same magnitude as expected . The depth migrated seismic images using the inverted monitor model are well focused, whereas the corresponding seismic images created using the inverted baseline model lack focusing. Hence, the inverted elastic monitor model is closer to the true solution after gas injection than the inverted elastic baseline model.
From the results there are no doubts that the FWI produces an elastic model that is more detailed than a conventional tomography model. The details in the corresponding seismic depth migrated images  are a result of the high frequencies used in the migration, in addition to the details in the migration velocity model. The monitor migration using the baseline model (Fig. 23) demonstrates that high frequencies and a correct overburden are not sufficient for producing detailed seismic images of the gas cloud. A correct migration model inside the gas cloud is crucial to achieve high resolution seismic images. This fact is also demonstrated by Warner et al. (2013) that used FWI to estimate an anisotropic acoustic migration model at the Tommeliten field that yielded a high resolution seismic image of an underburden gas cloud.
It is worth mentioning that the results presented here have not been used for any quantitative interpretation to estimate for instance CO 2 saturation in the Utsira formation. Several such studies exist in the literature, see for instance Chadwick et al. (2009), Chadwick et al. (2010 and Ghosh et al. (2015). Queißer and Singh (2013) performed 2D elastic FWI using a single line from the Sleipner area. Compared to the work flow presented here, they used 2D modeling and inversion, a finer numerical grid and 38Hz as the maximum dominant frequency. Their inverted baseline velocity model (Queißer and Singh, 2013, Fig. 6) includes clearly resolved layers at 300 m and 600 m depth, and less clearly resolved layers at 1000 m and 1100-1200 m depths. Their velocity model is in general patchy. The vertical slices through the 3D baseline velocity model (Figs. 14a and 15a) show clearly resolved layers down to approximately 1200 m depth. The layers at 300 m and 600 m depths correlates with the layers in the 2D model. From approximately 800 m depth, the 2D and 3D models do not correlate well, mainly due to the patchiness of the 2D model. In the inverted 2D time-lapse model (Queißer and Singh, 2013, Fig. 10) several gas layers are visible in the Utsira formation. These layers are mainly effected by a reduction in v p , but some increases in v p are also visible. The 3D time-lapse model (Figs. 14c and 15c) has also several layers with reduction in v p . The resolution of the layers in the 3D model is not as in the 2D model mainly due to the higher frequencies used to estimate the latter model. It is difficult to fairly compare the 2D and 3D time-lapse inversion results because of the differences in the way the inversions have been performed. However, in general, it is expected that the 3D inversion results are more reliable than the 2D inversion results as the physics of the elastic wave propagation is better accounted for in a full 3D inversion framework compared to a 2D framework.
The gas cloud at Sleipner has been interpreted several times using conventional time migrated seismic cubes (Eiken et al., 2000;Arts et al., 2004Arts et al., , 2008Chadwick et al., 2010). Even though the conventional interpretations of the gas cloud are made using different methods, it is interesting to compare the conventional results with the depth migrated results (Figs. 19-22). The horizontal extent of the gas cloud has been interpreted to be approximately 3.6 km × 1.0 km (Chadwick et al., 2009). By interpreting the seismic images the horizontal extent is found to be approximately 3.2 km × 1.0 km. The conventional interpretations show that the gas cloud has migrated following a north-south trend (Bickle et al., 2007;Chadwick et al., 2010, Fig . 1). The depth-migrated images  show the same trend. Another well-known feature in the conventional interpretation is the so-called fingering of the gas (Chadwick et al., 2009). The fingering feature is visible in the very top of the gas cloud (not shown here) and to some extent in the depth migrated image in Fig. 21.
The pathways that the gas has migrated from the injection point and up to the top of the Utsira have been discussed in the literature (e.g. Eiken et al., 2000;Arts et al., 2008;Cavanagh and Haszeldine, 2014;Furre and Eiken, 2014). The vertical slices (Figs. 17 and 18) confirm that the injected gas has migrated upwards from the injection point and to the top of the Utsira formation. The vertical events penetrating the gas cloud (Fig. 18b) may be interpreted as socalled gas chimneys, that are vertical gas-filled "pipelines" through the formation. Arntsen et al. (2007) did a synthetic study of gas chimneys and their results are to some extent comparable with the events visible in the monitor images. However, the reflector at approximately 925 m depth in the baseline image (Fig. 18a) shows a vertical discontinuity at approximately y = 4100 m, which might indicate that this is a fault. Moreover, a closer examination of the horizontal slice at 918.75 m depth (Fig. 21), an east-west trending event is clearly visible (middle arrow in Fig. 21) in the monitor image, and to some extent in the baseline image. It is worth noting that the horizontal position of this event correlates with the vertical event (at x = 1400 m) in Fig. 18b. The same is the case for the horizontal slice at 943.75 m depth (Fig. 22). In the baseline seismic image at 881.25 m depth (Fig. 20), the seismic response is severely changed in the areas around the circular event (marked with an arrow) between the baseline and monitor images. In addition, an east-west trending event is going from the same position. This event is also, to some extent, visible on the baseline image.
Hence, the vertical event that penetrates the gas cloud may be interpreted as a fault in the Utsira formation. When the gas was injected, the fault might have acted as a migration path for the gas as it ascended through the formation. The fault can be interpreted to be, at least, 500-800 m long from the horizontal slices. Furthermore, the seismic images show that the gas has migrated horizontally and been trapped by stratigraphic traps in the structures.

Conclusion
This paper has demonstrated an application of a three dimensional isotropic elastic TLFWI implementation using time-lapse seismic data from the Sleipner area, offshore Norway. The results show that TLFWI is capable of creating detailed elastic parameter models that reproduces the time-lapse field datasets. These models are used in a pre-stack depth migration method to create seismic images of the area before and after injection of CO 2 . The seismic images show events and discontinuities that can explain the migration pathways of the injected gas in the Utsira Formation.