Anisotropic anelastic seismic full waveform modeling and inversion : Application to North Sea offset VSP data

In sedimentary basin, elastic anisotropy can be described by a transverse isotropic medium and the attenuation in the seismic bandwidth can be approximated by a quasi-constant quality factor. Very few full waveform inversions were conducted for such realistic media to show the feasibility and the benefit of this approach. For illustration, we have chosen two offset VSP datasets from the North Sea displaying attenuated phases and where the medium is known to be transverse isotropic. By inverting elastic parameters, anisotropy, shear attenuation and source functions, we have been able to find an Earth model reproducing fairly the real data. By exploiting all the information in the seismograms, full waveform inversion allows us to localize and characterize the Brent gas reservoir target.


Introduction
Anisotropy is often observed due to the thin layering or aligned micro-structures, like small fractures.This type of anisotropy can be described by vertical (VTI) and horizontal transverse isotropy (HTI) when, respectively, layering is horizontal and the cracks are vertical.Moreover, the elastic approximation of Earth properties for seismic wave propagation is limited as waves undergo attenuation and dispersion that can be described by a quasi constant quality factor over the seismic bandwidth.Anelastic and anisotropic parameters need to be taken into account for a consistent seismic full waveform inverse problem.If the redundancy in the seismic surface data allows us to neglect many wave phases in the seismograms such as P to S wave conversions, the scarcity of data in borehole seismic, especially for offset VSP (OVSP), challenges us to interpret the full content of the seismograms.For that reason, the complexity of the wavefield (reflections, transmissions, phase conversions, etc. seen in the seismograms) needs to be accurately modeled.According to the review done by Virieux and Operto (2009) about full waveform inversion, the reconstruction of anisotropic parameters is probably one of the most undeveloped and challenging fields of investigation and they conclude that incorporating more sophisticated wave phenomena (attenuation, elasticity, anisotropy) in modeling and inversion is another field of investigation to be addressed.On one hand, the full waveform inversion for visco-elastic parameters in time domain have been proposed by Charara et al. (2000) and Barnes et al. (2004).The feasibility and practicality of this approach has been demonstrated on several noise free synthetic data: for a 1D offset VSP numerical inversion experiment (Barnes et al., 2004) and for a 2D crosswell numerical inversion experiment (Charara et al., 2004) and successfully tested on real data (Barnes and Charara, 2009).In addition, a feasibility study of the full waveform inversion for VTI parameters for a synthetic seismic crosswell data shows a reliable reconstruction of the Earth model even with noisy data (Barnes et al., 2008).
Based on these previous studies, we investigate the benefit of simultaneous full waveform inversion of elastic parameters, VTI Thomsen parameters and shear Q-factor when applied to two North Sea OVSP datasets.

Visco-anisotropic modeling in time domain
The linear visco-elasticity constitutive law, for a quasi constant quality factor, is modeled by the superposition of relaxation mechanisms, classically called Zener or standard visco-elastic bodies.The key concept of this technique is the replacement of the time convolution between relaxation rates and strain by a set of first order temporal partial differential equations, storing the strain history interactions with the medium trough new fields called strain memory variables (Carcione, 1990).The full set of visco-elastodynamic equations using this constitutive relation can be solved numerical by the finite difference method.Furthermore, the seismic modeling is based on the cylindrical system of coordinates by assuming azimuthal invariance of the propagation fields.The 3D Green function can be modeled as a 2D numerical problem by assuming axi-symmetry of the medium and the wavefields (Igel et al., 1996).The wave propagation equation is discretized using a 4th order staggered grid finite difference scheme for the cylindrical system of coordinates.This scheme allows one to model VTI anisotropy.

The inversion method and procedure
The inverse problem, solved by local optimization methods, can be expressed as the minimization of a misfit function (Tarantola, 1987).For the case of least squares, the misfit function is a scalar function defined over the model space as ( ) The minimization of the misfit function is performed using a conjugate gradient method.The iterative process of a non linear inversion is developed by Tarantola (1987); the expressions of the Fréchet derivatives for the elastic parameters and for the source function can be found in Tarantola (1986) and for the visco-elastic parameters in Charara et al. (2000).
In the present inversion, we invert for the P and S-wave velocities, the density, the Thomsen parameters ε and δ for VTI anisotropy, the shear quality factor Qs and the source time functions.All the parameter fields are 2D except for the Qs parameter which is considered 1D at this stage of the study for stabilization purpose.The model covariance matrix is filled with independent horizontal and vertical Laplacian correlations (Charara et al., 1996).The inversion procedure starts from the prior model: a stratified model inferred from travel time inversion and well logs.Then, four inversions are performed successively.From one inversion step to the next, the spatial correlations decrease while the frequency content of inverted data is increased.

North Sea Offset VSP data
The Lille-Frigg field is located on the eastern margin of the Viking Graben, about 22 km north-east of the Frigg field.The structure consists of a narrow N-S elongated horst with gas trapped in the Brent reservoir.The Brent formation is known to be very variable in thickness and facies.The well 25/2-C1H drilled in Lille-Frigg field encountered a fault dipping east and only the lower part of the Brent formation (Minsaas et al., 1994).
Two three component OVSP datasets were acquired with the aim of better defining the structural position of the fault plane and to ascertain the distance of the whole Brent formation from the well (Muller and Ediriweera, 1993).Based on the interpretation of these OVSPs, of logging data, high resolution borehole imagery and other seismic data (walkaways and other OVSPs), a side track 200 m west of the well was drilled finding the total sequence of the Brent formation (Minsaas et al., 1994).
The acquisition geometries of the two OVSPs are shown in Figure 1, they are denoted "West" and "East" and shots are aligned with the well.The well is vertical, the source offset is 2 km for both shots and the data were collected at depth levels ranging from 2300 m down to 3950 m.

Inversion results and discussion
The presented inversion results are for the 4th inversion (20Hz).The synthetic data obtained at the convergence are shown in Figure 2 as well as the corresponding observed data and residuals for the west OVSP.The fit is fair; the remaining misfit is 20%.Moreover, the major part of the residuals is unstructured noise.
The inverted source time function for the West offset VSP is shown in Figure 3.The high frequency content of the source is increasing during the inversion process as the frequency content in the observed data has been increased from the previous inversion step.
The field representing the Vs/Vp ratio is plotted in Figure 1 and is a good indicator of the gas bearing reservoir zone.The other estimated parameter fields are zoomed in Figure 4 displaying complementary structural and petrophysical information.The geological interpretation provided by (Minsaas et al., 1994) is plotted on the same figure as well as the side track of the well.
The inverted attenuation parameters are very smooth but show strong attenuation in the first layers below the water bottom necessary to explain the attenuated converted S-wave from the sea floor.The gas reservoir attenuation is not resolved due to the 1D inversion for the attenuation parameter.The anisotropic parameters cannot be well resolved due to the limited number of shots and therefore the limited range of incidence angles.The Thomsen ε and δ estimated fields provide poor information in absolute value (ε is shown in Figure 4) but some anomalies seems to be associated to real anisotropic layers.The inverted elastic parameters are of main importance in order to identify and delimitate the gas bearing reservoir as shown in Figure 4 (Vs/ Vp ratio and the Poisson ratio fields).However, the west/ east extension of the resolved region in the images is limited to the well vicinity (from 200 m to 500 m).This is due to the narrow illumination region both for transmitted (down-going) waves and reflected (up-going) waves although multiples or converted waves are inverted.Moreover, the starting model being stratified, the resolved region does not appear clearly.
The next stage is to define a better, likely 2D, prior model, to better understand the coupling of the shear attenuation inversion with source inversion, to introduce a compressive attenuation parameter and finally to provide a posterior resolution map for each estimated parameter field.

Conclusions
The simultaneous inversion of elastic, VTI Thomsen parameters and shear Q-factor allows us to extract more  information from the data.By using the adequate Earth properties, it decreases the modeling "noise" in the inversion process, even if the constraints over certain parameters are weak (and consequently providing a poor resolution).The complexity of the recorded wavefield is fairly reproduced by the synthetic data.In the resolved parts of the model, the results are consistent with the currently admitted structural interpretation, part of the structure being well retrieved.The present anisotropic visco-elastic full-wave inversion of offset VSPs data from the North Sea illustrates the impact of multiparameter inversion when using the adequate Earth properties; it shows the feasibility and the benefit of our proposed method.
where symbols are defined as follow: On the model domain, m is a model, pr ∆ = m m m is the difference between the model m and the a priori model m pr and C M denotes the model space covariance matrix; On the data domain, , i.e., the difference between the observed data d obs and the synthetic data g(m) obtained by simulation of the wave propagation; and, C D denotes the covariance matrix over the data space.The misfit function measures the discrepancy between observed and synthetic data, and also between the model and the prior model.The reference misfit used to normalize the current misfit is obtained for obs ∆ = d d and ∆ = m 0.

Figure 1 .
Figure 1.S-wave velocity over P-wave velocity ratio field.The acquisition geometries of the West and East OVSP are shown: stars near the surface denote the sources location while down triangles denote the receiver locations.The resolved region depends on the wave: the direct P-wave illuminates a triangular region from the source to the antenna, the down-going converted P-S at the antenna depth illuminates a region close to the well, the reflected waves illuminate region bellow the bottom of the well.The Vs/Vp ratio allows to clearly distinguish the gas bearing reservoir from other formations (purple region).

Figure 2 .
Figure 2. Horizontal component of the West OVSP data.Low-pass filtered observed data (corner frequencies of 20 and 40Hz), synthetic data and residuals are displayed respectively at left, on the middle and at right.The fit between observed and synthetic is fairly good providing low residuals (remaining misfit is 20%).Structured energy in the residuals shows that some reflected and transmitted S-waves are not fully explained.

Figure 3 .
Figure 3. Source time function for the 4 th inversion of the West OVSP data.The starting wavelet is denoted by the thick red line while the thin orange, green, blue and black lines denote respectively the inverted wavelet at iteration 5, 20, 50 and 100.

Figure 4 .
Figure 4. Images of the estimated elastic parameter fields in the well neighborhood and below the well.The down triangle denotes the receivers (1 over 2 are plotted).The side track is indicated by the thick magenta dashed line on the Poisson ratio image.The geological interpretation(Minsaas et al., 1994;Mittet et al., 1997) is indicated by thin lines: the black line denotes the major normal fault cutting the Brent formation, the dashed black line is a likely antithetic fault, the white line is the BCU (base cretaceous unconformity), the dotted black line is the coal seismic marker in the lower part of the Brent formation and the dark red solid line is the Brent bottom.The estimated fields provide complementary information on the structure, e.g. the density field allows to localize precisely the coal marker while the short wavelengths of the S-wave velocity field helps identifying the major fault.The Vs/Vp ratio or the Poisson ratio are good indicators of the presence of gas in the reservoir.According to these fields, the gas bearing region corresponds to the upper part of the Brent formation (in purple for Vs/Vp ratio and in red for Poisson ratio, at left of the fault).The structure bellow the well is difficult to interpret.Finally, the anisotropic parameters like the Thomsen ε, are not well resolved (narrow incidence angle range), only the zones illuminated by the down-going waves are reliable, for instance, the anisotropic layers located near 2980 m and 3100 m.