Estimation of Seismic Wave Attenuation from 3D Seismic Data: A Case Study of OBC Data Acquired in an Offshore Oilﬁeld

: Previous studies performed in Abu Dhabi oilﬁelds, United Arab Emirates, revealed the direct link of seismic wave attenuation to petrophysical properties of rocks. However, all those studies were based on zero offset VSP data, which limits the attenuation estimation at one location only. This is due to the difﬁculty of estimating attenuation from 3D seismic data, especially in carbonate rocks. To overcome this difﬁculty, we developed a workﬂow based on the centroid frequency shift method and Gabor transform which is optimized by using VSP data. The workﬂow was applied on 3D Ocean Bottom Cable seismic data. Distinct attenuation anomalies were observed in highly heterogeneous and saturated zones, such as the reservoirs and aquifers. Scattering shows signiﬁcant contribution in attenuation anomalies, which is unusual in sandstones. This is due to the complex texture and heterogeneous nature of carbonate rocks. Furthermore, attenuation mechanisms such as frictional relative movement between ﬂuids and solid grains, are most likely other important causes of attenuation anomalies. The slight lateral variation of attenuation reﬂects the lateral homogeneous stratigraphy of the oilﬁeld. The results demonstrate the potential of seismic wave attenuation for delineating heterogeneous zones with high ﬂuid content, which can substantially help for enhancing oil recovery. technique due to the location of OBC receivers on the seaﬂoor, was overcome by applying an optimized Gabor transform to create a collection of spectrum traces in the time domain, followed by a conversion from time to depth domain. The results showed high attenuation anomalies in highly saturated and heterogeneous zones, such as the reservoir zones and Simsima Formation. Such results show the advantage of using seismic wave attenuation for petroleum exploration and reservoir characterization. A slight lateral variation at each seismic line was revealed by the 2D attenuation proﬁles, which agrees with velocity proﬁles and results of full waveform inversion.


Introduction
Enhanced oil recovery (EOR) aims to improve oil and gas production by facilitating the flow of fluids. The three major EOR techniques are chemical flooding [1], miscible displacement [2] and thermal recovery [2]. The optimal application of each technique depends on petrophysical properties of reservoirs such as permeability, water saturations, porosity and fluid properties. Fracture network, control the hydraulic and many reservoir petrophysical properties. Therefore, receiving fractures properties is crucial for applying EOR techniques. Conventional methods used for investigating naturally fractured reservoirs are generally based on petro-physical well-logs, such as sonic, gamma ray, and resistivity logs or on core description. However, such data are obtained in the vicinity of the well location and at small scale. At large scale, fracture properties are obtained from surface seismic data by applying stiffness anisotropy and geometry-based methods. Chichinina and Martyushev [3] compared Thomson's anisotropy parameters in dry and saturated cracked rocks. They found out that in dry rock γ is much smaller than ε, while δ is larger than ε and becomes negative when the rock is fully saturated. By looking at the dominant direction of the extracted maximum curvature, Ali, et al. [4] successfully defined the dominant orientation of fractures in a reservoir of an offshore Abu Dhabi oilfield, and found that it is parallel to the NNE-SSW and NE-SW directions. However, such methods are sensitive to any geological feature causing seismic anisotropy other than fractures, such as lithological alignment.
Theoretical [5], experimental [6], and field studies [7] have demonstrated the potential and superiority of seismic wave attenuation for studying fractured rocks. Previous study in one of Abu Dhabi oilfields, showed the advantage of integrating seismic attenuation to characterize fractured carbonate reservoirs [7]. Bouchaala, et al. [8] was able to separate between open and cemented fractures by using seismic wave attenuation anisotropy. Yang, et al. [9] simultaneously took into account the effect of anisotropy and attenuation in the estimation of reflection coefficient, to be implemented in seismic inversion methods to predict fracture and fluids properties.
Seismic wave attenuation describes the energy losses caused to the wave as it propagates in the subsurface. It is mainly controlled by scattering and frictional losses caused by pore fluids. Scattering is an elastic deflection of seismic wave's energy in random directions resulting in generating small amplitude multiple reflections, caused by the heterogeneities of the subsurface [10,11]. The frictional displacement at the boundary of pore fluids and solid grains caused by seismic wave motion, induces a conversion of seismic wave's energy to heat. The wave induced fluid flow (WIFF) is one of the most known frictional attenuation mechanisms [12,13]. Sensitivity of attenuation to fluids and permeability makes it as a promising seismic tool for investigating the development of reservoir systems as monitoring enhanced oil recovery (EOR) during water, steam or CO 2 injections. Indeed, comparing attenuation changes before and after injection may provide information about the flooded zones better than conventional seismic attributes such as velocity and reflectivity [14,15]. Furthermore, attenuation may be used to study fractures which are an important source of scattering and fluid-related attenuation mechanisms.
Seismic wave attenuation is often quantified with the dimensionless parameter Q, called quality factor, which its inverse, is directly proportional to attenuation. Spectral ratio (SR) [16] and centroid frequency shift (CFS) [17] techniques, are the most popular techniques for estimating Q. Zero offset vertical seismic profiling (ZVSP) waveforms are the best suitable data for SR and CFS methods, because of the good vertical sampling provided by the downgoing waves of ZVSP data. Nevertheless, methods were also successfully applied on some surface seismic data, as in [18], where the slope of spectral ratio changes over frequency at different common midpoint gathers were calculated, then the intercept was defined from the linear regression of the slope versus the squared offset as Q −1 at zero offset. Rossi, et al. [19] developed an algorithm to jointly invert velocity and Q −1 based on frequency shift method, to estimate Q −1 from 3D surface seismic data and depict the bottom simulating reflector (BSR), which separates hydrates from free gas zones. Bouchaala and Guennou [20] also succeeded to resolve the BSR based on the magnitude variation of Q −1 obtained by applying the SR technique on Ocean Bottom Cable (OBC) data acquired at the Norwegian coasts. Although SR and CFS methods were applied on some surface seismic data, these methods still have limitations in terms of accuracy and depth-resolution, they are able to provide values of Q −1 between two reflectors only.
A workflow to estimate Q −1 from 3D OBC seismic data are proposed in this study, in order to estimate attenuation with high resolution and compare it with ZVSP attenuation.
The main steps followed in the workflow were as follows: -We used the linear filter proposed by Gabor [21] to create a collection of OBC spectrum traces in time domain. The bandwidth of Gabor filter was optimized by using a 3D walkaway VSP data, acquired in the same zone as surface seismic survey; Thereby, we succeed to sample the data at various known depths and made them suitable for the application of the CFS technique as in the case ZVSP; - The spectrum traces were converted from time to depth domain by using sonic log as a velocity model; -Estimation of Q −1 by applying CFS method on traces recorded by geophones with minimum offset at each geophone line; -Perform linear regression of Q −1 versus offset to calculate the intercept, which corresponds to zero offset attenuation.
The workflow developed in this study provides a high-resolution attenuation profiles with a stability and accuracy that can be checked and obtained with the assistance of VSP data. Such performance cannot be reached by the conventional methods, as the one used in the aforementioned studies.
In this study, seismic wave attenuation was not used for EOR monitoring, due to lack of seismic data during production. Nevertheless, this study provides a workflow for estimating seismic attenuation and shows the advantage of using seismic wave attenuation in delineating the zones with high oil content and fractures density, which is crucial for EOR studies.

Study Area and Field Data
The area is one of the offshore oilfields belonging to the emirate of Abu Dhabi, United Arab Emirates ( Figure 1). Its structure is an anticline with a gentle dip, having the highest magnitude in the north and northeast, where the average dip is about 2.5 • . The stratigraphy of the oilfield mainly consists of Mesozoic-Cenozoic carbonates. The reservoir zones are stacked sections of porous and shallow water carbonates of the Lower Cretaceous Thamama Group. A two-component 3D OBC seismic survey was carried out over the oilfield from 2000 to 2001. The full-fold area of the survey was 1505 km 2 . However, in this study we focused on an area of 3 km radius around the wellhead where the 3D VSP survey was acquired ( Figure 2). A series of parallel receiver cables were laid down on the seabed with sensor groups at every 50 m. The source lines were orthogonal to the receiver cables at 100 m line spacing and 25 m shot spacing on each line. For each group, hydrophone and vertical geophone data were recorded. The air gun used as source to generate data were placed 5 m below the sea surface in most of the acquisition area, and 3 m for most shallow water areas. In the shallowest area, the source was raised to 2.5 m. The raw 3D OBC data contain an acceptable level of noise with an average magnitude of 1 µbar on the hydrophones ( Figure 3a) and double on the geophones (Figure 3b). The receiver lines were orientated parallel to the direction of the strongest currents to reduce the ambient noise levels. The strongest noise in the data is the dispersive ground-roll events that have velocities of around 900 ms −1 . Other dispersive surface waves have velocities of around 1500 ms −1 , but these waves have less overlap with the reflection energy.
In order to properly sum the hydrophone and geophone data, we converted the amplitudes of hydrophone data from pressure to velocity to be consistent with the geophone amplitudes. This was performed by performing an integration in the frequency domain: . .
. X( f ) is the converted one, and ω = 2π f is the angular frequency.
The summation of hydrophone and geophone data has highly reduced the ground-roll events and removed the ghost effect as well, which is due to the reverberation at sea surface ( Figure 3c).
The 3D VSP survey was acquired in 2015 in the southwest part of the oilfield, where 18,364 shots were generated with an average spacing of 25 m and along a spiral pattern centered at the well head. The spiral arms were separated by an average radial interval of 50 m ( Figure 2). The VSP data were acquired within an area with a radius of 3000 m from the wellhead. Seventy receivers 20 m apart were laid down in a slightly deviated well from the depth of 1145 m. Sonic and density logs were also acquired in the well ( Figure 4).   In addition, also shown are 3D OBC receiver (geophones and hydrophones) lines (green dots) that recorded the data for the closest shot (green star) to the wellhead (black triangle) with a distance of 112 m. The geophones and hydrophones that have the smallest offset and closest to a VSP shot with same offset, at each receiver line are highlighted by red dots.

Methodology
By assuming a wavelet spectrum of a Gaussian shape and a frequency-independent Q model (Futterman, 1962), Quan and Harris (1997), estimated Q −1 from two waveforms recorded at two receivers placed at depth z 1 and z 2 as follows: where σ 2 1 is the variance of the wavelet at the shallower depth z 1 , ∆t is the travel time from z 1 to z 2 and f c is the centroid frequency of a recorded signal A(z, f ) defined as: The attenuation Q −1 estimated by Equation (2), combines scattering (Q −1 sc ) and intrinsic attenuation (Q −1 in ). As explained by [22], Q −1 sc can be estimated by appying the CFS technique on synthetic waveforms, generated from the convolution of the reflectivity with a unit amplitude spike. The reflectivity can be obtained by using Gouppillaud model [23], which is based on ray tracing and Z-transform performed in horizontally layered medium. The time interval ∆t in Gouppillaud model, is defined as the one-way travel time in each layer. The relationship between transmission (T k ) and reflection (R k ) responses of a wave propagating with an angular frequency ω in a medium of kth layer with an equal one-way travel time and a reflection coefficient r k at its top, is as follows: where r k can be calculated from the acoustic impedance obtained by multiplying the sonic and density logs (Figure 4), and τ k is the two-way travel time of the kth layer.

Workflow
A workflow ( Figure 6) was developed to overcome the limitations of estimating seismic attenuation from OBC data, and to validate the calculated attenuation profiles. According to Equation (2), application of the CFS technique requires that for a given seismic shot, the spectrum must be sampled at different depths, as VSP data which are acquired at receivers placed along the vertical direction ( Figure 7a). However, OBC data are acquired at receivers laid along the seafloor, which does not enable a spectrum record at different depths. This makes the application of CFS technique for a direct estimation of attenuation from OBC waveforms, not possible. In order to overcome this limitation, a collection of spectra at different depths were generated from each single seismic trace, by applying Gabor transform which is a short time Fourier transform with a Gaussian filter [21]. The Gabor transform of a signal f (t) is defined by the following formula: where FT indicates Fourier transform, W the Gaussian filter around the time, BW is the filter width and is the angular frequency. Subsequently, a conversion from time to depth domain (Figure 7b) by using a sonic log acquired at the 3D VSP well along with a density log (Figure 4). The use of the sonic log is justified by its high depth-resolution and by the slight lateral stratigraphy changes of the subsurface.  As shows Equation (5), Gabor transform is strongly dependent on the band filter width (BW), which makes the generated spectra strongly dependent on BW. In order to correctly choose BW, we selected the closest 3D OBC seismic shot to the wellhead, then extracted the waveforms generated by this shot and recorded at the geophones and hydrophones that have the smallest offset at each of six receiver lines ( Figure 5). According to the principal of propagation reciprocity, the waveforms extracted from the OBC data recorded at the selected geophones/hydrophones and from the 3D VSP data generated by the shots very close to the selected geophones/hydrophones ( Figure 5) should be similar, and hence similar frequency profiles also. Based on that, we applied Gabor transform on the 3D OBC data by choosing a BW that produces similar OBC and VSP centroid frequency profiles. The centroid frequency profiles of the VSP and OBC seismic waveforms show the best correlation at a filter width of 5.5 Hz (Figure 8). Eventually we obtained a Q −1 profile at the six offsets after applying the CFS technique (Equation (2)) on the Q −1 profiles of the waveforms acquired by the geophones and hydrophones with minimum offset at each of the six receivers. By analyzing Q −1 crossplots versus offset at several depths, we realized that the linear regression is the best model to fit the variation of Q −1 with offset. Therefore, the intercept of the linear regression at each depth corresponds to the magnitude of Q −1 at zero offset. Therefore, attenuation profile at each OBC shot position was obtained (Figure 9).  The comparison of zero offset Q −1 with the one estimated from zero offset VSP data shows a good similarity ( Figure 10). The good similarity of the Q −1 profiles are expected since they were obtained at similar location and in the same frequency range (3-70 Hz). Such similarity also suggests that a linear relationship of Q −1 with offset is reasonable. However, differences between the Q −1 magnitude of OBC and VSP profiles can be noticed, which makes our study more qualitative than quantitative. Furthermore, there is a shift between Q −1 peaks around the depth of 2000 m, which might be due to errors in velocity. Figure 10. Zero offset Q −1 profiles obtained from 3D VSP data (black) and from a linear fitting of Q −1 estimated from 3D OBC seismic data (red), generated by the closest shot (the green star in Figure 5) to the wellhead (the black triangle in Figure 5) and recorded by hydrophones that have the smallest offsets (the red dots in Figure 5).
As mentioned above, the scattering Q −1 sc is obtained by applying the CFS method on synthetic seismograms. The 1D density and sonic logs were used to estimate the reflectivity of the medium by using Equation (4). In order to obtain Q −1 sc in the same frequency range as Q −1 , a band pass filter (3-70 Hz) were applied on the generated seismograms. Overall, the Q −1 sc profile shows high magnitudes (Figure 11a) which explains the high magnitudes of Q −1 . The heterogeneous nature of carbonate subsurface of Abu Dhabi, explains the high magnitudes of Q −1 sc . The heterogeneity and the complex nature of the oilfield is notable through the high variation of velocity with depth (Figure 11b). Schoenberger and Levin [24] suggested that Q −1 = Q −1 in + Q −1 sc , meaning that by subtracting scattering from the total attenuation, the intrinsic can be obtained separately. However, it is difficult to prove this assumption, especially in complex media such as carbonates and in the case of non-zero offset data. Therefore, for sake of error we did not apply the separation between Q −1 sc and Q −1 in as suggested by [24]. Therefore, any interpretation related intrinsic attenuation in this study is based on the comparison between Q −1 and Q −1 sc .

Estimation of Attenuation Profiles
By performing Gabor transform with a BW of 5 Hz and estimating the zero offset Q −1 from linear regression as shown in the workflow, 3D attenuation volumes were obtained from hydrophone, geophone data, and their summation. The results are displayed by a succession of 2D attenuation profiles obtained at selected source lines (Figures 12 and 13).
We notice a good qualitative similarity between the Q −1 profiles obtained from hydrophone (Figure 12b-d), geophone (Figure 13a,c,e) and summation (Figure 13b,d,f) data. However, the magnitude of the attenuation estimated from the hydrophone data are higher than those obtained from the other two datasets. This might be due to the noise effect, which is more prominent in the geophone and summation data than in the hydrophone data. Furthermore, due to shallow sea bottom, the wavelength is larger than the water depth (6 m-20 m), which makes reverberation removal difficult. Slight lateral variation can be observed on the 2D attenuation profiles, while the vertical variation is much significant.
The high attenuation anomalies noticed in Simsima, Fiqa and Nahr Umr Formations are due to high scattering in these formations (Figure 11a). This is in agreement with [25], where the high attenuating nature of these formations in some oilfields of Abu Dhabi was reported. The reservoir zones display high attenuation and moderate scattering, which lets us suggest that the attenuation is dominated by the intrinsic.

Discussion
By performing Gabor transform and a linear regression of Q −1 with offsets, we obtained an attenuation profile at each source location used to generate the OBC data. Overall, the attenuation accuracy is reasonable, where the R-squared of the linear regression is more than 0.5 for all sources and more than 0.8 for many sources.
The attenuation magnitude shows a slight lateral variation along source lines, compared to the vertical direction. This is in agreement with seismic sections (Figure 14a) which show slight lateral variation in amplitudes. A recent velocity model obtained from a full waveform inversion (FWI) [26] generated from VSP data acquired along an azimuth of 56 • , also revealed a much smaller variation of velocity in the lateral direction than in the vertical one (Figure 14b). It can be noticed that the attenuation model obtained at the same azimuthal direction (56 • ), displays more lateral variation than the velocity model (Figure 14c). This demonstrates the higher sensitivity of attenuation than velocity to the variation of petrophysical and lithological properties. Such sensitivity might be useful for enhance-ment oil recovery studies, since it can help to delineate the zones with great hydrocarbon potential, which are usually noticed by their high intrinsic attenuation anomalies.
Seismic wave attenuation was suggested as a direct detector of oil and gas, in several previous studies carried out in sandstone [27,28]. However, for subsurfaces dominated by carbonate lithology, as in the case of this study, it is still too early to consider attenuation as a direct indicator of hydrocarbon. This is because of the complex texture of carbonate media, resulting in difficulty in the understanding of attenuation mechanisms and in their estimation with reasonable acuuracy. Nevertheless, in the case of this study, the reservoir zones are clearly identifiable by their high attenuation anomaly, whereas scattering has moderate magnitudes in these zones ( Figure 11). Hence, we suggest that in the reservoir zones, energy losses are dominated by intrinsic attenuation, which is the result of fluidrelated mechanisms, such as WIFF mechanism known to be significant in the reservoir zones. There are two types of WIFF mechanism: global and squirt flow. The former is important in highly porous and permeable rocks [29,30], while the latter is more present in fractured media [31,32]. The high fracture density [8] and porosity in our reservoirs suggest the co-existence of global and squirt flow mechanisms. Bouchaala, et al. [33] who studied seismic attenuation in a resevoir of one of Abu Dhabi oilfields, made the evidence of the existence of global flow mechanism in porous and free-fractures zones, and of squirt flow in dense and fractured zones.
Some formations above the reservoir zones also have high attenuation anomalies, such as Nahr Umr, Salabikh, and Simsima Formations (Figure 12b,d). However, opposite to reservoir zones where fluid-related mechanisms are important, the scattering is the dominant attenuation component in these formations (Figure 11a). The high magnitudes of scattering reflect the heterogeneous and complex lithology of these formations. Indeed, fine layering is the main cause of high scattering at the base of middle cretaceous Nahr Umr Formation which is made of variegated shales and occasional thin limestones [34]. Many full-bore formation micro-imager observations in Abu Dhabi oilfields show a pronounced thin bedding planes in Nahr Umr Formation, which is usual in shaly rocks [34,35]. The upper middle cretaceous Salabikh Formation is a sequence of three members made of highly variable succession of beds and consists of lime mudstones, soft dense argillaceous limestones with very calcareous shale interbeds [34]. Such bedding explains the high scattering in this formation. The Upper Cretaceous Simsima Formation is highly heterogeneous, and its lithology is dominated by depositional packstones, in a shallow water environment, representing a regressive sea level phase. Furthermore, Simsima Formation is an important aquifer in Abu Dhabi as well as an oil-bearing reservoir in the Shah Field [36,37]. Therefore, Simsima Formation is a favorable medium for scattering and WIFF mechanisms. A previous study performed in another oilfield located in Abu Dhabi has reported high attenuation anomalies in Simsima Formation.
The NW part of the study area displays negative attenuation values at depths between 1700 m and 1850 m and below 2400 m (Figure 12d). Mateeva [38] suggested that the most probable reason behind negative attenuation is the sensitivity of the CFS technique to the interference phenomenon caused by short wavelength upgoing waves, such as small amplitude multiples. Another reason might be the disturbance caused by noise. Decreasing the vertical resolution of Q −1 might be a solution to obtain more stable results, and to avoid the negative magnitudes of Q −1 . However, the resolution of Q −1 is important for its comparison with the petrophysical properties of the subsurface, which are typically acquired with a vertical resolution of 0.1524 m along the wellbore. Therefore, a compromise should be performed between the resolution and stable magnitudes of Q −1 . On this basis, the resolution of our attenuation profiles was set equals to 100 m. Furthermore, negative magnitudes of scattering were noticed, such as those around the depth of 1745 m where Q −1 sc = −0.04 (Figure 11a). Liner [39] explained negative scattering by negative normal reflection coefficient (R 0 ), which is the result of negative impedance contrast. Based on the continuity relationship R 0 + T 0 = 1, the normal transmission coefficient (T 0 ) has to be larger than one, meaning an amplification in the amplitude of transmitted waves which results in negative Q −1 sc . In this study, qualitative results were presented and the 3D variation of attenuation in the subsurface of an oilfield was shown. The results pointed out the advantages of using seismic wave attenuation for detecting and delineating, heterogeneous and fluid-saturated zones, such as reservoirs and aquifers. However, estimating attenuation from large offset 3D seismic data is still challenging, because at high offsets the interference phenomenon is important, which may reduce the accuracy of attenuation estimates.
Bouchaala, et al. [40] showed through a synthetic data with similar geometry to nonzero offset VSP data, that the error induced by the offset is noticeable at offsets higher than 2000 m. Nevertheless, for larger offsets they showed that the variation of the attenuation profile remains acceptable. The offsets used in this study vary from 20 m to 4000 m.
Furthermore, a correct separation between scattering and intrinsic attenuation can eliminate negative values of Q −1 and enhance interpretation. However, the methodology for estimating and separating scattering from 3D seismic data is still premature and challenging.

Conclusions
The workflow developed in this study allowed us to apply the CFS technique on OBC seismic data, in order to estimate seismic wave attenuation. The limitation of the CFS technique due to the location of OBC receivers on the seafloor, was overcome by applying an optimized Gabor transform to create a collection of spectrum traces in the time domain, followed by a conversion from time to depth domain. The results showed high attenuation anomalies in highly saturated and heterogeneous zones, such as the reservoir zones and Simsima Formation. Such results show the advantage of using seismic wave attenuation for petroleum exploration and reservoir characterization. A slight lateral variation at each seismic line was revealed by the 2D attenuation profiles, which agrees with velocity profiles and results of full waveform inversion.