Acoustic Wave Propagation in a Borehole with a Gas Hydrate-Bearing Sediment

<jats:p>A knowledge of wave propagation in boreholes with gas hydrate-bearing sediments, a typical three-phase porous medium, is of great significance for better applications of acoustic logging information on the exploitation of gas hydrate. To study the wave propagation in such waveguides based on the Carcione–Leclaire three-phase theory, according to the equations of motion and constitutive relations, a staggered-grid finite-difference time-domain (FDTD) scheme and a real axis integration (RAI) algorithm in a two-dimensional (2D) cylindrical coordinate system are proposed. In the FDTD scheme, the partition method is used to solve the stiff problem, and the nonsplitting perfect matched layer (NPML) scheme is extended to solve the problem of the false reflection waves from the artificial boundaries of the computational region. In the RAI algorithm, combined with six boundary conditions, the displacement potentials of waves are studied to calculate the borehole acoustic wavefields. The effectiveness is verified by comparing the results of the two algorithms. On this basis, the acoustic logs within a gas hydrate-bearing sediment are investigated. In particular, the wave field in a borehole is analyzed and the amplitude of a Stoneley wave under different hydrate saturations is studied. The results indicate that the attenuation coefficient of the Stoneley wave increases with the increase of gas hydrate saturation. The acoustic responses in a borehole embedded in a horizontally stratified hydrate formation are also simulated by using the proposed FDTD scheme. The result shows that the amplitude of the Stoneley wave from the upper interface is smaller than that from the bottom interface.</jats:p>

section_content: Natural gas hydrate is globally widespread in continental margin sediments and permafrost regions, it has been recognized as a potential energy source in the 21st century [1, 2]. The exploration and development technology of natural gas hydrate has been a research hotspot in recent years [3] [4] [5] [6] [7]. Acoustic logging is one of the promising technologies in the interpretation of formations parameters underground [5, 8], where the velocity from acoustic well logs can be used to distinguish gas hydrate layers along the depth of the borehole. Recent research shows that the attenuation of the first arrival can be used to estimate gas hydrate saturation [9], the precision of the result, however, cannot be guaranteed [8]. To solve the problem and enrich the applications of borehole acoustic logging data in gas-hydrate production, investigations on wave propagation in the fluid-filled borehole surrounded by a gas hydrate formation are essential. 

In fact, natural gas hydrate-bearing sediments are typical three-phase porous media (solid grain, gas hydrate, and water), which implies a particular topological configuration, namely the one where solid grain and gas hydrate form two continuous and interpenetrating solid matrices. Consequently, studies of wave propagation in this complex environment become increasingly important to obtain more detailed reservoir information. Leclaire proposed the percolation theory based on Biot's theory and analyzed the acoustic wave propagation in frozen porous media such as frozen soil or permafrost [10]. The theory predicts three compressional and two shear waves are produced when an acoustic wave propagates in frozen porous media. However, Leclaire assumed that there was no direct contact between solid grains and ice. On this basis, Carcione et al. put forward a Biot-type three-phase theory of wave propagation in frozen porous media (three-phase porous media) that allows for the interaction between the solid grain frame and pore solid, which is named the Carcione-Leclaire three-phase theory [11]. The Carcione-Leclaire theory predicts three compressional and two shear waves in a three-phase porous medium: the first kind of compressional (P1), shear (S1) waves, the second kind of compressional (P2), shear (S2) waves, and the third kind of compressional (P3) wave. Up to date, much work has been done to study the characteristics of gas hydrate-bearing sediments [12] [13] [14] [15] [16] [17] [18]. However, there is almost no research on the wave propagation in a borehole surrounded by a hydrate formation, therefore, how to better use the acoustic well logging information to evaluate the parameters of hydrate formation lacks sufficient theoretical guidance. 

RAI and FDTD are two effective methods to simulate the acoustic logs of porous media. Rosenbaum [19] first applied Biot's theory [20] [21] [22] [23] [24] [25] of porous media to the study of acoustic logging in a fluid-saturated porous formation and used the RAI method to simulate the full waveform of the borehole excited by a monopole source, which is the beginning of the theoretical study on acoustic logging in porous media. On this basis, several RAI algorithms have been applied to simulate the acoustic logging in a fluid-saturated porous formation, which is a valid method to deal with the borehole acoustic wavefield of a homogeneous medium [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]. For non-axisymmetric boreholes or inhomogeneous porous media, the RAI has limitations in simulating wave propagating problems. As a complement, the FDTD method can deal with the problem of heterogeneity [36]. Guan and Hu [37, 38] proposed a parameter averaging technique to deal with the interface of a discontinuous medium and applied it to the numerical simulation of the wave propagation in fluid-saturated porous media. They used the NPML technique [39] to eliminate the false reflected waves from the model boundaries in numerical simulations of wave propagation. Yan [40] et. al numerically simulate the acoustic fields excited by a point source in the borehole surrounded by a porous formation using the 3D staggered grid stress-velocity finite difference method. Ou [41] used the FDTD algorithm to simulate the wave propagation in a borehole within a porous medium and studied the characteristics of the Stoneley wave reflection. These studies show that much research has been performed on the acoustic logs with a two-phase porous medium. However, few attempts have studied the problems related to borehole wave propagations with a three-phase porous medium. 

In this paper, to model the acoustic responses excited by an axisymmetric point source on the borehole axis of a gas hydrate-bearing sediment (three-phase porous medium), we propose both a RAI algorithm and a 2D FDTD scheme. In the RAI algorithm based on the Carcione-Leclaire three-phase theory [11], we construct the displacement potentials of three compressional waves and two shear waves of three phases (solid grain frame, gas hydrate, and pore fluid) in the frequency wavenumber domain. In addition, the boundary conditions of the borehole wall of gas hydrate-bearing sediments are given. On these bases, combined with the equations of motion and constitutive relations, the borehole acoustic field of this multiphase porous medium is obtained. In the FDTD scheme, we use a partition method [11, 42, 43] to solve the stiff problem. Moreover, to truncate the computational region, the NPML scheme is extended to solve the problem of the false reflection waves from the artificial boundaries of the computational region. In addition, we use the proposed FDTD algorithm to describe the spatial distribution of the acoustic field and the waveforms of the monopole acoustic responses in a gas hydrate-bearing sediment. Further, to explore the influence of gas hydrate saturation on Stoneley wave, we simulate waveforms in the borehole with different gas hydrate saturations and calculate the attenuation coefficients. Finally, we simulate the acoustic responses in a horizontally stratified gas hydrate-bearing sediment with the proposed FDTD algorithm. The main contribution of this paper is to construct and simulate the borehole acoustic fields of natural gas hydrate-bearing sediments, which expand the knowledge of borehole acoustics in three-phase media. The research findings can provide a theoretical basis and guidance in applying acoustic logging information for the exploration and evaluation of natural gas hydrate.
section_name: Introduction

section_content: According to the Carcione-Leclaire three-phase theory [11], combined with the constitutive relations and equations of motion, the first-order velocity-stress equations in a gas hydrate-bearing sediment can be obtained:

where σ ij is a stress component of the solids; σ represents the stress of the fluid (here is water); the subscripts and superscripts 1, 2, and 3 denote the solid grains; pore fluid and gas hydrate, respectively; K 1, K 2, and K 3 refer to the bulk moduli; µ 1, µ 3 and µ 13 are the shear moduli; C 12 denotes the grain-fluid elastic coupling coefficient; C 23 is the fluid-hydrate elastic coupling coefficient; C 13 represents the grain-hydrate elastic coupling coefficient; and θ m, d ij define the strain tensors:

The equations of motion [11] can be expressed as

i

,j = ρ 12.
section_name: Carcione-Leclaire Three-Phase Theory

section_content: (1)

i

where ρ 11, ρ 22 and ρ 33 represent the mass densities of the effective solid grain, effective fluid, and effective hydrate, respectively; ρ 12, ρ 13 and ρ 23 refer to the coupling mass densities; is the temporal derivative of the particle velocity. Equations ( 4 )-( 6 ) can be rewritten as

where

i, v

,ib 12 v

(2)

In this manner, the following equations of motion are obtained, which are conducive to the establishing of the numerical simulation method:. v (1) i

(2) i

where γ ij and Π i are the elements in the matrices γ and Π. Equations ( 1 )-( 3 ) and ( 8 )- (10) constitute the first-order velocity-stress equations for wave propagation in gas hydrate bearing-sediments based on the Carcione-Leclaire theory. 

Carcione and Quiroga-Goode discussed the eigenvalues of the propagation matrix of Biot's acoustic equations in a low-frequency range [42]. Zhao applied the partition method to the staggered high-order finite-difference method [43]. Here, the same method is employed in our work because of the existence of the friction coefficients b 12, b 23 and b 13 [11, 44]. The stiff part can be solved analytically, then the non-stiff part of the velocitystress equations is treated by the staggered-grid finite-difference method. In this manner, the problem of the small time step caused by the existence of a friction coefficient is avoided. 

The stiff part of the velocity-stress differential equations can be written in the matrix form as

where In the following discussion, the treatment of a porous medium found in Ref. [44] is used for our case. The analytical solution of Equation ( 11 ) in a 2D situation can be expressed as

where

. ξ 1 and ξ 2 are two eigenvalues of exp(S∆t). 

A matrix form is used for the non-stiff part of the velocity-stress differential equations, i. e.,

where B = (σ

,i, σ

3. Finite Difference Scheme 3. 1. Discretization of Equations in the Carcione-Leclaire Three-Phase Theory Figure 1 presents a schematic diagram of the acoustic logging model. A fluid-filled borehole is embedded in an unbounded three-phase porous medium (gas hydrate-bearing sediment). With a cylindrical coordinate system (r, z, θ), the borehole axis lies along the central axis z of the cylinder, the monopole acoustic source is located at the center of the borehole which is placed at the origin of the coordinate system, the squares represent the receivers. As the medium is axisymmetric, all field quantities are independent of the θ-coordinate. The borehole radius a is 0. 1 m, ρ f denotes the fluid (water) density, and v f represents the fluid (water) velocity. First, Equations ( 1 )-(3) in a 2D axisymmetric cylindrical coordinate system can be rewritten as

∂σ (2) ∂t = C 12 ( ∂v (1)

)

where σ rr, σ θθ, σ zz and σ rz are the components of the bulk stress tensor. 

We define the following intermediate vector as the input for the staggered-grid finitedifference algorithm:

rz, σ (2), σ

rz ). The non-stiff part in a 2D axisymmetric cylindrical coordinate system can be written as

where v

(1)

r and v

(2)

r and v

z refer to the components of the particlevelocity vector of solid grain frame, fluid and gas hydrate, respectively. 

To implement a 2D finite-difference algorithm for the solutions of Equations ( 14 )-( 16 ) and ( 18 ), all quantities of the stress and velocity components should be discretized in staggered grids. We discretize the velocities and stresses using the staggered-grid finitedifference algorithm. The staggered-grid finite-difference form of the non-stiff part in a 2D axisymmetric cylindrical coordinate system can be expressed as (Here we only give the staggered-grid finite-difference schemes for the solid grain frame including particle velocity and stress components, which is similar to those of pore fluid and hydrate)

where n denotes the time; the subscripts i, j indicate the spatial position r, z of the field quantities, respectively; ∆t is the time interval; D r and D z are discretized differential operators in r and z directions, respectively. In the 2D case, the distribution of the field quantities in the staggered-grid finite-difference method is shown in Figure 2. NPML is used to eliminate the false reflection generated by artificial boundaries. In this work, the solid grain frame is taken as an example to analyze the FDTD algorithm for the first-order velocity-stress equations after using the NPML scheme. For the particle velocity components, we have

where

P n zzz = e -Ω z ∆t P n-1 zzz -

P n rzr = e -Ω r ∆t P n-1 rzr -

where L is the width of the PML, α = 10 -6 is a predefined level of wave absorption, and the coefficients a = 0. 25 and b = 0. 75 are used. Similar to the processing of Equations ( 25 ) and (26), Equations ( 21 )-( 24 ) can be discretized as

(3) (n+1/2) r(i,j)

(3) (n+1/2) r(i,j)

(3) (n+1/2) z(i,j)

(1) (n+1/2) z(i,j)

(2) (n+1/2) r(i,j)

(n+1/2) z(i,j)

(n+1/2) r(i,j)

(1) (n+1)
section_name: v

section_content: (3) (n+1/2) rz

where
section_name: +P

section_content: To unify the equations of acoustic wave propagation in both the outer three-phase porous medium and the inner fluid in the borehole, the limits of certain field quantities and medium parameters in the three-phase porous medium are set to obtain the equations for the fluid. In this limiting case, the stresses of the solid grain frame phase and the gas hydrate phase in the first-order velocity-stress equations of the three-phase porous medium are zero. The parameters in the Carcione-Leclaire three-phase theory change to

In addition, according to the wave equations in fluid, it can be directly defined

Therefore, the equations of a three-phase porous medium degenerate into governing equations in the fluid.
section_name: Discretization of the Equations in the Borehole and Borehole Wall

section_content: In the frequency wavenumber domain, the acoustic field in the borehole can be expressed as

where F denotes the frequency of the source; A is an undetermined coefficient, which can be regarded as the reflection coefficient of the wellbore; K 0 and I 0 are Bessel functions; r refers to the radial wave number of the fluid. The acoustic field in a three-phase porous medium can be expressed by the displacements of solid grain frame u 1, pore fluid u 2 and gas hydrate u 3. The three displacements produced by a point source can be written as

where ϕ m and η m represent the displacements of compressional waves (P1, P2 and P3) and shear waves (S1 and S2), respectively. The subscripts m = 1, 2, 3 refer to the solid grain frame, pore fluid and hydrate, respectively. The displacement potentials of the compressional waves in the frequency wave number domain can be expressed as

where B m, C m and D m are the undetermined amplitudes of the three compressional waves, respectively; α P1, α P2 and α P3 are the radial wave numbers of the three compressional waves, respectively. Moreover, in the frequency wavenumber domain, the displacement potentials of shear waves are given as

where E m and F m are the undetermined amplitudes of the two shear waves, respectively; α S1 and α S2 are the radial wave numbers of the two shear waves, respectively. In the above formulas, A, B m, C m, D m, E m and F m are undetermined coefficients, which are determined by the boundary conditions of the borehole wall. Assuming that there is a relationship between the compressional wave displacement of each phase:

then we can obtain

Assuming that there is a relationship between the shear wave displacement of each phase:

then we can get

The boundary conditions studied in this paper are the open hole boundary conditions, which can be written as

where the subscripts in and out refer to inside and outside the borehole wall, respectively, and the subscripts s, f, and h denote solid grains, pore fluid, and hydrate, respectively. The first equation represents the continuity of fluid flow in the normal direction of the borehole wall. The second and the third equations show that the pressure of fluid in the well in the normal direction of the borehole wall is equal to the normal stress of the unit body in the gas hydrate formation outside the well in the normal direction of the borehole wall. The fourth equation represents that the pressure of the fluid in the well in the normal direction of the interface should be equal to that of the unit body in the porous formation outside the well. The fifth and sixth equations represent the continuity of shear stresses in the tangential direction of the borehole wall. w p and w s are the flowing displacement of the compressional and shear wave, which can be expressed as

where φ s, φ f and φ h refer to the volume fraction of solid grains, pore fluid, and gas hydrate,

The specific solving process of each parameter in matrix M is shown in Appendix A. The reflection coefficient A in the borehole is given by

By substituting Equation (58) into Equation ( 46 ), the actual waveform received in the time domain can be obtained by the double Fourier transform.
section_name: RAI Algorithm

section_name: Numerical Modeling

section_content: Figure 3 gives the simulated waveforms of the FDTD and the RAI method at different locations from z = 0. 5 m to z = 4 m with different source dominant frequencies 6 kHz and 10 kHz in the gas hydrate formation. The physical parameters of the gas hydrate-bearing sediment are given in Table 1. The acoustic source pulse width T c = 0. 5 ms is employed in our simulation. The radius of the borehole is 0. 1 m. In our FDTD scheme, the time step is 2 × 10 -4 ms and the space steps in the r and z directions are both chosen as 5 × 10 -4 m. The radial length of the model is 0. 5 m, while the axial length is 4. 0 m. The width L of the NPML region is 0. 5 m. The source pulse function s(t) used in this paper is expressed as

where f 0 refers to the center frequency of the source, and t denotes the time. Figure 3b shows that there are three different wave groups of the waveforms: the first kind of compressional (P1) wave (A), the first kind of shear (S1) wave (B), and the Stoneley wave (C). The amplitude of the P1 wave is small, while the S1 wave and Stoneley wave are more evident. As depicted in Figure 3, the FDTD algorithm simulations (black solid lines) agree well with the RAI solutions (red dash-dot lines) at different locations from z = 0. 5 m to z = 4 m. Because of the discrete Fourier transform of the RAI method and the discrete FDTD grids, there are small differences in the amplitudes of the waveforms of the two methods [41]. 

As an example of using the finite-difference scheme proposed above, we simulate the acoustic field in a wellbore surrounded by a three-phase porous medium. The wave field distribution at instants of 0. 3, 0. 5, 0. 8, 1. 1, and 1. 4 ms is displayed in Figure 4. The black dotted line represents the borehole wall. To identify different wave groups in the wave field, a domain frequency f 0 = 20 kHz and a shorter pulse width T c = 0. 2 ms are employed in our simulation. The time and space steps are consistent with those in Figure 3. Figure 4a -e shows that there is no evident wave reflection from the outer boundaries, which indicates that the wave field is well absorbed in the NPML. In Figure 4a, we can observe the fastest pulse in the formation is the P1 wave and then the S1 wave. When the direct wave propagates to the borehole wall, part of the energy propagates outward, forming the refraction waves. Because there is no friction between the solid grain and the gas hydrate and the viscosity of the pore fluid is 1. 8 × 10 -3 Pa•s, four waves exist in the gas hydrate formation: two compressional waves (P1 and P2) and two shear waves (S1 and S2) [44]. As indicated in Figure 4, the energy of the Stoneley wave is mainly concentrated near the interface of the borehole wall, which is the same as that in a two-phase porous medium [37]. With the increase of time, the amplitude of the Stoneley wave attenuates. Furthermore, the compressional head wave in the borehole can also be observed.
section_name: Homogeneous Three-Phase Porous Media

section_content: Grain density (kg/m 3 ) 2650 Fluid density (kg/m 3 ) 1000 Gas hydrate density (kg/m 3 ) 900 Grain bulk modulus (GPa) 38. 7 Fluid bulk modulus (GPa) 2. 25 Gas hydrate bulk modulus (GPa) 8. 85 Gas hydrate shear modulus (GPa) 3. 32 Grain frame permeability κ so (m 2 )

1. 07 × 10 -13 Gas hydrate frame permeability κ ho (m 2 ) 5 × 10 -4

Figure 5 gives the simulated waveforms with a dominant frequency f 0 = 20 kHz and a pulse width T c = 0. 2 ms, the frequency range is 10 kHz-30 kHz. Figure 5 shows that there exist four different wave groups of the waveforms: the first kind of compressional (P1) wave, the first kind of shear (S1) wave, the pseudo-Rayleigh wave, and the Stoneley wave. The amplitudes of the pseudo-Rayleigh wave and Stoneley wave are much larger than those of the P1 and S1 waves. Because of the different physical characteristics of each component, there is a certain difference among the propagation velocities of the four waves. The pseudo-Rayleigh wave is a kind of interface wave, and its low-frequency cutoff velocity is usually the velocity of the S1 wave, and the velocity difference between the pseudo-Rayleigh wave and S1 wave is small. As can be observed in Figure 5, with the increase of the distance between the receiver and the source, the amplitudes of the four modes become smaller. The relative motion between the viscous fluid and the solid phase of the three-phase porous media causes viscous dissipation and dissipation of wave energy, the four waves are attenuated in the process of propagation. Moreover, P2 and S2 waves are not observed in the wave field because of their small amplitudes and large attenuation along the axis of the well. The estimation of gas hydrate saturation in a hydrate reservoir is a key index in hydrate exploration. The P1-wave velocity and attenuation were used to evaluate hydrate saturation in previous studies [45, 46]. However, there is no research on the analysis of the Stoneley wave in gas hydrate reservoirs. The Stoneley wave is usually used to evaluate the permeability of fluid-saturated porous media. Here, we further investigate the waveforms in the borehole with different gas hydrate saturations (Figure 6 ) with the domain frequency f 0 = 10 kHz. Gas hydrate saturations of 0. 4, 0. 5, and 0. 6 are considered. From Figure 6, it can be observed that the variation of the Stoneley wave amplitude with gas hydrate saturation is larger than those of the P1 wave and the S1 wave. Further, with the increase of hydrate saturation, the velocity of the Stoneley wave is almost unchanged, while the amplitude of the Stoneley wave decreases significantly. Under the condition of the same gas hydrate saturation, with the increase of the source-receiver distance, the Stoneley wave attenuates, and its amplitude decreases. Further, the attenuation coefficient of the Stoneley wave is calculated at different gas hydrate saturations. The amplitude of Stoneley wave A [47] can be expressed as

where M is sampling numbers; W i denotes the amplitude of the sampling point i. The attenuation coefficient of the Stoneley wave [47] is given by

where A n and A m represent the amplitudes of the received waveforms of the n-th and the m-th receivers, respectively; d s is the receiver spacing. After calculation, we find that when values of the hydrate saturation are 0. 4, 0. 5, and 0. 6, the attenuation coefficients of the Stoneley wave are 0. 35, 0. 52, and 0. 87, respectively. The results indicate that the attenuation coefficient of the Stoneley wave increases with the increase of gas hydrate saturation. Thus, we can use the attenuation of the Stoneley wave to estimate the gas hydrate saturation of the reservoir.
section_name: Parameter Value

section_content: Presented in Figure 7 is a schematic diagram of a horizontally stratified medium across a borehole on the r-z plan. The borehole is filled with water. As shown in Figure 7, there are three layers in the formation around the borehole, the layer between 1. 0 m and 2. 5 m is the gas hydrate-bearing sediment, the other two layers are elastic media. The parameters of the gas hydrate formation are consistent with those shown in Table 1, while the parameters in the other two elastic layers are the same as those of the solid grains in the hydrate formation. Fluid-hydrate, solid-hydrate, and fluid-solid interfaces exist in this model. We use the FDTD algorithm combined with the parameter averaging technique to simulate the wave propagation in the borehole. Shown in Figure 8 is the wave-field distribution at instants of 0. 5, 1. 0, 1. 5, and 2. 0 ms. The interfaces are marked with two blue lines. In our simulation, the source has a center frequency and a pulse width of f 0 = 10 kHz and a pulse width T c = 0. 5 ms, the frequency range is 6 kHz-14 kHz. The time and space steps are consistent with those in the homogeneous three-phase medium. In the figure, no evident wave reflection is observed from the outer boundaries, indicating that the wave-field is well absorbed in the NPMLs. As indicated in Figure 8, the energy of the Stoneley wave is reflected at the interface. Figure 9 displays the waveforms at different locations along the borehole. In Figure 9, the two blue lines indicate the reflected Stoneley waves at the two interfaces. It can be seen that the amplitude of the Stoneley wave from the upper interface is smaller than that from the bottom interface. This is because most of the wave energy is attenuated in the interlayer and transmitted from the interface. Further, because the layers of 0 m-1. 0 m and 2. 5 m-4. 0 m are two identical elastic media while the interlayer of 1. 0 m-2. 5 m is a gas hydrate-bearing sediment, it can be seen from the waveforms that the velocity of the Stoneley wave is the same when the distance between the receiver and the source is within 1 m and between 2. 5 m and 4 m. Besides, when the distance between the receiver and the source is between 1. 0 m and 2. 5 m, the velocity of the Stoneley wave is obviously smaller than the other two cases.
section_name: Horizontally Stratified Three-Phase Porous Medium

section_content: To extract hydrate information from acoustic data, we have developed both an RAI and a velocity-stress FDTD algorithm with a nonsplitting perfect matched layer to simulate acoustic wave propagation in the borehole surrounded by a gas hydrate-bearing sediment. The proposed FDTD algorithm can also be applied for simulations of acoustic logs in a heterogeneous three-phase porous medium. Clearly reflected Stoneley waves generated from the upper and lower interfaces of the different media can be observed in the waveforms. The above work can enrich our understanding of wave propagation in the borehole surrounded by a gas hydrate formation and provide a theoretical foundation for the exploration and development of hydrate. 

a 23 tortuosity for fluid flowing through the gas hydrate

a 13 tortuosity for solid grain flowing through the gas hydrate

a 31 tortuosity for gas hydrate flowing through the solid grain

Equations ( 50 )-(53) show the relationship between the unknown coefficients in displacements of compressional and shear waves of each phase. Therefore, in the frequency wavenumber domain, the displacement potential of compressional and shear waves can be expressed as )

. 

(A18)

Finally, we can get the values of each parameter in the M matrix, which are given by

Values in gradient and curl in Equations (A39)-(A41) are zero. Thus, there are six equations, which can be divided into two coupled equations. The first group is the gradient of the three equations, which is the compressional wave potential, and the value of this part is zero. Besides, ∇ 2 = -k. Take the real part of three wavenumbers as positive, and the real part of k P1 is less than that of k P2 which is less than that of k P3. The corresponding three k 2 P can get three l 1 and three l 2, which is given by l 11, l 12, l 13 and l 21, l 22 and l 23, corresponding to three kinds of compressional waves. 

The second group is the curl of the three equations, which is the shear wave potential, and the value of this part is zero. Besides, ∇ 2 = -k 2 S, the second group can be written in the matrix form as

Let the determinant of the left matrix be zero, we can obtain two values of k 2 S, substitute k 2 S into Equation (A43), two l 3 and two l 4 can be obtained, which can be expressed as l 31, l 32, l 41 and l 42, corresponding to two kinds of shear waves.
section_name: Conclusions

section_content: Funding: This research was funded by the National Natural Science Foundation of China (Grant Nos. 11974018, 11734017 ) and the Strategic Pilot and Technology Special of Chinese Academy of Sciences, China (Grant No. XDA14020303 ), and the Research Instrument and Equipment Development Project of Chinese Academy of Sciences-Key Technology Team Project (Grant No. GJJSTD20210008 ). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

section_content: Data Availability Statement: Not applicable.

section_name: Conflicts of Interest:

section_content: Author Contributions: The original idea was contributed by L. L. ; The manuscript preparation and simulation parts were finished by L. L. ; The follow-up guidance was given by X. Z. ; The supervisors of L. L. are X. W. and X. Z. ; The manuscript was improved and perfected by X. Z. and X. W. ; Y. J. participated in the discussion and gave meaningful suggestions. All authors have read and agreed to the published version of the manuscript.

section_content: The authors declare no conflict of interest. 

Substitute the displacement potential of compressional and shear waves into the boundary condition, respectively. For the first boundary condition at the borehole wall

where

For the second boundary condition,

The fluid pressure in the borehole P f is given by

The radial stress of solid grain frame in natural gas hydrate formation σ rrsout can be expressed as

rr + e

rr + e

θθ + e

zz )] + µ 13 [e

rr + e

θθ + e

zz )]

θθ + e

zz )µ 13 (e

θθ + e

zz ),

where e rr, e θθ and e zz are the strain tensors in three directions, which can be expressed as

where u z is the displacement in the z-direction, which is defined as

For the boundary condition

the radial stress of gas hydrate is given by

rr + e

θθ + e

(3)

rr + e

θθ + e

(1)

rr + e

(2)

rr + e

θθ + e

(3)

rr + e

θθ + e

zz )]

θθ + e

(3)

θθ + e

zz ). 

(A12)

For the boundary condition

the fluid pressure of the gas hydrate reservoir outside the borehole is

rr + e

θθ + e

(1)

rr + e

θθ + e

(2)

rr + e

θθ + e

(3)

For the boundary condition

The shear stresses of the solid grain frame and hydrate frame outside the well can be expressed as

In addition,

Among the unknowns of the above parameters, l 11, l 12, l 13 and l 21, l 22, l 23 correspond to three kinds of compressional waves, respectively. l 31, l 32 and l 41, l 42 correspond to two kinds of shear waves, respectively. They are constants determined by the medium parameters. 

The total stress component of the natural gas hydrate reservoir is given by

ij -

ij -

(A26)

The equation of motion can be expressed as

i + (ρ 13 + ρ 23 + ρ 33 )

. . 

where ρ ij refers to the coupling mass density. Substitute the equation of motion into the constitutive relation, we can obtain

i. 

(A28)

Similarly, the constitutive relations of pore fluid and hydrate phase are defined as

(A30)

The equations of motion of pore fluid and hydrate frame are

i + ρ 23. . 

i ),

i + ρ 23. . 

i ). 

(A32)

Substitute the equations of motion into the constitutive relations, we can obtain

i + ρ 23. . 

i -

i ),

(A33)

i,jj

i + ρ 23. . 

(A34)

Equations (A28), (A33) and (A34) form wave differential equations, the vector form in the frequency domain is given by

Any vector field can be decomposed into a sum of a P-wave field and two S-wave fields. If the z-axis is in a certain direction in space, we have Substituting Equation (A38) into (A35), (A36) and (A37), we can obtain (where the term χ is 0, which is because there is no SH wave):
section_name: Conflicts of Interest: