TY - JOUR
AU1 - Matsushima,, Jun
AU2 - Ali, Mohammed, Y
AU3 - Bouchaala,, Fateh
AB - Summary Intrinsic attenuation has the potential to clarify the physical properties of earth materials, whereas positive utilization of seismic scattering has been recognized as useful in investigating the heterogeneities thereof. However, it has not been easy to separate between intrinsic attenuation and scattering effects. The zero-offset vertical seismic profiling (VSP) survey is recognized to be ideal for attenuation estimation because the VSP survey observes direct waveforms propagating through subsurface formations. The most popular method of separating intrinsic and scattering attenuation generates numerical zero-offset VSP data from known velocity data, such as sonic velocity logs, and isolates the intrinsic attenuation by subtracting the synthetic scattering attenuation from the total attenuation. Numerical experiments have demonstrated the limitation of this conventional method due to its assumption that the intrinsic and scattering attenuation are independent each other. We should take into account the mutual interactions between the intrinsic and scattering attenuation. In order to overcome this limitation, we herein propose a novel method for separating intrinsic and scattering attenuation for zero-offset VSP data by reforming the modified median frequency shift (MMFS) method with seismic interferometry (SI) under the assumption that intrinsic and scattering attenuation are frequency independent and frequency dependent, respectively. The proposed method can simultaneously derive both intrinsic and scattering attenuation estimates from only VSP data without sonic and density well-log data. The numerical experiments of the present study also investigate the importance of parameter optimization in applying pre-processing filters in order to balance the resolving power and noise reduction effect. Finally, we apply the proposed method to zero-offset VSP data acquired in an onshore Abu Dhabi oil field in order to demonstrate the applicability of the proposed method in investigating heterogeneities of subsurface formations and discuss its limitations. Downhole methods, Interferometry, Seismic attenuation 1 INTRODUCTION The seismic attenuation estimated from seismic data includes two attenuation components: intrinsic attenuation and scattering attenuation. The total attenuation is widely assumed to be a linear summation of the intrinsic attenuation and the scattering attenuation (e.g. Spencer et al.1982). Note that this assumption holds for the prerequisite of extraordinary scattering, in which the seismic wave propagation should be described by the diffusion process (Warren 1972). Wegler (2004) discussed the applicability of the diffusion equation for the strong scattering case. Intrinsic attenuation, which transforms seismic energy into heat, has the potential to clarify the physical properties of earth materials because composite materials, such as pore spaces or fractures that are fully or partially saturated with fluids, exhibit viscoelastic properties that depend on the effective stress or fluid type (Winkler & Nur 1982). Although other possible intrinsic attenuation mechanisms, such as the effects of wetting on grain boundaries (Johnston et al.1979) or viscous shear relaxation (Walsh 1969), have been proposed, many researchers believe that the fluid flow induced by the seismic wave propagation mechanism is the most influential mechanism that can account for observed seismic attenuation over a broadband frequency range (e.g. Barton 2006). In contrast, scattering attenuation caused by acoustic impedance heterogeneities is a function of the wavelength of the seismic wave and the characteristic size of the heterogeneities (Wu & Aki 1988). Scattered seismic events have the potential to provide information about the orientation, spacing, and density of a fracture (Vlastos et al.2003; Wills et al.2006; Vlastos et al.2007). By conducting numerical wave propagation in media with discrete distributions of fractures, Vlastos et al. (2003) pointed out the importance of the spatial distribution of fractures and the fracture size relative to the wavelength. Willis et al. (2006) proposed the scattering index method to investigate the characteristics of fracture systems in subsurface formations using scattered coda waves. Through a comprehensive study of multiple scattering in a stochastic–deterministic model of an evolving fractured network, Vlastos et al. (2007) concluded that main factors affecting scattering attenuation are fracture density, fracture scale, and fracture orientation. They further implied that there is a potential to use scattering attenuation to characterize such fracture properties, and that scattering attenuation strongly depends on frequency and its dependency has a complex variation at different stages of fracture evolution. Fang et al. (2014) reported that if the wavelength of the seismic wave is comparable to the scale of the fractures, then scattered waves are significantly generated, allowing us to characterize reservoir fracture systems. Fractures can also cause significant scattering effects, which are most significant for open fractures. Although fracture orientation varies with historical stress regimes, open fractures are generated in the direction parallel to the current maximum horizontal stress direction (Burns et al.2007). Scattering effects generated by fractures have strong directional characteristics in accordance with fracture orientation. Burns et al. (2007) indicated that the orientation and density of the fracture systems can be evaluated using scattered waves. Vlastos et al. (2006) numerically investigated the effect of pore pressure changes on seismic wave propagation in a fracture network. They found the systematic effects of pore pressure changes on seismic waveforms and attenuation, and identified potential ways to estimate pore pressure changes from seismic data. Geological and geophysical studies of fractures systems indicate the fracture size on a scale ranging from microns (microcracks) to meters. Ali & Worthington (2011) reported the existence of complicated pore networks ranging in scale from micropores to megakarst in carbonate reservoirs through multi-stage diagenetic histories. However, it is unusual to clearly observe scattered waves caused by fractures because of the potentially weak signals of the scattered energy (Willis et al.2006). Until now, it has been widely believed that complicated seismic waves caused by small-scale heterogeneities degrade seismic reflection data, which has the potential to cause difficulties when delineating subsurface structures or characterizing physical properties by seismic methods. This difficulty is due to small-scale heterogeneities such as fractures. However, positive utilization of seismic scattering may enable us to characterize such heterogeneities. Separating total attenuation into intrinsic attenuation and scattering attenuation has not been easy. Lerche & Menke (1986) proposed a separation method under the assumption that the shear mode causes primarily intrinsic attenuation and the assumption that a scaling law exists between the random fluctuations of the bulk and shear modulus and density. A number of studies (e.g. Schoenberger & Levin 1974, 1978; Ganley & Kanasewich 1980; Maresh et al.2006; Matsushima et al.2006; Shaw et al.2008; Bouchaala et al.2016a, 2016b) calculated the scattering attenuation caused by multiples from sonic velocity and density well-log data using 1-D synthetic seismic modelling and then subtracted the scattering attenuation from the total attenuation to isolate intrinsic attenuation. O'Doherty & Anstey (1971) formulated the effects of short-path multiple reflections on the normal-incidence plane-wave reflection and transmission responses in a finely and vertically layered random medium. Shapiro et al. (1994) generalized the O'Doherty-Anstey formula to analytically describe the transmission response of the angle-dependent plane wave in a randomly multilayered media using a correlation function of sonic velocity and density well-log data. Although, in general, the scattering attenuation is estimated in a 1-D layered stratigraphy model and the effects of 3-D heterogeneities are underestimated (Gist 1994). Huang et al. (2009) estimated the 3-D scattering effect by combining petrophysical stochastic models and 3-D finite difference modelling. Although we should estimate the 3-D scattering effects induced by 3-D heterogeneities, Matsushima et al. (2011) pointed out the difficulty of constructing a 3-D heterogeneity model. Nevertheless, in the above studies, scattering attenuation estimation is based on the assumption that the intrinsic and scattering attenuation are independent each other, and furthermore attention has not been given to positive utilization of seismic scattering in order to characterize the heterogeneities of subsurface formations. In this paper, we first demonstrate the limitation of the conventional approach due to its assumption that the intrinsic and scattering attenuation are independent each other. We should take into account the mutual interactions between the intrinsic and scattering attenuation. In order to overcome this limitation, we propose a novel method for separating intrinsic and scattering attenuation for zero-offset vertical seismic profiling (VSP) data. We then investigate the effectiveness of the proposed method by applying the proposed method to synthetic data. Finally, we adopt the proposed novel method to real zero-offset VSP data acquired at an onshore Abu Dhabi oil field in order to examine the applicability and limitations of the proposed method in investigating heterogeneities of subsurface formations. 2 METHODOLOGY Matsushima et al. (2016) developed a stable method to estimate seismic attenuation from zero-offset VSP data by combining the modified median frequency shift (MMFS) method and seismic interferometry (SI) (hereinafter referred to as the SI+MMFS method). By modifying the SI+MMFS method, we can separate the intrinsic and scattering attenuation for zero-offset VSP data. In the following, the proposed method is described in part following Matsushima et al. (2016). In the SI+MMFS method, as a first step, the source position of zero-offset VSP data is redatumed to the receiver position by adopting SI. Then, the MMFS method can be applied to the redatumed VSP data. Fig. 1(a) shows a simplified zero-offset VSP measurement with a source located at the surface and three receiver locations R1, R2 and R3 at their depths of z1, z2 and z3, respectively. Based on these three traces received at three different depths, we can create new seismic traces TR1 and TR2 by the application of SI (Fig. 1b). Here, we use deconvolution interferometry (DCI; Vasconcelos & Snieder 2008) in the application of SI. The plane-wave solution including the exponential attenuation law can be expressed in the frequency domain as: \begin{equation}U\left( {z,\omega } \right) = S(\omega ){e^{ - \alpha z}}{e^{ikz}},\end{equation}(1) where z is the receiver depth from the surface, ω is the angular frequency, S(ω) is the source function at the surface in the angular frequency domain, |$\alpha {\rm{ }}( { {=} {{\pi f} \mathord{/ {\vphantom {{\pi f} {Qv}}} } {Qv}}} )$| is the attenuation coefficient (f: frequency, Q: quality factor, and v: velocity), k is the wavenumber, and i is the imaginary unit. Thus, the data recorded at three different receivers (R1, R2 and R3) can be expressed as U(z1, ω), U(z2, ω) and U(z3, ω), respectively. Figure 1. Open in new tabDownload slide (a) A simplified zero-offset VSP measurement with a source located at the surface and three receiver locations R1, R2 and R3 in a vertical borehole. In the figure, z1, z2 and z3 are the travel distances for the three receiver locations. (b) Two seismic traces created by measurement configuration using SI. Figure 1. Open in new tabDownload slide (a) A simplified zero-offset VSP measurement with a source located at the surface and three receiver locations R1, R2 and R3 in a vertical borehole. In the figure, z1, z2 and z3 are the travel distances for the three receiver locations. (b) Two seismic traces created by measurement configuration using SI. The application of DCI to U(z1, ω) and U(z2, ω) yields the first source–receiver pair trace: \begin{equation}{\rm{T}}{{\rm{R}}_1}\left\langle {{\rm{DCI}}} \right\rangle = \frac{{U({z_2},\,\omega )}}{{U({z_1},\,\omega )}} = {e^{ - \alpha (z2 - z1)}}{e^{ik(z2 - z1)}}.\end{equation}(2) Similarly, the application of DCI to U(z1, ω) and U(z3, ω) yields the second source-receiver pair trace: \begin{equation}{\rm{T}}{{\rm{R}}_2}\left\langle {{\rm{DCI}}} \right\rangle = \frac{{U({z_3},\,\omega )}}{{U({z_1},\,\omega )}} = {e^{ - \alpha (z3 - z1)}}{e^{ik(z3 - z1)}}.\end{equation}(3) In the present study, we do not apply a regularization in the calculation of deconvolution (Mehta et al.2007) to stabilize attenuation estimation on lower signal-to-noise ratio (SNR) data. This is because the quality of field zero-offset VSP data used in the present study is sufficiently high. Here, we assume the configuration of the multireceiver sonic logging measurement tool. The amplitude of the wave from a source at a receiver at separation distance d and angular frequency ω is expressed as: \begin{equation}A\left( {d,\omega } \right) = B{e^{ - \pi \omega \Delta t/2Q}}\end{equation}(4) where Δt is the travel time for the separation distance d and B is a parameter that is assumed to be frequency independent and includes various factors, such as the geometric spreading effect, the source and receiver function, and the coupling effect at the source and receiver positions. Eq. (4) can be rewritten after taking the logarithm and rearranging as: \begin{equation}\frac{{2\,{\rm ln}\,A\left( {d,\omega } \right)}}{\omega } = {-}{Q^{-1}}\left( z \right)\Delta t + \frac{{2\,{\rm ln}\,B}}{\omega }.\end{equation}(5) By taking the dependency on frequency or depth into account, we can rewrite eq. (5) as follows: \begin{equation}\Phi (z,\omega ) = - \mathop Q\nolimits^{ - 1} \left( z \right)\Delta t + L(\omega )\end{equation}(6) where |$\Phi ( {z,\omega } ) = \frac{{2\,{\rm ln}\,A( {d,\omega } )}}{\omega },$| and |$L( \omega ) = \frac{{2\,{\rm ln}\,B}}{\omega }.$| We apply the concept of median frequency shift (Frazer et al.1997) to Φ(z) as follows: \begin{eqnarray} \hat{\Phi }(z) = \mathop \omega \limits^{{\rm{median}}} \left\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \left\{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \left\{ {\Phi (z,\omega )} \right\}} \right\}} \right\}\!.\nonumber\\ \end{eqnarray} (7) Eq. (7) indicates that |$\hat{\Phi }(z)$| is independent of frequency. Therefore, the final absolute attenuation profile |$Q_{{\rm{average}}}^{ - 1}(z)$| calculated from |$\hat{\Phi }(z)$| is independent of frequency. Here, we define the fluctuation portion from the frequency-dependent component as follows: \begin{eqnarray} {\hat{\Phi }_F}(z,\omega ) \!&=&\! \left\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \left\{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \left\{ {\Phi (z,\omega )} \right\}} \right\}} \right\} - \hat{\Phi }(z).\nonumber\\ \end{eqnarray} (8) Eq. (8) indicates that |${\hat{\Phi }_F}(z,\omega )$| is dependent on frequency. Here, by assuming that intrinsic and scattering attenuation are frequency-independent and frequency-dependent, respectively, we can estimate the frequency-independent intrinsic attenuation from |$\hat{\Phi }(z)$| and the frequency-dependent scattering attenuation from |${\hat{\Phi }_F}(z,\omega )$|⁠. The appropriateness of this assumption will be investigated and discussed later herein. The resulting |$\hat{\Phi }(z)$| and |${\hat{\Phi }_F}(z,\omega )$| are estimates of |$ {-} \mathop Q\nolimits^{ - 1} (z)\Delta t(z)$| plus constant term |$\hat{L}$| and |$ - Q_F^{ - 1}(z,\omega )\Delta t(z)$| plus frequency-dependent term |$\hat{L}(\omega )$|⁠, respectively: \begin{equation} \hat{\Phi }(z) = - \mathop Q\nolimits^{ - 1} (z)\Delta t(z) + \hat{L} \end{equation}(9) \begin{equation}{\hat{\Phi }_F}(z,\omega ) = - Q_F^{ - 1}(z,\omega )\Delta t(z) + \hat{L}(\omega ).\end{equation}(10) Assuming |$\hat{L}$| and |$\hat{L}(\omega )$| to be constant over the depth of interest, we finally obtain |$\mathop Q\nolimits^{ - 1} (z)$| and |$Q_F^{ - 1}(z,\omega )$| from |$\hat{\Phi }(z)$| and |${\hat{\Phi }_F}(z,\omega )$| respectively, as follows: \begin{equation}{Q^{ - 1}}(z) = \frac{{\left( { - \hat{\Phi }(z) + \hat{\Phi }(\zeta ) + {Q^{ - 1}}(\zeta )\Delta t(\zeta )} \right)}}{{\Delta t(z)}}\end{equation}(11) \begin{equation}Q_F^{ - 1}(z,\omega ) = \frac{{\left( { - {{\hat{\Phi }}_F}(z,\omega ) + {{\hat{\Phi }}_F}(\zeta ,\omega ) + Q_F^{ - 1}(\zeta ,\omega )\Delta t(\zeta )} \right)}}{{\Delta t(z)}}\end{equation}(12) where Q− 1(ζ) and |$Q_F^{ - 1}(\zeta ,\omega )$| are arbitrarily values at an arbitrarily depth of ζ. For the i-th source-receiver pair at a certain depth, |${\hat{\Phi }_i}(z)$| and |${\hat{\Phi }_{F\_i}}(z,\omega )$|⁠, respectively, can be expressed as follows: \begin{equation} {\hat{\Phi }_i}(z) = - \mathop Q\nolimits^{ - 1} (z)\Delta {t_i}(z) + {\hat{L}_i} \end{equation}(13) \begin{equation}{\hat{\Phi }_{F\_i}}(z,\omega ) = - Q_F^{ - 1}(z,\omega )\Delta {t_i}(z) + {\hat{L}_i}(\omega ).\end{equation}(14) The amplitude spectra of TR1〈DCI〉 and TR2〈DCI〉 are substituted into A(d, ω) on the left-hand side in eq. (5), that is, |${\hat{\Phi }_1}(z)$| and |${\hat{\Phi }_2}(z)$|⁠, and |${\hat{\Phi }_{F\_1}}(z,\omega )$| and |${\hat{\Phi }_{F\_2}}(z,\omega )$| are calculated from TR1〈DCI〉 and TR2〈DCI〉 data, respectively. Furthermore, Δt1(z) and Δt2(z) are obtained by the first-arrival times of VSP data. The application of a linear fitting on eqs (13) and (14) enables us to obtain the attenuation values of Q− 1(z) and |$Q_F^{ - 1}(z,\omega )$|⁠, respectively, over the depth of interest. The calculated attenuation values of Q− 1(z) and |$Q_F^{ - 1}(z,\omega )$| are regarded as reference values of |$\mathop Q\nolimits^{ - 1} (\zeta )$| and |$Q_F^{ - 1}(\zeta ,\omega )$|⁠, respectively. Then, we obtain absolute attenuation values of Q− 1(z) and |$Q_F^{ - 1}(z,\omega )$| using eqs (11) and (12), respectively. In order to stabilize scattering attenuation estimation in eq. (12), we use the reference value of |$\mathop Q\nolimits^{ - 1} (\zeta )$| instead of |$Q_F^{ - 1}(\zeta ,\omega )$|⁠. In order to reduce the arbitrariness in the selection of |$\mathop Q\nolimits^{ - 1} (\zeta )$|⁠, an average attenuation value is calculated over reference depths of ζ ranging from z1 to zN using the following equations: \begin{equation} Q_{{\rm{average}}}^{ - 1}(z) = \sum\limits_{\zeta = {z_1}}^{{z_N}} {\frac{{Q_{{\rm{absolute}}}^{ - 1}(z,\zeta )}}{N}} \end{equation}(15) \begin{equation}Q_{F\_{\rm{average}}}^{ - 1}(z,\omega ) = \sum\limits_{\zeta = {z_1}}^{{z_N}} {\frac{{Q_{F\_{\rm{absolute}}}^{ - 1}(z,\zeta ,\omega )}}{N}} \end{equation}(16) where |$Q_{{\rm{absolute}}}^{ - 1}(z,\zeta )$| and |$Q_{F\_{\rm{absolute}}}^{ - 1}(z,\zeta ,\omega )$| are absolute attenuation estimates at an arbitrarily depth of ζ. The resulting values of |$Q_{{\rm{average}}}^{ - 1}(z)$| and |$Q_{F\_{\rm{average}}}^{ - 1}(z,\omega )$| are the final intrinsic and scattering attenuation estimates, respectively. 3 NUMERICAL TESTS Numerical tests were conducted in order to evaluate the performance of the proposed method on synthetic VSP data by comparing the proposed method with the conventional approach based on the spectral ration (SR) method proposed by McDonal et al. (1958). Three different 1-D numerical velocity and attenuation models with a constant density are shown in Figs 2(a), (d) and (g). In Figs 2(a) and (d), the heterogeneous zone ranging from 1950 to 2150 m consists of alternating base-velocity (2000 m s−1) and anomalous-velocity (2200 m s−1) layers with a thickness of 10 m. The Q values of three different numerical models shown in Figs 2(a), (d) and (g) are 200, 1010, and 200, respectively. Note that our intrinsic attenuation is based on an exponential decay with a frequency-independent Q model (Aki & Richards 1980). Therefore, the calculated attenuation values from three different models in Figs 2(a), (d) and (g) are considered as the total, scattering, and intrinsic attenuation, respectively. Figure 2. Open in new tabDownload slide (a, d, g) Simple 1-D attenuation and velocity model with a constant density. The alternating velocity zone in (a, d) ranging from 1950 to 2150 m consists of the alternation of base velocity (2000 m s−1) and anomalous velocity (2200 m s−1) layers with a thickness of 10 m. The Q values for three different numerical models shown in (a), (d), and (g) are 200, 1010, and 200, respectively. (b, e, h) Synthetic zero-offset VSP waveforms produced from the viscoelastic numerical models shown in panels (a), (d) and (g), respectively. (c, f, i) Application of first-arrival alignment to the waveforms shown in panels (b), (e), and (h), respectively. Figure 2. Open in new tabDownload slide (a, d, g) Simple 1-D attenuation and velocity model with a constant density. The alternating velocity zone in (a, d) ranging from 1950 to 2150 m consists of the alternation of base velocity (2000 m s−1) and anomalous velocity (2200 m s−1) layers with a thickness of 10 m. The Q values for three different numerical models shown in (a), (d), and (g) are 200, 1010, and 200, respectively. (b, e, h) Synthetic zero-offset VSP waveforms produced from the viscoelastic numerical models shown in panels (a), (d) and (g), respectively. (c, f, i) Application of first-arrival alignment to the waveforms shown in panels (b), (e), and (h), respectively. Numerical zero-offset VSP data were produced in the frequency domain using the 1-D viscoelastic waveform modelling method (Yang et al.2009). The Ricker wavelet with a dominant frequency of 50 Hz was used as a source signal at a depth of 0 m and a total of 101 positions were used as receiver stations with depths ranging from 1500 to 2500 m at intervals of 10 m. Figs 2(b), (e) and (h) show the synthetic shot gathers for three different numerical models shown in Figs 2(a), (d) and (g), respectively. In each shot gather, first-arrival waveforms were extracted with a Hanning window (length: 120 ms) centred at the first peak and were aligned at first arrivals (Figs 2c, f and i). In the present study, we applied a median filter and a trace mixing to reduce upgoing reflections. Here, we investigated the effect of the application of these two filters on attenuation estimates. Figs 3(a), (c) and (e) show the Fourier amplitude spectra for the waveforms shown in Figs 2(c), (f) and (i), respectively. Similarly, Figs 3(b), (d) and (f) show the Fourier amplitude spectra after the application of a seven-trace median filter and a seven-trace mixing filter to the waveforms shown in Figs 2(c), (f) and (i), respectively. Figure 3. Open in new tabDownload slide (a, c, e) Fourier relative spectral amplitude for the waveforms shown in Figs 2(c), (f), and (i), respectively. (b, d, f) Fourier relative spectral amplitude after the application of a seven-trace median filter and seven-trace mixing filter to the waveforms shown in Figs 2(c), (f), and (i), respectively. Figure 3. Open in new tabDownload slide (a, c, e) Fourier relative spectral amplitude for the waveforms shown in Figs 2(c), (f), and (i), respectively. (b, d, f) Fourier relative spectral amplitude after the application of a seven-trace median filter and seven-trace mixing filter to the waveforms shown in Figs 2(c), (f), and (i), respectively. First, we applied the SR method to three different synthetic models with two different pre-processing conditions. The frequency range for attenuation estimation was 10–100 Hz. Figs 4(a)–(c) show the attenuation results for the amplitude spectra shown in Figs 3(a), (c), and (e), respectively. Then, we obtained the isolated intrinsic attenuation (Fig. 4d) by subtracting the scattering attenuation (Fig. 4b) from the total attenuation (Fig. 4a). The attenuation results shown in Fig. 4(c), which were estimated from the data in Fig. 3(e), correlate well with the given attenuation value. Comparison of the attenuation results in Fig. 4(c) with the isolated intrinsic attenuation results shown in Fig. 4(d) reveals the significant difference between these results, meaning that the estimation of scattering attenuation is not appropriate for separating intrinsic and scattering attenuation. Similarly, Figs 5(a)–(c) show the attenuation results for the amplitude spectra shown in Figs 2(c), (f) and (i), respectively. Fig. 5(d) shows the isolated intrinsic attenuation obtained by subtracting the scattering attenuation (Fig. 5b) from the total attenuation (Fig. 5a). Comparison of the true intrinsic attenuation (Fig. 5c) and the isolated intrinsic attenuation (Fig. 5d) reveals these results to be significantly different. However, the degree of difference between Figs 5(c) and (d) is smaller than that between Figs 4(c) and (d). This is due to the smoothing effect of a trace mixing filter. Figure 4. Open in new tabDownload slide (a, b, c) Attenuation results obtained using the SR method for the amplitude spectra shown in panels (a), (c) and (e), respectively. (d) Isolated intrinsic attenuation by subtracting the scattering attenuation (b) from the total attenuation (a). Figure 4. Open in new tabDownload slide (a, b, c) Attenuation results obtained using the SR method for the amplitude spectra shown in panels (a), (c) and (e), respectively. (d) Isolated intrinsic attenuation by subtracting the scattering attenuation (b) from the total attenuation (a). Figure 5. Open in new tabDownload slide (a, b, c) Attenuation results obtained using the SR method for the amplitude spectra after the application of a seven-trace median filter and a seven-trace mixing filter, shown in Figs 3(b), (d) and (f), respectively. (d) Isolated intrinsic attenuation by subtracting the scattering attenuation (b) from the total attenuation (a). Figure 5. Open in new tabDownload slide (a, b, c) Attenuation results obtained using the SR method for the amplitude spectra after the application of a seven-trace median filter and a seven-trace mixing filter, shown in Figs 3(b), (d) and (f), respectively. (d) Isolated intrinsic attenuation by subtracting the scattering attenuation (b) from the total attenuation (a). Second, we applied the proposed method to the numerical data described above. Note that the proposed method can directly derive both intrinsic and scattering attenuation estimates from only the total attenuation model shown in Fig. 3(a) or (b). Figs 6(a) and (c) show the intrinsic attenuation results for two different numerical models in Figs 3(a) and (b), respectively. Similarly, Figs 6(b) and (d) show the scattering attenuation results for these two different models, respectively. Figs 7(a)–(e) show Φ(z, ω) in eq. (6), |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| in eq. (7), |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| sorted in increasing order of amplitude in the ω direction at each z, |$\hat{\Phi }(z)$| in eq. (7), and |${\hat{\Phi }_F}(z,\omega )$| in eq. (8), respectively. Fig. 6(a) shows the large intrinsic attenuation values at the both the top and bottom of the alternation zone, and, within the alternation zone, we can also see the slightly fluctuating values and their averages correspond approximately to the given attenuation value (i.e. 1/Q = 0.005). On the other hand, in Fig. 6(b), the alternating zone around a frequency of 40 to 60 Hz and its cyclic pattern correlate well with that of the given model, as shown in Fig. 2(a). By comparing Figs 6(a) and (c), the application of both a median filter and a trace mixing filter smoothed the intrinsic attenuation estimates, making these estimates approach the given attenuation values. By comparing Figs 6(b) and (d), the application of both a median filter and a trace mixing filter smoothed the scattering attenuation estimates, diminishing the cyclic pattern of the alternation zone. We discuss these results in detail later herein. Figure 6. Open in new tabDownload slide (a, c) Intrinsic (frequency-independent component) attenuation results obtained by the proposed method for two different pre-processing conditions: no pre-processing filter and pre-processing using a median filter and a trace mixing filter. (b, d) Scattering (frequency-dependent component) attenuation results obtained using the proposed method for two different pre-processing conditions: no pre-processing filter and pre-processing using a median filter and a trace mixing filter. Figure 6. Open in new tabDownload slide (a, c) Intrinsic (frequency-independent component) attenuation results obtained by the proposed method for two different pre-processing conditions: no pre-processing filter and pre-processing using a median filter and a trace mixing filter. (b, d) Scattering (frequency-dependent component) attenuation results obtained using the proposed method for two different pre-processing conditions: no pre-processing filter and pre-processing using a median filter and a trace mixing filter. Figure 7. Open in new tabDownload slide Each step of the application of the proposed method to numerical data: (a) Φ(z, ω) in eq. (6), (b) |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| in eq. (7), (c) |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| sorted in increasing order of amplitude in the direction of ω at each depth, (d) |$\hat{\Phi }(z)$| in eq. (7), (e) |${\hat{\Phi }_F}(z,\omega )$| in eq. (8). Figure 7. Open in new tabDownload slide Each step of the application of the proposed method to numerical data: (a) Φ(z, ω) in eq. (6), (b) |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| in eq. (7), (c) |$\{ {\Phi (z,\omega )-\mathop z\limits^{{\rm{median}}} \{ {\Phi (z,\omega )-\mathop \omega \limits^{{\rm{mean}}} \{ {\Phi (z,\omega )} \}} \}} \}$| sorted in increasing order of amplitude in the direction of ω at each depth, (d) |$\hat{\Phi }(z)$| in eq. (7), (e) |${\hat{\Phi }_F}(z,\omega )$| in eq. (8). 4 APPLICABILITY TO FIELD DATA We applied the proposed method to field zero-offset VSP data in order to investigate the applicability of the proposed method. The field zero-offset VSP data were acquired at an onshore Abu Dhabi oil field. A vibroseis source was used to generate sweep signals in the frequency range of 6–130 Hz with a time sampling rate of 2 ms. A total of 461 levels were acquired by a clamped three-component geophone over a depth range of from 121 to 3258 m with a receiver spacing of approximately 7.5 m. Here, we used the vertical component of acquired VSP data in the depth range of from 1200 to 3200 m (Fig. 8a). We applied a nine-trace median filter to isolate downgoing waves, and first-arrival waveforms were then extracted with a Hanning window (length: 100 ms) centred at the first peak and were aligned at first arrivals (Fig. 8b). Trace mixing filters with four different traces (one, three, five, and seven traces) were then applied. Figs 9(a)–(c) show the relative Fourier amplitude spectra, the centroid frequency (calculated from 20 to 130 Hz) for the isolated first-arrival waveforms, and the interval velocity profile derived from first-arrival time selection on VSP data, respectively. We then adopted the proposed method at a receiver separation of 7.5 m. Figs 10(a)–(d) show the intrinsic attenuation results calculated from the frequency range of 20–130 Hz for four different trace mixing filters with one, three, five, and seven traces, respectively. On the other hand, Figs 11(a)–(d) show the scattering attenuation results for four different trace mixing filters with one, three, five, and seven traces, respectively. In the estimation of scattering attenuation, the same frequency range as in the case of intrinsic attenuation was applied. However, the scattering attenuation results in the frequency range of from 40 to 80 Hz are shown in Fig. 11. The selection of the frequency range was based on the intensity of the signal in the relative Fourier amplitude spectra in Fig. 9(a). Fig. 12(a) shows a seismic time section across the VSP survey well with geological formation and lithology information. In the seismic section, the tops of several key formations are indicated by the red line with depth information. Figure 8. Open in new tabDownload slide (a) Vertical component of acquired zero-offset VSP data at depths of 1200–3200 m. (b) After the application of first-arrival waveform extraction with a Hanning window (length: 100 ms), first-arrival alignment, a nine-trace median filter, and a seven-trace mixing filter to the waveforms in (a). Figure 8. Open in new tabDownload slide (a) Vertical component of acquired zero-offset VSP data at depths of 1200–3200 m. (b) After the application of first-arrival waveform extraction with a Hanning window (length: 100 ms), first-arrival alignment, a nine-trace median filter, and a seven-trace mixing filter to the waveforms in (a). Figure 9. Open in new tabDownload slide (a) Relative spectral amplitude for the isolated first-arrival waveforms in Fig. 8(b). (b) Centroid frequency as a function of depth (calculated from 20 to 130 Hz) for the isolated first-arrival waveforms in Fig. 8(b). (c) Interval velocity derived from first-arrival time based on zero-offset VSP data. Figure 9. Open in new tabDownload slide (a) Relative spectral amplitude for the isolated first-arrival waveforms in Fig. 8(b). (b) Centroid frequency as a function of depth (calculated from 20 to 130 Hz) for the isolated first-arrival waveforms in Fig. 8(b). (c) Interval velocity derived from first-arrival time based on zero-offset VSP data. Figure 10. Open in new tabDownload slide (a) Intrinsic attenuation results obtained from field zero-offset VSP by the proposed method without the application of a trace mixing filter. (b) Same as (a) except for the application of a trace mixing filter with three traces. (c) Same as (a) except for the application of a trace mixing filter with five traces. (d) Same as (a) except for the application of a trace mixing filter with seven traces. Figure 10. Open in new tabDownload slide (a) Intrinsic attenuation results obtained from field zero-offset VSP by the proposed method without the application of a trace mixing filter. (b) Same as (a) except for the application of a trace mixing filter with three traces. (c) Same as (a) except for the application of a trace mixing filter with five traces. (d) Same as (a) except for the application of a trace mixing filter with seven traces. Figure 11. Open in new tabDownload slide (a) Scattering attenuation results obtained from field zero-offset VSP by the proposed method without the application of a trace mixing filter. (b) Same as (a) except for the application of a trace mixing filter with three traces. (c) Same as (a) except for the application of a trace mixing filter with five traces. (d) Same as (a) except for the application of a trace mixing filter with seven traces. Figure 11. Open in new tabDownload slide (a) Scattering attenuation results obtained from field zero-offset VSP by the proposed method without the application of a trace mixing filter. (b) Same as (a) except for the application of a trace mixing filter with three traces. (c) Same as (a) except for the application of a trace mixing filter with five traces. (d) Same as (a) except for the application of a trace mixing filter with seven traces. Figure 12. Open in new tabDownload slide (a) Seismic time section across the VSP survey well with geological formation and lithology information. The tops of several key formations are indicated by the red line with depth information. (b) Categorization based on three characteristics: seismic facies inferred from the seismic section in (a), intrinsic attenuation inferred from Fig. 10, and scattering attenuation inferred from Fig. 11. Figure 12. Open in new tabDownload slide (a) Seismic time section across the VSP survey well with geological formation and lithology information. The tops of several key formations are indicated by the red line with depth information. (b) Categorization based on three characteristics: seismic facies inferred from the seismic section in (a), intrinsic attenuation inferred from Fig. 10, and scattering attenuation inferred from Fig. 11. In Fig. 11, blue and red indicate positive and negative scattering attenuation, respectively. Negative scattering attenuation indicates that seismic waves gain energy with increasing propagation distance. Mateeva (2003) pointed out that the scattering effects most significantly cause negative attenuation in the attenuation estimation for VSP data. Matsushima et al. (2011) mentioned that negative attenuation is caused by destructive interference of high-frequency components at a propagation distance and constructive interference at a longer propagation distance. In the depth range of from around 1900 to 2200 m, the negative scattering attenuation zone is indicated by the double-headed arrow in Fig. 11. Furthermore, this negative scattering attenuation zone correlates with the zone denoted by ‘constant’ in the centroid frequency profile in Fig. 9(b), where the centroid frequency is maintained approximately constant. On the other hand, in two depth ranges (around 1500–1800 m and around 2300–2800 m), positive scattering attenuation zones are indicated by the double-headed arrows in Fig. 11. In positive scattering attenuation zones, seismic waves easily lose high-frequency components with increasing propagation distance. The positive scattering attenuation zone correlates with the zone denoted by ‘slightly decreasing’ around 1500–1800 m and the zone denoted by ‘decreasing’ around 2300 to 2800 m in the centroid frequency profile in Fig. 9(b), where the centroid frequency slightly or significantly decreases. Furthermore, the positive scattering attenuation zone around 2300–2800 m also correlates with the low to moderate intrinsic attenuation estimates indicated by ‘Small variation’ and ‘Medium variation’ in Fig. 10. Furthermore, the low-scattering attenuation zone in the depth range around 2800–3000 m, as indicated by the double-headed arrow in Fig. 11. This zone corresponds to the transition zone from decreasing to increasing centroid frequency, as indicated by the double-headed arrow in Fig. 9(b). In addition, in Fig. 9(b), several local peaks in the centroid frequency appear, as indicated by the arrow. These peaks correspond to large velocity contrasts, as shown in Fig. 9(c). One possible reason for generating such peaks might be the influence of residual reflection events after wavefield extraction procedures. Finally, in Fig. 12(b), we categorized the VSP survey into five sections based on three characteristics: seismic facies inferred from the seismic section in Fig. 12(a), intrinsic attenuation inferred from Fig. 10, and scattering attenuation inferred from Fig. 11. A discussion of the interpretation of these results is provided below. 5 DISCUSSION 5.1 Limitation of the conventional approach and applicability of the proposed method Although most existing studies assume that the intrinsic and scattering attenuation are independent each other, in the following, we show the limitation of this assumption. As shown in Fig. 4, the results of numerical experiments demonstrate that the shape of the total attenuation profile in Fig. 4(a) is largely different from that of the scattering attenuation in Fig. 4(b), that is, the number of peaks and troughs are different. A possible reason for this is the difference in the effect of multiple reflections, that is, the effect of multiple reflections in the viscoelastic medium with higher attenuation (i.e. Q = 200) is much smaller than that with extremely low attenuation (i.e. Q = 1010). In other words, the mutual interactions between the intrinsic and scattering attenuation are not taken into account in the conventional approach. On the other hand, the intrinsic attenuation values (Fig. 6a) by the application of the proposed method, with the exception of the alternating velocity zone, correspond approximately to the given attenuation value (i.e. 1/Q = 0.005), whereas, in the alternating velocity zone, the intrinsic attenuation estimates exhibit a slight fluctuation but their average corresponds approximately to the given attenuation value, except at the top and bottom of the alternation zone. When a seismic wave is transmitted at a single reflector, transmission loss occurs but its frequency component does not change before or after transmission loss. Matsushima et al. (2016) reported that, in attenuation estimation, the conventional SR method uses the frequency component change between two receiver depths, whereas the MMFS method uses the amplitude change, that is, the MMFS method can detect transmission loss even in the case where there is no frequency component change. Due to the decrease in interference of waves at the top and bottom reflectors in the alternation zone, the degree of frequency component variation at these two reflectors is smaller than for other reflectors within the alternating velocity zone. Therefore, both the top and bottom reflectors yield frequency-independent (or weakly dependent) transmission responses, whereas the reflectors within the alternating velocity zone yield frequency-dependent transmission responses, which depend on the thickness of each layer of the alternation zone. In the alternating velocity zone of the numerical model shown in Fig. 2, the same magnitude of reflection coefficients is alternately changed to positive and negative states at a constant depth interval of 10 m, leading to the constructive interference of sinusoidal waves with the wavelength of 40 m (corresponding frequency of approximately 50 Hz). The significant frequency response of scattering attenuation in the frequency range of around 40–60 Hz, as shown in Fig. 6(b), can be explained by this constructive interference. When we change the thickness of each layer of alternating velocity zone from 10 to 7 m, the significant frequency response of scattering attenuation shifts to a higher frequency (approximately 65–85 Hz) (figure not shown). Based on the above considerations, the alternating velocity zone exhibits the significant frequency dependency, as shown in Fig. 6(b), leading to clear separation between the intrinsic and scattering attenuation, and the top and bottom reflectors in the alternating velocity zone exhibit weak frequency dependency, leading to unclear separation between the intrinsic and scattering attenuation. In Fig. 6, comparing the results obtained with and without a pre-processing filter reveals that the application of a median filter and a trace mixing filter smoothed both the intrinsic and scattering attenuation results. In particular, for the intrinsic attenuation estimates (Fig. 6c), such smoothing effects can reduce the slight fluctuation of intrinsic attenuation estimates observed at the alternating velocity zone in Fig. 6(a), providing slightly fluctuating intrinsic attenuation estimates of the given attenuation values. On the other hand, for the scattering attenuation estimates (Fig. 6d), such smoothing effects reduce the high resolution power of scattering attenuation estimates in Fig. 6(b). Thus, when applying pre-processing filters and determine their parameters, we should consider the optimization of both the resolution power and noise reduction effect. 5.2 Assumption of linear summation of the intrinsic and scattering attenuation As describe above, the total attenuation is widely assumed to be a linear summation of the intrinsic attenuation and the scattering attenuation. However, this assumption does not hold in the numerical experiments of the present study (e.g. Fig. 4), meaning that our numerical models do not satisfy the condition of strong scattering case assumed in Warren (1972). In the numerical experiments, the alternating velocity zone consists of the alternation of base-velocity (2000 m s−1) and anomalous-velocity (2200 m s−1) layers. In order to test the condition in the weak-scattering case, we change the anomalous velocity from 2200 to 2002 m s−1. As a result, the shape of the total attenuation profile becomes much more similar to that of the scattering attenuation (figures not shown) compared with the strong-scattering case shown in Fig. 4. However, the assumption of linear summation of the intrinsic and scattering attenuation does not hold even in the weak scattering case. Although the seismic scattering regime is based on the relationship among the characteristic heterogeneity size, the propagation distance, and the wavelength (Aki & Richards 1980), ‘strong scattering’ has not been clearly defined. Wu & Aki (1988) categorized the scattering attenuation into several domains. The product of the wavenumber k (k = 2π/λ, where λ is the wavelength) and the characteristic heterogeneity length a can represent the effects of heterogeneity size on seismic wave propagation as follows: Quasi-homogeneous regime (ka < 0.01): the medium is assumed to be effectively homogeneous. Rayleigh scattering regime (0.01 < ka < 0.1): the medium is assumed to exhibit single scattering (Born approximation). Mie scattering regime (0.1 < ka < 10): scattering effects are most significant. Forward scattering regime (ka > 10): the medium can be treated by ray theory. In the numerical experiments of the present study, the characteristic heterogeneity length a is 10 m (equivalent to the thickness of each layer in the alternating base-velocity and anomalous-velocity layer zone), the corresponding wavelength is 20–220 m, as calculated from both velocities (2000 and 2200 m s−1), and the frequency range used for attenuation analysis was 10–100 kHz. Thus, the ka values vary from 0.29 to 3.1. Following the criterion of Wu & Aki (1988) described above, the numerical models used herein are classified as being valid in the ‘Mie scattering regime’, where significant scattering may occur. In order to adequately estimate the scattering attenuation from velocity or acoustic impedance data, it is necessary to validate the assumption of strong scattering assumed by Warren (1972). To this end, in the future, we intend to investigate in detail the condition under which the assumption holds. 5.3 Applicability to field data and interpretation We applied the proposed method to real data. In order to investigate the uncertainties of the attenuation results caused by different parameters in pre-processing filters, we applied a trace mixing filter with four different traces, as shown in Figs 10 and 11. For the application of a median filter, although in the present paper we show the attenuation results for the case of using nine traces, we also applied median filters with five and 13 traces and found that the attenuation results from three median filter parameters (five, nine, and 13 traces) exhibit a similar trend. The difference between these results is recognized as being due to the characteristics of the median filter. Mars et al. (1999) reported that a median filter tends to attenuate the desired signal at low frequencies. However, there is trade-off between low-frequency-component preservation and high-frequency-noise elimination, that is, we can preserve the lower frequencies using a longer filter, whereas we cannot eliminate the noise at high frequencies. Mars et al. (1999) also pointed out the importance of amplitude scaling in order to enhance the flattened arrival. Although, in the present study, we did not optimize the amplitude scaling, optimization of the application of a median filter may improve the accuracy of the attenuation results. Similarly, as for waveform extraction using a Hanning window, although in the present paper, we show the attenuation results obtained using a window length of 100 ms, we applied three different window lengths (60, 80, and 100 ms) and found that there were no significant differences between the obtained results. Attenuation estimation methods require extraction of direct-arrival waveforms to isolate these waveforms from unnecessary events such as residual reflection events after applying wavefield separation filters, but such waveform extraction disturbs the attenuation estimation due to spectral leakage, which is defined as the leaking of energy from its true frequency component into other frequency components, distorting the spectral distribution (Matsushima et al.2014). Therefore, the difference in attenuation results between different window lengths can be attributed to the effects of spectral leakage or residual reflection events. Nevertheless, the selection of parameters for pre-processing does not largely change the attenuation results, as shown in Figs 10 and 11. Based on the results for the both intrinsic and scattering attenuation shown in Figs 10 and 11, and by taking into account the geological formation, lithology, and seismic facies, we categorized the results into five groups, as shown in Fig. 12: Depth range from 1200 to 1431 m: Based on the observed seismic facies (semi-parallel to parallel, continuous), we can assume horizontal layering, causing scattering attenuation due to intrabed multiples. By carefully comparing Figs 10 and 11 in this depth range, the variation pattern of the scattering attenuation correlates well with that of the intrinsic attenuation. In general, the depth intervals in horizontal layering are not constant, but rather are random, leading to less frequency dependency of scattering attenuation. Consequently, such weakly frequency-dependent scattering attenuation can be recognized as intrinsic attenuation. On the other hand, the strong scattering attenuation observed in Fig. 11 may be caused by the local frequency dependency associated with local layering. Depth range from 1431 to 1940 m: Although we can see high values of both the intrinsic and scattering attenuation, unlike the case of (i), their variation patterns do not correlate well with each other. Furthermore, the scattering attenuation reveals the dominance of the positive values. Seismic facies exhibit discontinuous and semi-chaotic reflection features, meaning that geological heterogeneities make seismic propagation complicated and consequently contaminate the seismic reflection section. Huang et al. (2009) attributed the strong scattering attenuation to small-scale geological heterogeneities in the horizontal direction (ka is nearly equal to 10). They constructed heterogeneous petrophysical models to investigate the effects of horizontal discontinuities on the scattering attenuation, and finally pointed out the importance of horizontal heterogeneities. Through detailed analyses of core data and FMI/FMS (Formation Micro Imager/Formation Micro Scanner) logs, Ali & Worthington (2011) mentioned that the Simsima reservoir (upper part of this depth range) is highly fractured. Furthermore, by combining well data (cores, FMI and sonic logs) and curvature analysis using 3-D seismic data from the top Simsima reservoir, BEICIP-FRANLAB (2002) indicated that there are fracture swarms or sub-seismic faults with very few fractures outside the clusters, and the width of these fractured zones range from 50–100 m and are spaced 100–400 m apart. Assuming that the dominant wavelength of seismic data is around 50 m, such heterogeneities caused by fracture swarms or sub-seismic faults are categorized as the Mie scattering regime. Therefore, in this group, geological heterogeneities in the horizontal direction may be more dominant than in other groups. Depth range from 1940 to 2308 m: Similar to the case of (ii), although we can see high values of both the intrinsic and scattering attenuation, their variation patterns do not correlate well with each other. Thus, the intrinsic attenuation may be considered to be associated with the physical properties of each layer. Furthermore, the scattering attenuation is primarily negative. As mentioned above, negative scattering attenuation indicates that seismic waves gain energy with increasing propagation distance. Although there are many possible factors that cause negative attenuation, one possible factor is the superposition of the forward scattering of high-frequency components on direct arrivals. According to the classification by Wu & Aki (1988), the forward scattering regime is dominant when ka > 10, that is, the heterogeneity size can be considered to be larger than that in the case of (ii). If the heterogeneity size is larger than that in the case of (ii), for which heterogeneities are significant and cause large fluctuations in seismic reflected data, such larger heterogeneities in this group cause smaller fluctuations in seismic reflected data and enhance the tendency of forward scattering, which may result in negative attenuation. The combination of negative scattering attenuation and high (positive) intrinsic attenuation may cause the constant trend in the centroid frequency profile in Fig. 9(b). Depth range from 2308 to 2733 m: This group exhibits weak reflection in the seismic section, low intrinsic attenuation with moderate variation, and moderate positive scattering attenuation. As shown in Fig. 9(b), the centroid frequency significantly decreases within this group, implying that small-scale heterogeneities may play an important role in causing scattering attenuation. Although it is widely believed that seismic events caused by small-scale heterogeneities such as fractures cannot be detected in the seismic frequency range, the proposed method may be able to indirectly estimate small-scale heterogeneities (e.g. fractured zone). Depth range from 2733 to 3051 m: This group exhibits weak reflection in the seismic section, low intrinsic attenuation with low variation, and low scattering attenuation. Based on this observation, this group may be interpreted as a homogeneous zone. Lubis & Harith (2014) mentioned that the significant diagenesis process and complex depositional environment makes pore systems in carbonates far more complicated than in clastics, leading to the variations in acoustic impedance between layers. Ali & Worthington (2011) mentioned that carbonate reservoirs contain open or partially open fractures with a wide range of scales, varying spatial distributions and hydraulic properties. When the wavelength is much larger than the size of fractures, the medium is assumed to be effectively homogeneous (quasi-homogeneous regime), and its acoustic impedance becomes lower than without fractures. Decker et al. (2015) pointed out that high acoustic impedance contrasts typically exist between carbonate structures and surrounding rocks, leading to a strong reflective interface that can generate multiple reflections. In this case, the mechanism of scattering attenuation is due to such intrabed multiples. Bouchaala et al. (2016a) mentioned that the scattering attenuation estimated in carbonate reservoirs is much higher than that in siliciclastic media. Furthermore, as described above, fracture swarms cause various scales of heterogeneities relative to the wavelength, leading to various scattering regimes. On the other hand, it is widely believed that intrinsic attenuation is caused by relative motion between viscous fluid and rock framework. Johnston et al. (1979) illustrated several proposed intrinsic attenuation mechanisms such as crack lubrication facilitating friction, Biot fluid flow, and squirting flow. Furthermore, Winkler & Nur (1982) proposed an inter-crack flow model to explain the reason why shear wave attenuation is larger than compressional wave attenuation. Bouchaala et al. (2016a) estimated both compressional and shear attenuation by using sonic logging waveform data acquired in a carbonate reservoir located in Abu Dhabi, United Arab Emirates. They observed that shear wave attenuation is larger than compressional wave attenuation, implying the intercrack flow model as a possible attenuation mechanism. Further investigation of possible attenuation mechanisms by using sonic logging waveform data and laboratory measurements is needed. 5.4 Limitations of the proposed method and future research The proposed method does not require sonic and density well-log data or the assumption of a 1-D layered medium and can simultaneously derive both intrinsic and scattering attenuation estimates from only VSP data, whereas the conventional approach uses well-log velocities and densities to calculate scattering attenuation caused by multiples. As one potential limitation, the proposed method is based on the assumption that intrinsic and scattering attenuation are frequency independent and frequency dependent, respectively. As for the assumption of frequency-independent attenuation, Patton (1988) pointed out the validity of a frequency-independent Q in a narrow frequency range. Note that the attenuation is implicitly assumed to consist of only intrinsic attenuation. Nevertheless, if the intrinsic attenuation is dependent on frequency in the frequency range of interest, the proposed method will cause an error in attenuation estimates. On the other hand, the scattering effect including multiple reflections in a layered medium is generally dependent on frequency, as described by Wu & Aki (1988). O'Doherty & Anstey (1971), reported that the transmission response is critically dependent on the nature of the stratigraphy. However, as mentioned above, transmission loss at a single reflector, which should be categorized as scattering attenuation, is independent of frequency and thus is categorized as intrinsic attenuation. In the future, we intend to perform numerical modelling to quantitatively estimate the effects of heterogeneities on seismic reflection sections, intrinsic attenuation, and scattering attenuation, allowing us to conduct quantitative interpretation. Another possible area for future research is to apply the proposed method to sonic logging waveform data. Suzuki & Matsushima (2013) reported that the spatial sampling of the velocity logging data at 0.15 m intervals is not sufficiently fine to characterize heterogeneities when adopting the conventional approach using well-log velocities and densities to calculate scattering attenuation caused by multiples. In this case, the proposed method, which does not require well-log velocities and densities, has the potential to provide high resolution power of scattering attenuation for characterizing small-scale heterogeneities such as fractures. 6 CONCLUSIONS We proposed a novel method for separating intrinsic and scattering attenuation for zero-offset VSP data by modifying the SI+MMFS method based on the assumption that intrinsic attenuation is frequency independent and scattering attenuation is frequency dependent. The proposed method can simultaneously derive both intrinsic and scattering attenuation estimates from only VSP data without sonic and density well-log data. The numerical results obtained herein demonstrate the limitation of the conventional approach, which uses well-log velocities and densities to calculate scattering attenuation based on the assumption that the intrinsic and scattering attenuation are independent each other. The numerical results also demonstrated the superiority of the proposed method as compared to the conventional approach and the importance of optimizing parameters in the application of pre-processing filters to balance the resolution power and noise reduction effect. Finally, we applied the proposed method to real zero-offset VSP data acquired at an onshore Abu Dhabi oil field to conduct an investigation of the obtained attenuation results through correlation with seismic facies and discussed the limitations of the proposed method. Acknowledgements We would like to thank one anonymous reviewer for his/her valuable and constructive comments that were helpful to greatly improve the quality of this paper. We also would like to thank the Oil-Subcommittee of the Abu Dhabi National Oil Company (ADNOC) and its operating companies for sponsoring this project and providing data used in this study. REFERENCES Aki K. , Richards, P.G., 1980 . Quantitative Seismology: Theory and Methods , W.H. Freeman & Co . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Ali M.Y. , Worthington M.H., 2011 . Seismic modelling of a fractured carbonate reservoir in Abu Dhabi, United Arab Emirates , GeoArabia , 16 89 – 106 . OpenURL Placeholder Text WorldCat Barton N. , 2006 . Rock Quality, Seismic Velocity, Attenuation and Anisotropy , Taylor & Francis . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC BEICIP-FRANLAB , 2002 . Update of the Fracture Model and Fracture Permeability Computation on Shah Field , BEICIP-FRANLAB . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bouchaala F. , Ali M.Y., Matsushima J., 2016a . Estimation of seismic attenuation in carbonate rocks using three different methods: application on VSP data from Abu Dhabi oilfield , J. Appl. Geophys. , 129 79 – 91 . Google Scholar Crossref Search ADS WorldCat Bouchaala F. , Ali M.Y., Matsushima J., 2016b . Attenuation modes from vertical seismic profiling and sonic waveform in a carbonate reservoir, Abu Dhabi, United Arab Emirates , Geophys. Prospect. , 64 1030 – 1047 . Google Scholar Crossref Search ADS WorldCat Burns D.R. , Willis M.E., Toksöz M.N., Vetri L., 2007 . Fracture properties from seismic scattering , Leading Edge , 26 ( 9 ), 1186 – 1196 . Google Scholar Crossref Search ADS WorldCat Decker L. , Janson X., Fomel S., 2015 . Carbonate reservoir characterization using seismic diffraction imaging , Interpretation , 3 ( 1 ), SF21 – SF30 . Google Scholar Crossref Search ADS WorldCat Fang X.D. , Fehler M.C., Zhu Z.Y., Zheng Y.C., Burns D.R., 2014 . Reservoir fracture characterization from seismic scattered waves , Geophys. J. Int. 196 ( 1 ), 481 – 492 . Google Scholar Crossref Search ADS WorldCat Frazer L.N. , Sun X.H., Wilkens R.H., 1997 . Changes in attenuation with depth in an ocean carbonate section: Ocean Drilling Program sites 806 and 807, Ontong Java Plateau , J. geophys. Res. , 102 2983 – 2997 . Google Scholar Crossref Search ADS WorldCat Ganley D.C. , Kanasewich E.R., 1980 . Measurement of absorption and dispersion from check shot surveys , J. geophys. Res. , 85 5219 – 5226 . Google Scholar Crossref Search ADS WorldCat Gist G.A. , 1994 . Seismic attenuation from 3-D heterogeneities: A possible resolution of the VSP attenuation paradox , in 64th Ann. Internat. Mtg., Soc. of Expl. Geophys ., pp. 1042 – 1045 . OpenURL Placeholder Text WorldCat Huang J.W. , Bellefleur G., Milkereit B., 2009 . Seismic modeling of multidimensional heterogeneity scales of Mallik gas hydrate reservoirs, Northwest Territories of Canada , J. geophys. Res. , 114 B07306 , doi:10.1029/2008JB006172 . Google Scholar Crossref Search ADS WorldCat Johnston D.H. , Toksöz M.N., Timur A., 1979 . Attenuation of seismic waves in dry and saturated rocks: II. Mechanisms , Geophysics , 44 691 – 711 . Google Scholar Crossref Search ADS WorldCat Lerche I. , Menke W., 1986 . An inversion method for separating apparent and intrinsic attenuation in layered media , Geophys. J. R. astron. Soc., 87 333 – 347 . OpenURL Placeholder Text WorldCat Lubis L.A. , Harith Z.Z.T., 2014 . Pore type classification on carbonate reservoir in offshore Sarawak using rock physics model and rock digital images , Earth Environ. Sci. , 19 ( 1 ), doi:10.1088/1755-1315/19/1/012003 . OpenURL Placeholder Text WorldCat Maresh J. , White R.S., Hobbs R.W., Smallwood J.R., 2006 . Seismic attenuation of Atlantic margin basalts: observations and modeling , Geophysics , 71 B211 – B221 . Google Scholar Crossref Search ADS WorldCat Mars J. , Rector W.J. III, Lazaratos S.K., 1999 . Filter formulation and wavefield separation of cross-well seismic data , Geophys. Prospect. , 47 ( 5) , 611 – 636 . Google Scholar Crossref Search ADS WorldCat Mateeva A. , 2003 . Thin horizontal layering as a stratigraphic filter in absorption estimation and seismic deconvolution , PhD thesis , 108 ., Colorado School of Mines , Colorado . OpenURL Placeholder Text WorldCat Matsushima J. , 2006 . Seismic wave attenuation in methane hydrate-bearing sediments: Vertical seismic profiling data from the Nankai Trough exploratory well, offshore Tokai, central Japan , J. geophys. Res. , 111 B10101 , doi:10.1029/2005JB004031 . Google Scholar Crossref Search ADS WorldCat Matsushima J. , Suzuki M., Kato Y., Rokugawa S., 2011 . Estimation of ultrasonic scattering attenuation in partially frozen brines using magnetic resonance images , Geophysics , 76 ( 1 ), T13 – T25 . Google Scholar Crossref Search ADS WorldCat Matsushima J. , Suzuki M., Matsugi I., Kato Y., Rokugawa S., 2014 . Attenuation estimation using sweep signals in ultrasonic laboratory measurements , Geophysics , 79 ( 3 ), V115 – V130 . Google Scholar Crossref Search ADS WorldCat Matsushima J. , Ali M.Y., Bouchaala F., 2016 . Seismic attenuation estimation from zero-offset VSP data using seismic interferometry , Geophys. J. Int. 204 1288 – 1307 . Google Scholar Crossref Search ADS WorldCat McDonal F.J. , Angona F.A., Mills R.L., Sengbush R.L., Van Nostrand R.G., White J.E., 1958 . Attenuation of shear and compressional waves in Pierre shalt , Geophysics , 23 421 – 439 . Google Scholar Crossref Search ADS WorldCat Mehta K. , Snieder R., Graizer V., 2007 . Extraction of near-surface properties for a lossy layered medium using the propagator matrix , Geophys. J. Int. 169 271 – 280 . Google Scholar Crossref Search ADS WorldCat O'Doherty R.F. , Anstey N.A., 1971 . Reflections on amplitudes , Geophys. Prospect. , 19 430 – 458 . Google Scholar Crossref Search ADS WorldCat Patton S.W. , 1988 . Robust and least-squares estimation of acoustic attenuation from well-log data , Geophysics , 53 1225 – 1232 . Google Scholar Crossref Search ADS WorldCat Shaw F. , Worthington M.H., White R.S., Andersen M.S., Petersen U.K., 2008 . Seismic attenuation in Faroe Islands basalts , Geophys. Prospect. , 56 5 – 20 . Google Scholar Crossref Search ADS WorldCat Schoenberger M. , Levin F.K., 1974 . Apparent attenuation due to intrabed multiples , Geophysics , 39 278 – 291 . Google Scholar Crossref Search ADS WorldCat Schoenberger M. , Levin F.K., 1978 . Apparent attenuation due to intrabed multiples, II , Geophysics , 43 730 – 737 . Google Scholar Crossref Search ADS WorldCat Shapiro S.A. , Zien H., Hubral E., 1994 . A generalized O'DohertyAnstey formula for waves in finely-layered media , Geophysics , 59 1750 – 1762 . Google Scholar Crossref Search ADS WorldCat Spencer T.W. , Sonnad J.R., Butler T.M., 1982 . Seismic Q—Stratigraphy or dissipation , Geophysics , 47 16 – 24 . Google Scholar Crossref Search ADS WorldCat Suzuki H. , Matsushima J., 2013 . Quantifying uncertainties in attenuation estimation at methane-hydrate-bearing zones using sonic waveform logs , Geophysics , 78 ( 5 ), D339 – D353 . Google Scholar Crossref Search ADS WorldCat Vasconcelos I. , Snieder R., 2008 . Interferometry by deconvolution, Part 1: Theory for acoustic waves and numerical examples , Geophysics , 73 ( 3 ), S115 – S128 . Google Scholar Crossref Search ADS WorldCat Vlastos S. , Liu E., Main I.G., Li X.Y., 2003 . Numerical simulation of wave propagation in media with discrete distributions of fractures: effects of fracture sizes and spatial distributions , Geophys. J. Int. 152 ( 3 ), 649 – 668 . Google Scholar Crossref Search ADS WorldCat Vlastos S. , Liu E., Main I.G., Schoenberg M., Narteau C., Li X.Y., Maillot B., 2006 . Dual simulations of fluid flow and seismic wave propagation in a fractured network: effects of pore pressure on seismic signature , Geophys. J. Int. 166 ( 2 ), 825 – 838 . Google Scholar Crossref Search ADS WorldCat Vlastos S. , Liu E., Main I.G., Narteau C., 2007 . Numerical simulation of wave propagation in 2-D fractured media: scattering attenuation at different stages of the growth of a fracture population , Geophys. J. Int. 171 ( 2 ), 865 – 880 . Google Scholar Crossref Search ADS WorldCat Walsh J.B. , 1969 . New analysis of attenuation in partially melted rock , J. geophys. Res. , 74 4333 – 4337 . Google Scholar Crossref Search ADS WorldCat Warren N. , 1972 . Q and structure , Moon 4 430 – 441 . Google Scholar Crossref Search ADS WorldCat Wegler U. , 2004 . Diffusion of seismic waves in a thick layer: theory and application to Vesuvius volcano , J. geophys. Res. , 109 B07303 , doi:10.1029/2004JB003048 . Google Scholar Crossref Search ADS WorldCat Willis M. , Burns D., Rao R., Minsley B., Toksoz M., Vetri L., 2006 . Spatial orientation and distribution of reservoir fractures from scattered seismic energy , Geophysics , 71 ( 5 ), O43 – O51 . Google Scholar Crossref Search ADS WorldCat Winkler K. , Nur A., 1982 . Seismic attenuation—effects of pore fluids and frictional sliding , Geophysics , 47 1 – 15 . Google Scholar Crossref Search ADS WorldCat Wu R. , Aki K., 1988 . Seismic wave scattering in three-dimensionally heterogeneous earth—introduction , Pure appl. Geophys. , 128 1 – 6 . Google Scholar Crossref Search ADS WorldCat Yang Y. , Li Y., Liu T., 2009 . 1D viscoelastic waveform inversion for Q structures from the surface seismic and zero-offset VSP data , Geophysics , 74 ( 6 ), WCC141 – WCC148 . Google Scholar Crossref Search ADS WorldCat © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - A novel method for separating intrinsic and scattering attenuation for zero-offset vertical seismic profiling data
JF - Geophysical Journal International
DO - 10.1093/gji/ggx391
DA - 2017-12-01
UR - https://www.deepdyve.com/lp/oxford-university-press/a-novel-method-for-separating-intrinsic-and-scattering-attenuation-for-5VatAsuiWc
SP - 1655
VL - 211
IS - 3
DP - DeepDyve
ER -