Wave Properties of Gas-Hydrate Bearing Sediments Based on Poroelasticity

Natural gas hydrates have the properties of ice with a microporous structure and its concentration in sediments highly affects the wave velocity and attenuation. Previous studies have performed investigations based on the measurements of laboratory data, sonic-log data, and field data, whereas the variation trend of wave dissipation with increasing hydrate concentration at different frequencies is still unclear. We consider two different models to study this problem, both based on the Biot-Rayleigh double-porosity theory. In the first model, hydrate is part of the pore infill, unbonded from the grains, and brine saturates the remaining pore space. In the second model, hydrate forms a second skeleton and cements the grains. We obtain the P-wave velocity dispersion and attenuation as a function of the inclusion radius, porosity, and hydrate content. The analysis shows that the predictions of both models agree with the experimental data. At sonic log frequencies, the second model predicts much more attenuation, due to wave-induced local fluid flow (mesoscopic loss), and the behavior is such that below a given hydrate concentration the dissipation increases and then decreases beyond that concentration.


INTRODUCTION
The concentration of natural gas hydrate in sediments affects their acoustic properties (Guerin and Goldberg, 2002). In particular, the wave velocity highly increases, even with small concentrations (Chand and Minshull, 2003;Dvorkin et al., 2003;Guerin and Goldberg, 2005;Lee and Waite, 2008). This property is used to determine the distribution of hydrate in sediments (Tinivella and Carcione, 2001). Moreover, the P-wave quality factor obtained, by attenuation tomography, for instance, can also be used to monitor the presence of gas hydrate . Samples from marine sediments or from permafrost regions show that hydrates are microporous (Kuhs et al., 2004). It has been observed that seismic attenuation decreases due to a stiffening effect when the hydrate concentration increases (Dvorkin et al., 2003;Rossi et al., 2007;Gei and Carcione, 2003). Vertical-seismic-profile data in Nankai areas showed no significant attenuation (Matsushima, 2006). However, the sonic log data in the Mallik area show that attenuation increases linearly with hydrate saturation (Guerin and Goldberg, 2002;. Actually, attenuation depends not only on hydrate concentration, but also on other reservoir properties, such as microstructure (Ecker et al., 2000;Priest et al., 2009). Figure 1 shows different cases (Ecker et al., 1998(Ecker et al., , 2000Waite et al., 2009;Zhan and Matsushima, 2018): (1) hydrate grows freely in the pores away from the grains (pore-filling hydrate); (2) hydrate forms a second skeleton (load-bearing hydrate); (3) hydrate cements the grains (contact-cementing hydrate); (4) hydrate evenly deposited on the grain surface (envelope-cementing hydrate). Gas hydrate can be part of the pore infill or rock skeleton (Liu et al., 2018;Lin et al., 2019). Dvorkin and Nur (1993) proposed the BISQ model to effectively combine the global Biot-flow mechanism with the local squirt-flow mechanism. This model can be used to describe seismic attenuation in combination with frozen porous-media theories (Leclaire et al., 1994;Guerin and Goldberg, 2005). The effective-medium theory has been implemented by Helgerud et al. (1999) and Dvorkin et al. (1999). Carcione and Tinivella (2000) presented a Biot-type theory with two solids and one fluid, where hydrate forms an additional skeleton and water is the pore fluid. The concentration of gas hydrate and saturation of free gas has been estimated from seismic velocity and attenuation based on this theory Carcione et al., 2005). Lee (2002a,b) extended the classical Biot-Gassmann theory to predict the S-wave velocity. Chand et al. (2006) use the self-consistent approximation, differential effective medium theory and consider the Biot and squirt-flow attenuation mechanisms to estimate the hydrate concentration in the Mallik area. An effective-medium model, considering different microstructures, has been considered in Liu et al. (2017) and Qu et al. (2016) use a model based on pennyshaped and infinite cracks. Sell et al. (2016) combine tomography and 3D modeling to simulate the acoustic properties with digital rock physics, whereas Sell et al. (2018) show that there is an interfacial water film between hydrate and grains through 3D micro-tomography, and a new conceptual squirt-flow model is proposed. Tuan et al. (2019) uses a homogenization theory for multiphase composites to predict wave velocities that agree with laboratory and sonic-log data.
Based on the effective-grain model (Leurer, 1997), Best et al. (2013) apply the squirt-flow mechanism to explain the observed attenuation, and Marín-Moreno et al. (2017) consider squirtflow mechanisms and gas-bubble damping factors, showing that there are different behaviors of the attenuation as a function of frequency above and below a given hydrate saturation (6%) (Sahoo et al., 2019).
In this work, the Biot-Rayleigh theory of wave propagation in double-porosity media is used here to describe the gas-hydrate distribution in the sediment and study wave dispersion and attenuation (Ba et al., 2011). The model predictions are then compared with laboratory data (Priest et al., 2005(Priest et al., , 2009) and sonic-logging data (Zhan and Matsushima, 2018).

ROCK-PHYSICS MODELS
As shown in Figure 2, solid hydrate can be part of the pore infill (Model 1, Figure 2A) or constitute a skeleton and cement the grains (Model 2, Figures 2B,C). The skeleton or frame properties are obtained with the effective-medium theory according to Ecker et al. (2000), and the wet-rock P-wave properties with the Biot-Rayleigh theory. Figures 3, 4 show the modeling workflow and the description of the two models, respectively. Hydrate is abundant in the seafloor and permafrost (e.g., Marín-Moreno et al., 2017), for instance, in Qilian region (China) (Qu et al., 2016), the Mount Elbert in Alaska , and in the Eastern Nankai Trough offshore Japan (Konno et al., 2015). These cases showed the characteristics of Model 1. Lei et al. (2018) and Chen et al. (2020) observed the distribution characteristics of Model 2 by CT imaging.

Model 1: Hydrate as Pore Infill
At low concentrations, hydrate can be assumed as pore infill (Lee et al., 2010), such that a fluid is formation water and the other is a mixture of solid hydrate and free gas (hydrate/gas shortly), with ρ (1) where K w and ρ w are the water bulk modulus and density, respectively, and ρ (2) where K h and K g are the bulk moduli of gas hydrate and free gas, respectively, and c h and c g are the respective volume ratios (c h +c g = 1). ρ h and ρ g are the densities of gas hydrate and free gas, respectively. The saturation of free gas is with S h + S g + S w = 1. The bulk and shear moduli of the m-phase mineral mixture without hydrate (moduli of the mineral mixture, referred to as mineral moduli) are obtained with Hill's average: Where f i is the volumetric fraction of the i-th component and K i and G i are the respective bulk and shear moduli. The dry-rock moduli can be obtained by combining Eqs 1-6 and the modified Hashin-Shtrikman theory given in Supplementary Appendix A. By substituting a plane-wave solution into the Biot-Rayleigh equations (see Supplementary Appendix B) of Ba et al. (2011Ba et al. ( , 2012, we obtain FIGURE 1 | Microstructure in gas hydrate-bearing sediments. The white, orange, and blue colors indicate the hydrate, grain, and pore fluid. where ω is the angular frequency, k is the complex P-wave wave number,  where φ 1 and φ 2 are the porosities of the water-and hydratesaturated pores, respectively, φ 10 and φ 20 are the local porosities in the two regions, ρ f , η 1 and κ 1 are the density, viscosity and permeability of the host phase fluid, respectively. b 1 and b 2 are Biot dissipation coefficients, ρ 11 , ρ 12 , ρ 13 , ρ 22 , and ρ 33 are density coefficients, and A, Q 1 , Q 2 , R 1 , and R 2 are stiffnesses (see Supplementary Appendix C).
The mixture of hydrate/gas forms a spherical inclusion of radius of R 0 . It is related to the scale of hydrate as a solid component or a hydrate (fluid type)-saturated porous solid, which are embedded in the host rock frame. The phase velocity and quality factor are (Carcione, 2014).

Examples
The properties of the components the of hydrate-bearing sediment are given in Table 1, and other properties are shown in Table 2 (Best et al., 2013). The free gas properties in Table 1 are given in Ecker et al. (2000), and the free gas ratio is set as c g = 0.02. The energy losses caused by local fluid flow depends on the inclusion radius of the mixture of hydrate/gas in the pores. Four radii are selected to illustrate the physics. Figure 5 shows the P-wave phase velocity ( Figure 5A) and dissipation factor ( Figure 5B) as a function of frequency for different inclusion radii. Two peaks occur, corresponding to the local and global flow mechanisms. The first moves towards the low frequencies with increasing radius. The radius mainly controls the location  of the peak and has no significant effect on the global flow dispersion/attenuation, which occurs at high frequencies.
Porosity is an important reservoir property that affects the dry-rock moduli. The porosity of core samples in gas hydratebearing sediments of the Nankai Trough ranges between 35 and 43% (Zhan and Matsushima, 2018), while the porosity in Qilian Mountain permafrost is much smaller, with values between 1 and 5% (Lin et al., 2019). Figure 6 shows the P-wave phase velocity ( Figure 6A) and dissipation factor ( Figure 6B) as a function of frequency for different porosities. Velocity decreases and the local and global-flow attenuations increase with increasing porosity. At 10% porosity, these loss mechanisms start to be evident and at higher porosities the two peaks split, with the global-flow one more affected. Figure 7 shows the P-wave phase velocity ( Figure 7A) and dissipation factor ( Figure 7B) as a function of frequency for different hydrate concentrations. Model 1 is assumed for hydrate content less than 40%. Velocity increases with hydrate concentration. The local-and global flow peak amplitude increase with concentration. The locations of peaks are hardly affected. Basically, Figure 7B shows that attenuation depends on hydrate content and frequency.

Model 2: Hydrate as an Additional Skeleton
In this case, the pore fluid is a mixture of water and free gas, such that the effective bulk modulus is where S g +S w = 1. The density of the pore fluid is The porosity, after the gas-hydrate deposition can be obtained as (Helgerud et al., 1999;Ecker et al., 2000) and the volume percentage of hydrate in the solid matrix is According to Eq. 6, the mineral moduli considering the presence of hydrate are where K n and G n are the bulk modulus and shear modulus of the mineral mixtures without hydrate, respectively. Clay is a solid component in the theory of Ba et al. (2016), i.e., and the host skeleton is composed of grains and intergranular pores, forming two solids and one fluid (water/gas). The clay squirt-flow mechanism induced by acoustic waves is discussed in that paper. Compared with hydrate as a pure solid, the hydrate micropores and intergranular pores are connected to favor squirt flow. Here, gas hydrate is modeled as clay, and it is considered to also cause wave-induced fluid flow and dissipation (Best et al., 2013;Marín-Moreno et al., 2017;Sahoo et al., 2019). X-ray CT analysis has observed the contact and cementation of gas hydrate and mineral grains in the core samples, supporting the setup of skeleton type model (Lei et al., 2018).
The dry-rock moduli are also determined with the modified Hashin-Shtrikman theory given in Supplementary Appendix A and the wave propagation by Eq. 7, redefining the variables and properties as in Supplementary Appendix C (Ba et al., 2011;Sun et al., 2015). φ 10 and φ 20 are the local porosities of the main skeleton and inclusions, and φ 1 and φ 2 are the absolute porosities of the main skeleton and inclusions, respectively. The hydrate skeleton forms a spherical inclusion with a radius of R 0 .

Examples
The inclusion radius of hydrate, which contains micropores, also affects the squirt flow. We choose the same four inclusion radii for the comparison with Model 1. Figure 8 shows the P-wave phase velocity ( Figure 8A) and dissipation factor (Figure 8b) as a function of frequency for different inclusion radii. If the radius increases, the local fluid-flow attenuation peak moves to the low frequencies. The global fluid-flow peak, occurring at high frequencies, is much weaker, almost negligible. Figure 9 shows the P-wave phase velocity ( Figure 9A) and dissipation factor ( Figure 9B) as a function of frequency for A B FIGURE 5 | Model 1. P-wave velocity (A) and dissipation factor (B) as a function of frequency for different inclusion radii of hydrate/gas. The hydrate concentration is 20%, and the porosity is 40%. different porosities. As in Model 1, increasing porosity implies decreasing P-wave velocity and increasing local-flow attenuation. Beyond 40%, two attenuation peaks can be observed, with the weaker one at high frequencies being the global-flow one. Figure 10 shows the P-wave phase velocity ( Figure 9A) and dissipation factor ( Figure 9B) as a function of frequency for different hydrate concentrations. The P-wave velocity increases with the concentration and the local-flow peak dissipation has a maximum value at 30% and then decreases, while the global-flow peak is weak. The first peak is mainly located at seismic frequencies.

Laboratory Data
We have used the experimental data from Priest et al. (2005Priest et al. ( , 2009 and Best et al. (2013), who performed the resonance-column experiments on synthetic hydrates, to further analyze the two models. Two methods, termed "excess water" and "excess gas, " are implement to generate solid hydrate (Priest et al., 2009). P-wave velocity and attenuation at an effective pressure of 500 kPa were measured in sand samples (Priest et al., 2005;Priest et al., 2009;Best et al., 2013). The moduli, permeability, and critical porosity are given in Tables 1, 2. In section "Laboratory Data, " we choose the same grain coordination number n = 4 as Best et al. (2013), to compare the theoretical curves proposed by those authors.

Excess-Water Method
The continuously injected water reacts with a certain amount of methane gas to form small quantities of hydrates, which are not in contact with the grains. In this case, we have the assumption of Model 1 (Best et al., 2013;Tuan et al., 2019), with hydrate concentrations less than 40% (Priest et al., 2009;Zhao et al., 2015). Figures 11, 12 compare the theoretical and experimental P-wave velocities and dissipation factors as a function of hydrate concentration around 200 Hz, respectively. Porosities are 40 A B FIGURE 7 | Model 1. P-wave velocity (A) and dissipation factor (B) as a function of frequency for different hydrate contents. The inclusion radius is 1 cm and the porosity is 40%. and 42%, respectively. The inclusion radius is 0.8 and 4 cm, respectively, based on fitting. We assume a small free-gas volume ratio (0.002), considering that the free gas does not fully react with water during the experiment. As the hydrate content increases, it gradually cements the mineral grains and we have Model 2. Chen et al. (2020) applied microtomography (CT) to analyze the evolution of hydrate pore habit during hydrate formation. Therefore, at low concentrations, the P-wave velocity increases slowly with hydrate content (Ecker et al., 1998).
Both the pore-filling model and cementing model proposed by Priest et al. (2009) are based on effective-medium theories, because Gassmann equation cannot describe the characteristics of rocks with a double-porosity structure. Therefore, Model 1 can better describe the P-wave velocity, as shown in Figure 11. As the hydrate content exceeds 30%, Model 1 does not fit the P-wave velocity and the assumptions of Model 2 hold. Best et al. (2013) used the effective-medium theory and the correspondence principle to model attenuation (poro-filling curve shown in A B FIGURE 8 | Model 2. P-wave velocity (A) and dissipation factor (B) as a function of frequency for different inclusion radii of hydrate/gas. The hydrate concentration is 20% and the porosity is 40%. Figure 12). Their prediction results are lower than the measured attenuation values, while Model 1 provides a relatively good fit.
The key parameter is the size of the hydrate-inclusion radius.

Excess-Gas Method
A high amount of methane gas is injected to react with water, and then gas hydrate is formed at the grain contacts or at the surface of the grains. Thus, we have the conditions of Model 2 (Best et al., 2013;Tuan et al., 2019). Figures 13, 14 compare the theoretical and experimental P-wave velocities and dissipation factors as a function of hydrate concentration around 200 Hz, respectively (the porosities are 42%). The hydrate inclusion radius are 2 and 0.2 cm, respectively, and as in section "Excess-Water Method, " the free-gas volume ratio is 0.002.
The hydrate cementing the grains significantly increases the dry-rock moduli and the P-wave velocity. Although the P-wave velocity of the cementing model proposed by Priest et al. (2009) does not fit the experimental data (Figure 13), it follows the trend A B FIGURE 9 | Model 2. P-wave velocity (A) and dissipation factor (B) as a function of frequency for different rock porosities. The inclusion radius is 1 cm, and hydrate concentration is 20%.
at low hydrate content. The results from Model 2 show that the P-wave velocity increases sharply at low hydrate saturations. It is evident that Model 2 better describes the P-wave velocity data. In Figure 14, the load-bearing curve proposed by Best et al. (2013) agrees with the data when the hydrate content is lower than 10%, but the match of the attenuation peak of attenuation is not good.
Instead, Model 2, shows a better agreement.

Sonic-Logging Data
In the hydrate-bearing sediments, estimations of hydrate concentration based on a single model may not be appropriate, because of the different spatial distributions of the hydrate. For instance, Liu et al. (2020) confirmed the coexistence of different gas-hydrate distributions from seismic velocity.
We have used the sonic-log data obtained by Zhan and Matsushima (2018) in Nankai Trough of Japan for the analysis, data required by the two models. In their work, the sonic-log (magnetic monopoles) frequency is 14 kHz, the average porosity is 40%, the grain coordination number is 8.5, the seawater viscosity is 0.0018 Pa s, and the effective pressure is about 500 kPa. The hydrate inclusion radius is 0.18 cm. Figure 15 shows the results. The sonic-log data in Figure 15A shows that A B FIGURE 10 | Model 2. P-wave velocity (A) and dissipation factor (B) as a function of frequency for different hydrate contents. The inclusion radius is 1 cm and the porosity is 40%.
the P-wave velocity gradually increases with hydrate content, showing a linear relation. Zhan and Matsushima (2018) used the Marín-Moreno and Guerin models to predict P-wave velocity and attenuation, respectively. When the hydrate saturation is higher than 20%, the Guerin model is consistent with the P-wave velocity of sonic-logging data. The Marín-Moreno model considers 20% pore-filling hydrate and 80% contact cementing hydrate, giving a better prediction. Without considering the diverse hydrate morphologies, the Models 1 and 2 are applied to predict the velocities in the lower and higher bounds, respectively. Figure 15B, shows that the experimental dissipation factor is distributed on a wide range, which may be related to different geometrical distributions of the hydrate. Model 1 predicts low values, while those of Model 2 are more consistent with the data. The attenuation of P-wave increases first and then decreases with hydrate saturation increasing. Rossi et al. (2007) estimated a lower dissipation factor when hydrate is present at lower frequencies than those involved in Figure 15B. The increasing    attenuation at low concentrations can be due to a local-flow mechanism or to scattering loss, since the data corresponds to sonic frequencies.

CONCLUSION
We propose two models to calculate the wave velocity and attenuation of gas-hydrate bearing sediments, based on the Biot-Rayleigh double-porosity theory and two different distributions of the hydrate in the porous medium. The difference is that in one case hydrate is part of the pore infill and in the other constitutes an additional (load-bearing) skeleton. The local inclusion radius of the Biot-Rayleigh theory is related to the hydrate/gas mixture in the first model and to the hydrate frame in the second. The models predict two relaxation peaks, namely the localflow one (mesoscopic loss) and the classical global flow one, also predicted by the Biot theory. As expected, P-wave velocity decreases and attenuation increases with porosity. The localflow relaxation (attenuation) peak moves to the low frequencies with increasing inclusion radius. More hydrate content has the effect to increase the attenuation at low values, reach a maximum loss and then decrease it. Modeling results are compared with data obtained from the excess-water and excess-gas experiments,