A Borehole Acoustic Calculation Approach with Gas Hydrate Saturation Inversion in Gas Hydrate-Bearing Sediments

Abstract

The inversion of gas hydrate saturation is a critical procedure in the evaluation of hydrate reservoirs. In this paper, a theoretical model for a borehole acoustic wavefield excited by multipole sources is established for the first time. On this basis, the attenuation of the dipole flexural waves is obtained, and in combination with the results of sensitivity analysis, an approach for inverting natural gas hydrates using the attenuation characteristics of the dipole flexural wave is proposed. The results of the sensitivity analysis demonstrate that the attenuation of the dipole flexural wave is sensitive to gas hydrate saturation. Numerical results derived from synthetic logging data are provided to illustrate the viability of the inversion of gas hydrate saturation. Even when significant noise is introduced into the receiver signal arrays, the inversion method remains stable and accurately assesses gas hydrate saturation. The correctness and effectiveness of the proposed approach are substantiated through the processing of numerical simulation data. This work provides a potent processing approach for evaluating reservoir hydrate saturation utilizing acoustic well-logging information.

1. Introduction

In recent years, the focus of oil and gas exploration and development has transitioned from conventional reservoirs to complex unconventional ones [1,2]. The discovery of gas hydrates in permafrost regions of polar continental shelves and beneath seafloor sediments has piqued great interest. While opinions may vary regarding the concentration of gas hydrates in hydrate-bearing sediments and the overall methane content within gas hydrate accumulations, there is a consensus that the global methane reserves in gas hydrate deposits possess tremendous energy potential exceeding that of traditional oil, gas, and coal reserves [3,4,5,6]. Thus, gas hydrate epitomizes an unconventional energy asset poised to be significant in the forthcoming decades [7,8,9,10,11,12,13].

Geophysical well-logging is pivotal in pinpointing and quantifying gas hydrate reservoirs along with their constituents. The acoustic method is an important method of studying reservoir characteristics. Acoustic logging can proficiently measure the elastic parameters, porosity, and saturation of adjacent formations and is an indispensable technique for evaluating gas hydrate-bearing sediments. Gas hydrate-bearing sediments are three-phase porous media composed of a solid frame saturated with a fluid (usually water) and hydrates. The Carcione–Leclaire theory [14,15,16] is important for studying acoustic wave propagation in such a three-phase porous medium. The theory predicts three compressional and two shear waves produced when an acoustic wave propagates in a three-phase porous medium such as a gas hydrate-bearing sediment.

Currently, there is limited research on the acoustic wavefield within three-phase porous media in borehole environments. Due to the small difference between the compressional wave velocities of such a three-phase porous medium and the equivalent two-phase medium, scholars primarily utilize the compressional wave velocity obtained from acoustic logging to estimate the saturation of natural gas hydrates on the basis of the equivalent two-phase porous medium model. The saturation of gas hydrate in hydrate-bearing sediments was estimated by Shankar and Riedel using a combination of compressional wave velocity and resistivity logging [17]. Gao et al. [18] applied the improved Biot–Gassmann model proposed by Lee [19] to invert the saturation of gas hydrate in the hydrate-bearing reservoir using the acoustic wave velocity. Lin et al. [20] established a relationship between gas hydrate saturation and geological velocity based on the particle contact model [21] and the simplified three-phase equation [22]. Zhang et al. [23] studied the velocity dispersion and attenuation characteristics of acoustic waves in gas hydrate-bearing sediments using the effective medium theory [24] and used this model to predict the gas hydrate content in three drilling sites in the Shenhu area of the South China Sea. The authors of this paper [25] estimated gas hydrate saturation based on the velocity of the first kind of compressional wave of the Carcione–Leclaire three-phase Biot-type theory [26,27].

The aforementioned studies indicate that numerous scholars have researched the characteristics of gas hydrate reservoirs based on different models and utilized various models to estimate the hydrate saturation. The existence of gas hydrates essentially affects the physical properties of sediments. These alterations encompass key parameters such as permeability, porosity, and various geophysical characteristics including velocity, attenuation, and resistivity. [28]. Within these geophysical anomalies, the velocity and attenuation of both P- and S-waves have played a crucial role in characterizing the distribution and occurrence of hydrates. Furthermore, they have been instrumental in evaluating the hydrate saturation in gas hydrate-bearing sediments [29,30,31,32,33]. Logging curves exhibit large guide wave amplitudes, making them easy to identify, and they can, to some extent, reflect the reservoir’s characteristics. However, the attenuation characteristics of guide waves in a borehole of three-phase porous media differ significantly from those in two-phase media. The lack of research on the borehole acoustic field of multiphase porous media hampers the effective utilization of guide wave information. The authors of this paper [34] studied the monopole borehole acoustic wavefield of a three-phase porous medium using the staggered-grid finite-difference time-domain (FDTD) scheme and the real axis integration (RAI) algorithm.

Due to the inability of traditional monopoles to receive shear head wave signals in slow formations, people usually use dipole sources to invert formation information [35,36,37]. This technology has not been limited to measuring only shear wave velocity but has also seen wider application and development. Its rapid development and extensive application are set to make significant contributions to hydrocarbon exploration in the field of acoustic logging. Currently, dipole acoustic logging is widely used in two-phase porous media. Nevertheless, the intricate nature of three-phase porous media poses significant challenges for research. Studies are scarce on the acoustic fields in a borehole and the responses elicited by multipole sources. This research gap hinders the potential application of multipole logging techniques for comprehensive reservoir evaluation.

In this paper, we present a representation approach for the borehole acoustic wavefield within three-phase porous media excited by multipole sources. Combining the acoustic potential function of the borehole and natural gas hydrate-bearing sediment with the boundary conditions of the borehole wall, the borehole acoustic waveforms can be obtained in the calculation. Further, the attenuation of the dipole flexural wave can be calculated. Building on this foundation, we propose an innovative methodology for estimating the saturation of natural gas hydrates in reservoirs using dipole flexural wave attenuation. The theoretical analysis investigates the sensitivity of the dipole flexural wave to gas hydrate saturation in a natural gas hydrate-bearing sediment (a three-phase porous medium). Based on this analysis, flexural wave attenuation is employed to inverse the saturation of natural gas hydrates. The results demonstrate that this method can effectively estimate the saturation of natural gas hydrates.

2. Theoretical Model of the Borehole Acoustic Wavefield

Gas hydrate-bearing sediments represent a three-phase porous medium in which solid grains coexist with gas hydrates and water within the pore spaces. Drawing upon Biot’s foundational theory, both Leclaire [14] and Carcione et al. [15] put forth a percolation theory, delving into the characteristics of acoustic wave propagation in such three-phase porous media. Their theory posits that, when an acoustic wave propagates in a three-phase porous medium, three compressional (P1, P2 and P3) and two shear (S1 and S2) waves are produced. Here, we delve into the analysis of the acoustic field within a borehole encompassed by such a porous medium.

The borehole wavefield in the frequency–wavenumber domain can be formulated as

$$\phi (k_z,\omega )=\frac{iF}{n!{\rho }_f{\omega }^2}\left[{\epsilon }_n{K}_n\left({\alpha }_{f}r\right)\right.\left.+A_f(k_z,\omega )I_n\left({\alpha }_{f}r\right)\right]\cos n(\theta -{\theta }_0),$$

(1)

where F symbolizes the frequency spectrum of the source. ${\epsilon }_n$ refers to the Neumann coefficient. The value of ${\epsilon }_n$ is equal to one when n is zero; when n is not zero, ${\epsilon }_n$ is equal to two. ${\alpha }_f$ pertains to the radial wavenumber associated with borehole fluid, r refers to the borehole radius, and $A_f$ is an undetermined coefficient, interpreted as the reflection coefficient of the borehole wall. $K_n$ and $I_n$ are the modified Bessel functions, and the term “n” designates an n-th order sound source; in dipole acoustic logging, n = 1. In addition, $i=\sqrt{-1}$. The genuine time-domain waveform can be retrieved by employing a double Fourier transform. Both the dispersion and attenuation of mode waves hinge on the poles of $A_f(k_z,\omega )$ in the complex number domain.

In a three-phase porous medium, the acoustic wavefield is represented by the displacements of the solid grain frame ${\overrightarrow{u}}_1$, pore fluid ${\overrightarrow{u}}_2$, and gas hydrate ${\overrightarrow{u}}_3$. The displacements induced by monopole sources can be written as follows:

$$\begin{array}{l}{\overrightarrow{u}}_1=\nabla {\phi }_1+\nabla \times \nabla \times ({\eta }_1\widehat{z})+\nabla \times ({\chi }_1\widehat{z}),\\ {\overrightarrow{u}}_2=\nabla {\phi }_2+\nabla \times \nabla \times ({\eta }_2\widehat{z})+\nabla \times ({\chi }_2\widehat{z}),\\ {\overrightarrow{u}}_3=\nabla {\phi }_3+\nabla \times \nabla \times ({\eta }_3\widehat{z})+\nabla \times ({\chi }_3\widehat{z}),\end{array}$$

(2)

where ${\phi }_m$, ${\eta }_m$, and ${\chi }_m$ ($m$ = 1, 2, 3) denote the displacement potentials for compressional (P1, P2, and P3), SV (SV1, SV2) and SH (SH1 and SH2) waves, respectively. The displacement potentials for the compressional waves of each phase are segmented into three distinct parts, corresponding to the P1, P2, and P3 waves.

In the frequency–wavenumber domain, the displacement potentials for these compressional waves can be articulated as

$$\begin{array}{l}{\phi }_1(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[B_1(k_z,\omega )K_n({\alpha }_{\mathrm{P}1}r)+C_1(k_z,\omega )K_n({\alpha }_{\mathrm{P}2}r)+D_1(k_z,\omega )K_n({\alpha }_{\mathrm{P}3}r)]\cos n(\theta -{\theta }_0),\\ {\phi }_2(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[B_2(k_z,\omega )K_n({\alpha }_{\mathrm{P}1}r)+C_2(k_z,\omega )K_n({\alpha }_{\mathrm{P}2}r)+D_2(k_z,\omega )K_n({\alpha }_{\mathrm{P}3}r)]\cos n(\theta -{\theta }_0),\\ {\phi }_3(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[B_3(k_z,\omega )K_n({\alpha }_{\mathrm{P}1}r)+C_3(k_z,\omega )K_n({\alpha }_{\mathrm{P}2}r)+D_3(k_z,\omega )K_n({\alpha }_{\mathrm{P}3}r)]\cos n(\theta -{\theta }_0),\end{array}$$

(3)

where $B_m$, $C_m$, and $D_m$ are the unknown amplitudes of three (one fast and two slow) compressional waves, respectively; ${\alpha }_{\mathrm{P}1}$, ${\alpha }_{\mathrm{P}2}$, and ${\alpha }_{\mathrm{P}3}$ are the radial wavenumbers of these compressional waves. In addition, the displacement potentials of SV waves in the frequency–wavenumber domain can be expressed as

$$\begin{array}{l}{\eta }_1(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[E_1(k_z,\omega )K_n({\alpha }_{S\mathrm{V}1}r)+F_1(k_z,\omega )K_n({\alpha }_{\mathrm{SV}2}r)]\cos n(\theta -{\theta }_0),\\ {\eta }_2(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[E_2(k_z,\omega )K_n({\alpha }_{\mathrm{SV}1}r)+F_2(k_z,\omega )K_n({\alpha }_{\mathrm{SV}2}r)]\cos n(\theta -{\theta }_0),\\ {\eta }_3(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[E_3(k_z,\omega )K_n({\alpha }_{\mathrm{SV}1}r)+F_3(k_z,\omega )K_n({\alpha }_{\mathrm{SV}2}r)]\cos n(\theta -{\theta }_0),\end{array}$$

(4)

where $E_m$ and $F_m$ represent two unknown amplitudes of the two kinds of SV waves, and ${\alpha }_{\mathrm{SV}1}$ and ${\alpha }_{\mathrm{SV}2}$ depict the radial wavenumbers corresponding to these SV waves. The displacement potentials of SH waves can be expressed as

$$\begin{array}{l}{\chi }_1(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[G_1(k_z,\omega )K_n({\alpha }_{\mathrm{SH}1}r)+H_1(k_z,\omega )K_n({\alpha }_{\mathrm{SH}2}r)]\sin n(\theta -{\theta }_0),\\ {\chi }_2(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[G_2(k_z,\omega )K_n({\alpha }_{\mathrm{SH}1}r)+H_2(k_z,\omega )K_n({\alpha }_{\mathrm{SH}2}r)]\sin n(\theta -{\theta }_0),\\ {\chi }_3(r,k_z,\omega )=\frac{iF}{2{\rho }_f{\omega }^2}[G_3(k_z,\omega )K_n({\alpha }_{\mathrm{SH}1}r)+H_3(k_z,\omega )K_n({\alpha }_{\mathrm{SH}2}r)]\sin n(\theta -{\theta }_0),\end{array}$$

(5)

where $G_m$ and $H_m$ represent two unknown amplitudes of the SH waves, while ${\alpha }_{\mathrm{SH}1}$ and ${\alpha }_{\mathrm{SH}2}$ depict the radial wavenumbers corresponding to these SH waves. Among the mentioned formulas, $A_f$, $B_m$, $C_m$, $D_m$, $E_m$, $F_m$, $G_m$, and $H_m$ function as the undetermined coefficients. Their precise values are ascertained by adhering to the borehole wall’s boundary conditions. Assuming a connection exists between the displacements of compressional waves for each phase:

$${\phi }_2=l_1{\phi }_1=l_2{\phi }_3,$$

(6)

subsequently, we can acquire

$$\begin{array}{l}{B}_2=l_{11}{B}_1,C_2=l_{12}{C}_1,D_2=l_{13}{D}_1,\\ B_3=\frac{l_{11}}{l_{21}}{B}_1,C_3=\frac{l_{12}}{l_{22}}{C}_1,D_3=\frac{l_{13}}{l_{23}}{D}_1.\end{array}$$

(7)

Assuming a correlation exists between the displacements of SV waves for each phase:

$${\eta }_2=l_3{\eta }_1=l_4{\eta }_3,$$

(8)

then we have

$$E_2=l_{31}{E}_1, F_2=l_{32}{F}_1, E_3=\frac{l_{31}}{l_{41}}{E}_1, F_3=\frac{l_{32}}{l_{42}}{F}_1.$$

(9)

Assuming a correlation exists between the displacement of SH waves for each phase:

$${\chi }_2=l_5{\chi }_1=l_6{\chi }_3,$$

(10)

then we can get

$$G_2=l_{51}{G}_1, H_2=l_{52}{H}_1, G_3=\frac{l_{51}}{l_{61}}{G}_1,{ \mathrm{H}}_3=\frac{l_{52}}{l_{62}}{H}_1.$$

(11)

The focus of this paper is on exploring the open-hole boundary conditions, which can be formulated as:

$$\begin{array}{l}{u}_{rfin}={\varphi }_s{u}_{rs}+{\varphi }_f{u}_{rf}+{\varphi }_h{u}_{rh},\\ -{\varphi }_s{P}_{fin}={\sigma }_{rrs},\\ -{\varphi }_h{P}_{fin}={\sigma }_{rrh},\\ P_{fin}=P_f,\\ {\tau }_{rzs}=0,\\ {\tau }_{rzh}=0,\\ {\tau }_{r\theta s}=0,\\ {\tau }_{r\theta h}=0,\end{array}$$

(12)

where the subscripts s, f, and h correspond to solid grains, pore fluid, and hydrate, respectively. The subscripts in and out differentiate between the regions within and beyond the borehole wall. The first part of Equation (12) describes the continuous normal displacement of fluid both inside and outside the well. The next two equations indicate that the pressure of the fluid within the well, perpendicular to the borehole wall, is equal to the normal stress exerted by the solid grains and natural gas hydrates in the gas hydrate-bearing sediment. The fourth equation indicates that the pressure of the fluid within the well, in the vertical borehole direction, must be equivalent to the pressure of the fluid in the porous formation outside the well in that direction. The next two equations express the continuity of shear stresses in the rz direction. The seventh and eighth equations assert the continuity of shear stresses in the $r\theta$ direction. The specific expressions in the above equations can be expressed by

$$k_f=\frac{\omega }{v_f},$$

(13)

$${\alpha }_f=\sqrt{k_z^2-k_f^2},$$

(14)

$${\alpha }_{\mathrm{Pi}}=\sqrt{k_z^2-k_{\mathrm{Pi}}^2}, \mathrm{i}=1, 2, 3,$$

(15)

$${\alpha }_{\mathrm{SVj}}=\sqrt{k_z^2-k_{\mathrm{SVj}}^2}, \mathrm{j}=1, 2,$$

(16)

$${\alpha }_{\mathrm{SHk}}=\sqrt{k_z^2-k_{\mathrm{SHk}}^2}, \mathrm{k}=1, 2.$$

(17)

$${\alpha }_{\mathrm{SHk}}=\sqrt{k_z^2-k_{\mathrm{SHk}}^2}, \mathrm{k}=1, 2.$$

(18)

$$\begin{array}{l}{\sigma }_{rrs}=\left(K_1-\frac{2}{3}{\mu }_1\right)(-k_{\mathrm{P}1}^2{B}_1{K}_n\left({\alpha }_{\mathrm{P}1}a\right)-k_{\mathrm{P}2}^2{C}_1{K}_n\left({\alpha }_{\mathrm{P}2}a\right)\\ -k_{\mathrm{P}3}^2{D}_1{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)+C_{12}(-k_{\mathrm{P}1}^2{B}_1{l}_{11}{K}_n\left({\alpha }_{\mathrm{P}1}a\right)\\ -k_{\mathrm{P}2}^2{C}_1{l}_{12}{K}_n\left({\alpha }_{\mathrm{P}2}a\right)-k_{\mathrm{P}3}^2{D}_1{l}_{13}{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)\\ +\left(C_{13}-\frac{1}{3}{\mu }_{13}\right)(-k_{\mathrm{P}1}^2{B}_1\frac{l_{11}}{l_{21}}{K}_n\left({\alpha }_{\mathrm{P}1}a\right)-k_{\mathrm{P}2}^2{C}_1\frac{l_{12}}{l_{22}}{K}_n\left({\alpha }_{\mathrm{P}2}a\right)\\ -k_{\mathrm{P}3}^2{D}_1\frac{l_{13}}{l_{23}}{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)+2{\mu }_1\frac{\mathsf{\partial }{\mu }_{r1}}{\mathsf{\partial }r}+{\mu }_{13}\frac{\mathsf{\partial }{\mu }_{r3}}{\mathsf{\partial }r},\end{array}$$

(19)

$$\begin{array}{l}{\sigma }_f=C_{12}(-k_{\mathrm{P}1}^2{B}_1{K}_n({\alpha }_{\mathrm{P}1}a)-k_{\mathrm{P}2}^2{C}_1{K}_n({\alpha }_{\mathrm{P}2}a)\\ -k_{\mathrm{P}3}^2{D}_1{K}_n({\alpha }_{\mathrm{P}3}a))\cos n(\theta -{\theta }_0)+K_2(-k_{\mathrm{P}1}^2{B}_1{l}_{11}{K}_n({\alpha }_{\mathrm{P}1}a)\\ -k_{\mathrm{P}2}^2{C}_1{l}_{12}{K}_n({\alpha }_{\mathrm{P}2}a)-k_{\mathrm{P}3}^2{D}_1{l}_{13}{K}_n({\alpha }_{\mathrm{P}3}a))\cos n(\theta -{\theta }_0)\\ +(C_{23}(-k_{\mathrm{P}1}^2{B}_1\frac{l_{11}}{l_{21}}{K}_n({\alpha }_{\mathrm{P}1}a)-k_{\mathrm{P}2}^2{C}_1\frac{l_{12}}{l_{22}}{K}_n({\alpha }_{\mathrm{P}2}a)\\ -k_{\mathrm{P}3}^2{D}_1\frac{l_{13}}{l_{23}}{K}_n({\alpha }_{\mathrm{P}3}a))\cos n(\theta -{\theta }_0),\end{array}$$

(20)

$$\begin{array}{l}{\sigma }_{rrh}=\left(C_{13}-\frac{1}{3}{\mu }_{13}\right)(-k_{\mathrm{P}1}^2{B}_1{K}_n\left({\alpha }_{\mathrm{P}1}a\right)-k_{\mathrm{P}2}^2{C}_1{K}_n\left({\alpha }_{\mathrm{P}2}a\right)\\ -k_{\mathrm{P}3}^2{D}_1{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)+C_{23}(-k_{\mathrm{P}1}^2{B}_1{l}_{11}{K}_n\left({\alpha }_{\mathrm{P}1}a\right)\\ -k_{\mathrm{P}2}^2{C}_1{l}_{12}{K}_n\left({\alpha }_{\mathrm{P}2}a\right)-k_{\mathrm{P}3}^2{D}_1{l}_{13}{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)\\ +\left(K_3-\frac{2}{3}{\mu }_{13}\right)(-k_{\mathrm{P}1}^2{B}_1\frac{l_{11}}{l_{21}}{K}_n\left({\alpha }_{\mathrm{P}1}a\right)-k_{\mathrm{P}2}^2{C}_1\frac{l_{12}}{l_{22}}{K}_n\left({\alpha }_{\mathrm{P}2}a\right)\\ -k_{\mathrm{P}3}^2{D}_1\frac{l_{13}}{l_{23}}{K}_n\left({\alpha }_{\mathrm{P}3}a\right))\cos n\left(\theta -{\theta }_0\right)+2{\mu }_3\frac{\mathsf{\partial }{\mu }_{r1}}{\mathsf{\partial }r}+{\mu }_{13}\frac{\mathsf{\partial }{\mu }_{r3}}{\mathsf{\partial }r},\end{array}$$

(21)

$${\sigma }_{rzs}={\mu }_1(\frac{\mathsf{\partial }{u}_{rs}}{\mathsf{\partial }z}+\frac{\mathsf{\partial }{u}_{zs}}{\mathsf{\partial }r})+\frac{1}{2}{\mu }_{13}(\frac{\mathsf{\partial }{u}_{rh}}{\mathsf{\partial }z}+\frac{\mathsf{\partial }{u}_{zh}}{\mathsf{\partial }r}),$$

(22)

$${\sigma }_{rzh}={\mu }_3(\frac{\mathsf{\partial }{u}_{rh}}{\mathsf{\partial }z}+\frac{\mathsf{\partial }{u}_{zs}}{\mathsf{\partial }r})+\frac{1}{2}{\mu }_{13}(\frac{\mathsf{\partial }{u}_{rs}}{\mathsf{\partial }z}+\frac{\mathsf{\partial }{u}_{zs}}{\mathsf{\partial }r}),$$

(23)

$${\sigma }_{r\theta s}={\mu }_1(\frac{\mathsf{\partial }{u}_{\theta s}}{\mathsf{\partial }r}-\frac{u_{\theta s}}{r}+\frac{1}{r}\frac{\mathsf{\partial }{u}_{rs}}{\mathsf{\partial }\theta })+\frac{1}{2}{\mu }_{13}(\frac{\mathsf{\partial }{u}_{\theta h}}{\mathsf{\partial }r}-\frac{u_{\theta h}}{r}+\frac{1}{r}\frac{\mathsf{\partial }{u}_{rh}}{\mathsf{\partial }\theta }),$$

(24)

$${\sigma }_{r\theta h}={\mu }_3(\frac{\mathsf{\partial }{u}_{\theta h}}{\mathsf{\partial }r}-\frac{u_{\theta h}}{r}+\frac{1}{r}\frac{\mathsf{\partial }{u}_{rh}}{\mathsf{\partial }\theta })+\frac{1}{2}{\mu }_{13}(\frac{\mathsf{\partial }{u}_{\theta s}}{\mathsf{\partial }r}-\frac{u_{\theta s}}{r}+\frac{1}{r}\frac{\mathsf{\partial }{u}_{rs}}{\mathsf{\partial }\theta }),$$

(25)

where $K_1$, $K_2$, and $K_3$ are the bulk moduli related to the solid grain, pore fluid, and gas hydrate; $C_{12}$, denotes the grain–fluid elastic coupling coefficient; $C_{23}$ refers to the fluid–hydrate elastic coupling coefficient; and $C_{13}$ represents the grain–hydrate elastic coupling coefficient.

The count of boundary conditions aligns with the number of undetermined coefficients both within and outside the well. By amalgamating the displacement derived from compressional and shear wave potentials with the boundary conditions specific to each phase, the ensuing linear equations can be formulated:

$$\left(\begin{array}{cccccccc}{m}_{11}& m_{12}& m_{13}& m_{14}& m_{15}& m_{16}& m_{17}& m_{18}\\ m_{21}& m_{22}& m_{23}& m_{24}& m_{25}& m_{26}& m_{27}& m_{28}\\ m_{31}& m_{32}& m_{33}& m_{34}& m_{35}& m_{36}& m_{37}& m_{38}\\ m_{41}& m_{42}& m_{43}& m_{44}& m_{45}& m_{46}& m_{47}& m_{48}\\ m_{51}& m_{52}& m_{53}& m_{54}& m_{55}& m_{56}& m_{57}& m_{58}\\ m_{61}& m_{62}& m_{63}& m_{64}& m_{65}& m_{66}& m_{67}& m_{68}\\ m_{71}& m_{72}& m_{73}& m_{74}& m_{75}& m_{76}& m_{77}& m_{78}\\ m_{81}& m_{82}& m_{83}& m_{84}& m_{85}& m_{86}& m_{87}& m_{88}\end{array}\right)\left(\begin{array}{c}{A}_f\\ B_m\\ C_m\\ D_m\\ E_m\\ F_m\\ G_m\\ H_m\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\\ b_4\\ 0\\ 0\\ 0\\ 0\end{array}\right),$$

(26)

By equating the determinant of the ${\left[m_{ij}\right]}_{8\times 8}$ to $D(k_z,\omega )$, and

$$N\left(k_z,\omega \right)=\left(\begin{array}{cccccccc}{b}_1& m_{12}& m_{13}& m_{14}& m_{15}& m_{16}& m_{17}& m_{18}\\ b_2& m_{22}& m_{23}& m_{24}& m_{25}& m_{26}& m_{27}& m_{28}\\ b_3& m_{32}& m_{33}& m_{34}& m_{35}& m_{36}& m_{37}& m_{38}\\ b_4& m_{42}& m_{43}& m_{44}& m_{45}& m_{46}& m_{47}& m_{48}\\ 0& m_{52}& m_{53}& m_{54}& m_{55}& m_{56}& m_{57}& m_{58}\\ 0& m_{62}& m_{63}& m_{64}& m_{65}& m_{66}& m_{67}& m_{68}\\ 0& m_{72}& m_{73}& m_{74}& m_{75}& m_{76}& m_{77}& m_{78}\\ 0& m_{82}& m_{83}& m_{84}& m_{85}& m_{86}& m_{87}& m_{88}\end{array}\right).$$

The borehole reflection coefficient $A_f$ can be written as

$$A_f(k_z,\omega )=\frac{N(k_z,\omega )}{D(k_z,\omega )}.$$

(27)

By substituting Equation (27) into Equation (1), we can deduce the real waveform in the time domain by employing the Fourier transform. The expressions for m_ij (i, j = 1, 2, 3, 4, 5, 6, 7, 8) and b_i (i = 1, 2, 3, 4) are given in Appendix A.

3. Influence of Gas Hydrate Saturation on the Flexural Wave

The characteristics of the borehole wavefield are significantly influenced by reservoir parameters. Parameters like permeability and borehole characteristics notably impact the dipole flexural wave. In contrast to fluid-saturated porous media, the saturation of the pore solid in three-phase porous media remains a primary research focus. This is crucial for inverting the gas hydrate saturation. Yet, there is a noticeable gap in the research concerning the acoustic attributes of three-phase porous media and the impact of hydrate saturation on the characteristics of flexural waves. In our prior work, we introduced an approach centered on component waves, which integrates the branch-cut integral with the residue theorem, factoring in residues at leaky poles [38]. Sensitivity serves as a metric to highlight how a parameter affects the phase velocity and attenuation of waves across varying frequencies. In this section, we calculate the sensitivity of the flexural wave’s phase velocity and attenuation to diverse parameters, aiming to discern the impact of these parameters on the flexural wave. The sensitivity is presented as [39]:

$$S_p(v)=\frac{p\mathsf{\partial }v}{v\mathsf{\partial }p}, S_p(Q^{-1})=\frac{p\mathsf{\partial }{Q}^{-1}}{Q^{-1}\mathsf{\partial }p},$$

(28)

where $S$ indicates the sensitivity, $p$ represents the medium parameter, $v$ signifies the velocity, and $Q^{-1}$ refers to the attenuation. A sensitivity value greater than zero suggests that the characteristic increases with an increment in the variable, while a value less than zero implies a decrease.

Figure 1 shows the sensitivity of the attenuation of the flexural wave to hydrate saturation. The base value of the solid frame permeability ${\kappa }_{s0}=1.07\times {10}^{-13}$ m² [14,15]. Porosities of 30–80% fundamentally cover the scope observed in hydrate-bearing sediments [28]; thus, the basic value of the porosity is chosen to be 0.5. The gas hydrate saturation is chosen to be 0.5, the pore fluid viscosity is $1.8\times {10}^{-3}$ Pa·s [14,15,16,30], and the borehole fluid velocity is 1500 m/s. The physical parameters [14,15,30] are presented in

Table 1. Physical parameters of different phases [14,15,30].

Parameter	Value
Grain density (kg/m3)	2650
Water density (kg/m3)	1000
Gas hydrate density (kg/m3)	900
Grain bulk modulus (GPa)	38.7
Grain shear modulus (GPa)	39.6
Water bulk modulus (GPa)	2.25
Gas hydrate bulk modulus (GPa)	8.58
Gas hydrate shear modulus (GPa)	3.32

, respectively. As can be observed from Figure 1, the flexural wave attenuation is sensitive to gas hydrate saturation, and the flexural wave attenuation decreases as the gas hydrate saturation increases.

4. The Inversion Method and Application to Field Data

Based on the sensitivity analysis of flexural wave attenuation, it can be seen that the gas hydrate saturation will have a certain influence on the flexural wave attenuation within a borehole of a gas hydrate-bearing sediment. Therefore, in this section, a method is proposed to invert the gas hydrate saturation by using flexural wave attenuation.

Assuming the amplitude of the mode wave is A₀ at a certain point, and the amplitude of the wave is A₁ after it has propagated L meters along the axis, then the relationship between the amplitude and the imaginary part of the pole wavenumber is given by

$$\frac{A_0}{L{A}_1}=\exp \left(\mathrm{Im}\left(k_{\mathrm{pole} }\right)\right).$$

(29)

As can be observed from the above equation, when the imaginary part is zero, A₀ = A₁, which means that the amplitude of the wave does not decay when propagating axially. When the imaginary part is not zero, A₀ > A₁, indicating that the amplitude of the wave decreases with axial propagation, and the larger the imaginary part value of the pole wavenumber, the faster the wave attenuation. In the field of acoustic instruments, decibel per meter (dB/m) is often used to define the attenuation, and the conversion relationship between the attenuation in dB/m and the imaginary part of the pole wavenumber is

$$\mathrm{Attenuation} \mathrm{in} \mathrm{dB}/\mathrm{m}=20\log e^{\mathrm{Im}\left(k_{\mathrm{pole} }\right)}=8.686\times \mathrm{Im}\left(k_{\mathrm{pole} }\right),$$

(30)

The inversion process can be achieved by minimizing the following objective function, constructed as the root-mean-square misfit error between the attenuation curve of the model and the measured attenuation data for the processing frequency range:

$$E(S_h)={\displaystyle \sum _{\mathsf{\Omega }}\left|\right|{\alpha }_m(S_h;\omega )-{\alpha }_d(\omega )|{|}^2},$$

(31)

where $\mathsf{\Omega }$ is the frequency range of the calculated objective function, ${\alpha }_m$ denotes the attenuation curve calculated using the theoretical model, and ${\alpha }_d$ refers to the processed attenuation data, given by

$${\alpha }_d(\omega )=\frac{1}{L}20\log [\frac{A_0(\omega )}{A_1(\omega )}],$$

(32)

and the Euclidean norm of the mismatch error can be used to simultaneously invert the gas hydrate saturation of the reservoir.

An example of a gas hydrate-bearing sediment proves that Equation (13) can be used to invert the gas hydrate saturation of the reservoir. Figure 2a,b show the borehole waveforms at distances of 2 m and 3 m, which are simulated by the staggered-grid finite-difference time-domain method in our previous work [34]. A rectangular window was used to intercept the flexural wave. Figure 3 represents the frequency spectrum of the flexural wave at a distance of 2 m and 3 m. The theoretical attenuation curve of ${\alpha }_m$ is obtained by changing the gas hydrate saturation of the reservoir; the parameters employed for calculating ${\alpha }_d$ are enumerated in Table 1. Figure 4 shows the inversion result calculated using the flexural wave attenuation; the calculation frequency band $\mathsf{\Omega }$ is chosen to be 5.35 kHz to 7.75 kHz. The gas hydrate saturation obtained by inversion is 0.492, with an error of 1.67%. Figure 5 shows the attenuation of the flexural wave of the theoretical calculation result and the processed attenuation data calculated by the simulated waveforms. As observed in Figure 5, in this frequency range, the theoretical calculation result is in good agreement with the processed attenuation data.

To verify the reliability of this inversion method for signals with a signal-to-noise ratio (SNR), we artificially add noise to the data of the numerical simulation. Here, the low-frequency noise is processed by Gaussian noise through a low-pass filter in the 0–6.5 kHz frequency band. Each receiver is added with linearly independent noise signals, and the signal-to-noise ratio (SNR) is five. The central frequency of the monopole sources is selected as 6 kHz, and the borehole waveforms after adding noise are shown in Figure 6. It can be seen that noise seriously interferes with the waveforms. The waveform data after adding noise are used to invert the gas hydrate saturation by the inversion method proposed in this paper; the frequency spectrum is shown in Figure 7. As shown in Figure 8, after adding noise, the inversion result of gas hydrate saturation is 0.495, and the error is 1%. Figure 9 shows a comparison between the attenuation of the theoretical calculation and the attenuation of the synthetic data.

To further verify the reliability of this inversion method for signals with a lower signal-to-noise ratio (SNR), we artificially add stronger noise to the data of numerical simulation. Here, the SNR is 3. The central frequency of the monopole sources is selected as 6 kHz, and the borehole waveforms after adding noise are shown in Figure 10. It can be seen that noise seriously interferes with the waveforms. The waveform data with added noise are used to invert the gas hydrate saturation by the inversion method proposed in this paper; the frequency spectrum is shown in Figure 11. As shown in Figure 12, after adding noise, the inversion result of gas hydrate saturation is 0.51, and the error is 2.33%. This result further proves the effectiveness of the inversion method. Figure 13 shows a comparison between the attenuation of the theoretical calculations and the attenuation of the processed data.

5. Conclusions

In this paper, based on the three-phase porous media theory, the borehole wavefield of a three-phase porous medium with multipole sources is calculated for the first time. The sensitivity of the flexural wave attenuation to gas hydrate saturation in a gas hydrate-bearing sediment is calculated and analyzed. A method using flexural wave attenuation to invert gas hydrate saturation is proposed. On this basis, the gas hydrate saturation of a gas hydrate-bearing sediment is estimated. The findings of our analysis are as follows.

(1) The attenuation of the flexural wave from dipole acoustic logging is sensitive to the gas hydrate saturation. Specifically, the attenuation of the flexural wave decreases with increasing hydrate saturation. The flexural wave attenuation can therefore be used to quantify the gas hydrate saturation.

(2) Based on the flexural wave analysis, we have developed an inversion method to simultaneously determine the gas hydrate saturation. We demonstrate that the attenuation of the flexural wave can effectively invert the gas hydrate saturation of natural gas hydrate-bearing sediments.

(3) In the inversion section, our research has some shortcomings. We have only theoretically demonstrated the feasibility of using dipole flexural waves to invert the saturation of natural gas hydrates, ignoring the complexity of the reservoir and the presence of logging tools that can have an impact on the results. In addition, real logging curves are often more complex than what we synthesized. However, we were unable to obtain real on-site data, which prevented further analysis. Nevertheless, our research provides a theoretical foundation for the inversion of hydrate saturation using the attenuation of the flexural wave.