The Influence of Gas Hydrate Morphology on Reservoir Permeability and Geophysical Shear Wave Remote Sensing

We show that direct estimates of the permeability of hydrate‐bearing geological formations are possible from remote measurements of shear wave velocity (Vs) and attenuation (Qs−1). We measured Vs, Qs−1 and electrical resistivity at time intervals during methane hydrate formation in Berea sandstone using a laboratory ultrasonic pulse‐echo system. We observed that Vs and Qs−1 both increase with hydrate saturation Sh, with two peaks in Qs−1 at hydrate saturations of around 6% and 20% that correspond to changes in gradient of Vs. We implemented changes in permeability with hydrate saturation into well‐known Biot‐type poro‐elastic models for two‐ and three‐phases for low (Sh < 12%) and high (Sh > 12%) hydrate saturations respectively. By accounting for changes in permeability linked to hydrate morphology, the models were able to describe the Vs and Qs−1 observations. We found that the first Qs−1 peak is caused by a reduction of permeability during hydrate formation associated with a transition from pore‐floating to pore‐bridging hydrate morphology; similarly, the second Qs−1 peak is caused by a permeability reduction associated with a transition from pore‐bridging hydrate morphology to an interlocking network of hydrate in the pores. We inverted for permeability using our poro‐elastic models from Vs and Qs−1. This inverted permeability agrees with permeability obtained independently from electrical resistivity. We demonstrate a good match of our models to shear wave data at 200 Hz and 2 kHz frequencies from the literature, indicating the general applicability of the models.

Accurate estimates of the amount, distribution and permeability of methane hydrate deposits are essential for predicting greenhouse gas fluxes from hydrate to the atmosphere and their impact on future climate change. Such estimates are also needed for assessing the role of hydrates in geohazards, and their energy resource potential. It is well known that the presence of hydrate in sediments changes elastic wave velocity and attenuation, electrical resistivity and intrinsic permeability; the extent of change depends on the amount (saturation) and distribution (morphology) of hydrate (Attias et al., 2016;Bauer et al., 2015;Chen et al., 2014;Crutchley et al., 2010;Ecker et al., 1998;Han et al., 2019;Helgerud et al., 1999;Kwon et al., 2011;Lee, 2008;Lee & Collett, 2013;Nimblett & Ruppel, 2003;Ojha et al., 2010;Pecher et al., 1996;Sahoo, Chi, et al., 2018;Singhroha et al., 2020;Spangenberg & Kulenkampff, 2006;Vadakkepuliyambatta et al., 2015).
Methane hydrate exhibit different spatial distribution (i.e., morphology or habit) in sediment pores (e.g., Dai et al., 2012;Holland et al., 2008). In sandy host medium, hydrate morphology can be broadly defined as cementing or non-cementing hydrate (e.g., Ecker et al., 1998). Cementing hydrate forms at the interfaces between sediment grains and bonds them together, increasing the strength of the host sediment. Non-cementing hydrate grows away from the sediment grain contacts (e.g., Ecker et al., 1998) and be further sub-divided into pore-floating or pore-bridging. Pore-floating (or pore-filling) hydrate contained wholly within the pore fluids in the sediment pore space without bridging neighboring sediment grains, and only directly changes the bulk average pore fluid properties. Pore-bridging (or load-bearing or frame-supporting) hydrate bridges neighboring sediment grains (i.e., contacts more than one grain in the sand frame), so influences both pore fluid and host sediment properties.
These different hydrate morphologies were initially deduced from the effect of hydrate morphology on geophysical properties in the 90 s (e.g., Ecker et al., 1998). Some of these terms can be bit confusing for broad geophysical readers; for example, load-bearing and cementing can be a bit confusing as cementing hydrate can also bear load. Similarly, pore-filling can also refer to load bearing or cementing hydrate or grain coating hydrate as all these can fill some part of pore space. Recent studies showing high resolution imaging of hydrate morphologies have helped to describe these morphologies in a more descriptive way (e.g., Chaouachi et al., 2015;Hu et al., 2014;Kerkar et al., 2014;Tohidi et al., 2001). So based on these imaging studies of hydrate and to help the readers of broad background we use the term "pore-floating" and "pore-bridging." Recent high resolution synchrotron imaging has shown how pore-bridging hydrate can coalesce to form an extensive network of hydrate or so-called "inter-pore hydrate framework morphology" . These different hydrate morphologies affect the permeability in different ways (e.g., Wang et al., 2021). The pore-floating hydrate is part of the pore fluid and leads to only a small reduction in permeability of the effective medium, while pore-bridging hydrate causes a larger reduction in permeability as it gives rise to more tortuous pathways for fluid flow. The inter-pore hydrate framework causes the most significant reduction in permeability compared to the first two morphologies. The inter-pore hydrate framework behaves as a single solid phase which interlocks with the sediment solid phase, and results in a significant increase in shear coupling between the two solid phases, and hence an increase in shear wave velocity . Another aspect of gas hydrate which affects the physical properties is homogenous or heterogeneous distribution of hydrates in sediments at larger scale within the area of interest (e.g., Dai et al., 2012). All the morphologies described above can be either homogenously or heterogeneously distributed within the area of interest.
Elastic wave velocity is used most commonly to remotely quantify gas hydrate deposits from seismic geophysical surveys (Ecker et al., 2000;Helgerud et al., 1999;Lee et al., 1993;Yuan et al., 1996). Rock physics models are used to relate the effect of gas hydrate content and distribution to wave velocity. However, the most commonly occurring hydrate morphologies (pore-floating and pore-bridging) only have a significant effect on velocity at hydrate saturations above 25%-40% (e.g., Waite et al., 2009). Elastic wave attenuation may be used as an alternative geophysical property for gas hydrate quantification (Bellefleur et al., 2007;Dvorkin & Uden, 2004;Guerin & Goldberg, 2002;Jyothi et al., 2017). Hydrate quantification from attenuation has uncertainties due to our limited understanding of attenuation mechanisms in sediments in general, and of the role of hydrate in particular (e.g., Best et al., 2013;Dvorkin & Uden, 2004;Guerin & Goldberg, 2002;Marín-Moreno et al., 2017). It is often challenging to relate in-situ attenuation measurements to gas hydrate content and distribution due to uncertainties like spatial averaging effects 10.1029/2021JB022206 3 of 19 (e.g., Bellefleur et al., 2007;Matsushima, 2006). Laboratory studies may offer useful insights, but there are few published laboratory measurements of wave attenuation on hydrate bearing sediments (e.g., Best et al., 2013;Priest et al., 2006;Sahoo et al., 2019).
The intrinsic permeability of hydrate bearing sediments is an important property that affects hydrate content and distribution, methane flux into seawater, geophysical response and gas production (e.g., Ren et al., 2020). Sediment permeability is controlled by connectivity, tortuosity and size distribution of pores and fractures (if present). Most of the geophysical rock physics models based on Biot's poro-elastic and squirt flow formulations (Biot, 1956a(Biot, , 1956b also use permeability as an input parameter (e.g., Carcione & Tinivella, 2000;Chand et al., 2006;. Studies of permeability of hydrate bearing sediments have used careful laboratory fluid flow measurements and X-ray computed tomography on synthetic hydrates and preserved natural hydrates cores (e.g., Konno et al., 2013). The permeability of hydrate bearing sediment has also been estimated from nuclear magnetic resonance (NMR) methods during wireline logging (Jain et al., 2019;Kleinberg et al., 2003;Lee, 2008;Li, Zhao, et al., 2014;Uchida & Tsuji, 2004). However, it is difficult to estimate permeability in the absence of cores or NMR logging data. A pressure gradient is created when elastic waves pass through a partially saturated porous medium, which in turn causes relative fluid flow between the solid and fluid phases. The ease of fluid flow depends on the permeability. This effect raises the possibility of estimating permeability from seismic measurements, which can provide information at much larger scales than core or log analysis. Klimentos and McCann (1990) studied the effect of permeability on P-wave velocity and attenuation using laboratory measurements and found attenuation is more sensitive to permeability than velocity. Best et al. (1994) made laboratory measurements on several sandstone samples and found that P and S wave attenuation are more responsive than velocity to changes in permeability. There are only a handful of studies on the effect of fluid permeability on elastic properties, and none specific to hydrate (Akbar et al., 1993;Berryman et al., 2016;Goloshubin et al., 2008;Klimentos & McCann, 2002;Kwon & Ajo-Franklin, 2013;Martin, 1996;Pride et al., 2003Pride et al., , 2004Rubino et al., 2012) and all have concluded the relationship to be very complex and to require more work.
So, to summarize, there are studies on how gas hydrate affects permeability and separate studies on the effect of hydrate on geophysical properties. Geophysical properties of fluid saturated porous media depend on fluid permeability. Hence, we aim to combine these models and discover whether the permeability of hydrate bearing sediments can be inferred from geophysical properties. To this end, we analyzed experiment data to see if the permeability of gas hydrate bearing sandstone can be inferred from shear wave measurements. We performed consecutive cycles of methane hydrate formation and dissociation in a Berea sandstone sample, and obtained novel laboratory measurements of shear wave attenuation; we also measured electrical resistivity and used it to calculate permeability. We calculated the hydrate saturation independently using a thermodynamic method with pressure and temperature measurements as input. In this study, we provide insight into the likely dominant shear wave propagation mechanisms. We were able to link the observed changes in S-wave attenuation to inferred permeability changes by adapting Biot poro-elastic theory with our knowledge of different hydrate morphologies and how hydrate saturation affects permeability. The resulting models provide a way to interpret the permeability of hydrate deposits from S-wave seismic survey and well log data.

Hydrate Formation
We formed methane hydrate in a cylindrical sample of Berea sandstone (4.97 cm diameter, 2.06 cm height; 22% porosity, 448 mD permeability) inside a stainless-steel pressure cell ( Figure S1 in Supporting Information S1). Berea sandstone is as an inert, stable and well-characterized porous medium typical of geological hydrocarbon reservoirs. We oven-dried the sample at 60°C, and placed in a rubber sleeve to prevent direct contact between the mineral oil used as confining fluid and the rock sample. A temperature sensor was also attached to the outer surface of rubber sleeve. We then placed the sample including the rubber sleeve inside the pressure cell and then subjected it to hydrostatic confining pressure. We connected the pore fluid pipe to (a) a syringe pump containing a 35 g/L NaCl solution in deionized-deaerated water, (b) a vacuum pump and (c) a CH4-bottle pressurized at 12 MPa ( Figure S1 in Supporting Information S1). We can switch between 4 of 19 these inlets as needed during hydrate formation. Then we applied a vacuum up to 1 Pa on the pore fluid line to remove air from the pore space. We then injected brine (35 g/L NaCl) to occupy 83.5% of the pore volume using a syringe pump, giving excess water conditions for hydrate formation (Priest et al., 2009;. We left sample in these conditions for three days for the brine to distribute in the pore space. We then injected methane gas from a methane gas cylinder via a valve and kept the cylinder open to reach a stable pore fluid pressure ( Figure S1 in Supporting Information S1). We then sealed the pore-fluid system by closing valve A. Then, we formed hydrate by cooling the sample to 5°C to take the system into gas hydrate stability conditions and above the freezing point of water (Figure S1b in Supporting Information S1). We measured pore fluid pressure and temperature (sample and room), P-and S-wave velocity and attenuation, and electrical resistivity during hydrate formation (Figure 1 and Figures S1 and S2 in Supporting Information S1). We calculated the hydrate saturation using a thermodynamic method with measured pressure and temperature changes as input . We inferred hydrate formation was complete by observing the onset of constant values of pore pressure ( Figure S1b in Supporting Information S1) and velocity (P and S wave) following their decrease and increase during hydrate formation (e.g., Waite et al., 2004). Then, we increased the temperature to take the sample out of hydrate stability condition and dissociate hydrate. Four consecutive cycles of hydrate formation and dissociation allowed us to form uniform hydrate distributions (Choi et al., 2014;Spangenberg et al., 2015). The differential pressure (confining minus pore pressure) affects elastic wave properties (e.g., Green & Wang, 1994;Prasad & Manghnani, 1997). So, when pore pressure changed during hydrate formation or dissociation, we changed the confining pressure accordingly to maintain a constant differential pressure in each cycle. We kept a differential pressure of 10 MPa in Cycles 1 and 2, and 55 MPa in Cycles 3 and 4. The presence of micro cracks in sandstone affects elastic wave velocity and attenuation, and these effects are insignificant at high differential pressure. Therefore, we use 55 MPa differential pressure as microcracks are closed at such high differential pressures (e.g., Green & Wang, 1994;Prasad & Manghnani, 1997) and hence any changes in attenuation can be attributed to hydrate formation and dissociation processes only.

Elastic and Electrical Measurements
We used the pulse-echo system and signal analysis described in detail in Sahoo et al. (2019), after Winkler and Plona (1982). The sample was sandwiched between two cylindrical buffer rods of the same diameter as the sample, that electrically isolate the sample and provide a time delay for sample elastic wave reflections to be isolated from extraneous multiple reflections. An ultrasonic transducer with nominal frequency of 1 MHz was placed in contact with the free end of one buffer rod. We used the same transducer to send and receive the reflections from the Perspex-rock interfaces. We separated reflections from top and bottom of the sample using time domain gating and then converted them into the frequency domain using a fast Fourier transform. We then took a spectral ratio of these two signals at each frequency. Our method then diverges from the method of Winkler and Plona (1982) by employing a 1-D frequency-domain forward model of the experimental system. This forward model is then used in a nonlinear optimization scheme to obtain a best fit complex velocity at each discrete frequency in the FFT spectra obtained in the previous step. Velocity and attenuation are then calculated from the complex velocity. Our calculation method is accurate for inelastic materials unlike the Winkler and Plona (1982) approximate analytical solution, which is strictly only applicable to elastic (zero loss) materials. Winkler and Plona (1982) discussed that realistic materials are inelastic (non-zero loss) and their reflection coefficient is complex (due to loss or attenuation); however, they did not use complex reflection coefficient leading to some error. The ultrasonic transducer is of a finite size and this causes diffraction-like beam spreading; we correct for this using method of Papadakis (1972) and Benson and Kiyohara (1974). S-wave velocity (V s ) and attenuation (Q s −1 ) data have an accuracy of ±0.3% and ±10% respectively, using an ultrasonic pulse-echo system (Winkler & Plona, 1982).
We measured electrical resistivity using the system of North et al. (2013). The rubber sleeve enclosing the sample was perforated with 16 electrodes arranged in two rings around the sample ( Figure S1 in Supporting Information S1). The relative error in resistivity measurements is <0.1% under typical operating conditions (at frequencies 1-500 Hz) for homogenous and isotropic samples in the electrical resistivity range 1-100 Ωm (North et al., 2013). Some examples of recorded waveforms are shown in Figure S2 in Supporting Information S1, showing clearly identifiable top and base sample echoes with no internal reflections. A heterogeneous sample ), (e) shear modulus and (f) electrical resistivity during consecutive cycles of methane hydrate formation in Berea sandstone. Hydrate dissociation data is not shown here. Shear modulus is calculated from shear wave velocity and density. Error in electrical resistivity is <0.1% (North et al., 2013). The accuracy of velocity measurement is up to ±0.3% and attenuation is ±0.2 dB cm −1 (Best et al., 1994). P-wave data were previously reported in Sahoo et al. (2019). would give rise to internal reflections, not observed in our data ( Figure S2 in Supporting Information S1). The experimental methods, along with P-wave velocity and attenuation results, were reported in Sahoo et al. (2019). We also measured all these properties during hydrate dissociation (Sahoo, 2018), not discussed here. We report novel findings concerning the V s and Q s −1 results at a frequency of 648 kHz, although a full data set exists for 400-800 kHz (Sahoo et al., 2019). Biot (1956aBiot ( , 1956b formulated a rock physics model for wave propagation through a fluid saturated porous medium which calculates frequency dependent compressional and shear wave velocity and attenuation. This model is a two-phase model where one phase is the porous medium and the second phase is the pore fluid. This model accounts for inertial motion between these two-phases with viscous fluid flow between the two-phases giving rise to elastic wave energy dissipation. Leclaire et al. (1994) extended Biot's global fluid flow model for partially frozen porous media with three-phases. Their model considers relative inertial motions and resulting viscous losses between (a) the host porous medium, (b) ice and (c) unfrozen water in the pore space (see Supporting Information S1). This model was applied to hydrate-bearing sediments by replacing ice with hydrate (e.g., Carcione & Tinivella, 2000;Zhan & Matsushima, 2018). Leclaire et al. (1994) assumed there is no contact between the ice and the host porous medium that are always separated by a film of water. Carcione and Tinivella (2000) considered some contact between hydrate and mineral grains, and also considered the viscous dissipation between solid and fluid which was not considered in Leclaire's model.  introduced squirt flow, and a coupling shear modulus between mineral grains and hydrates. Like Best et al. (2013), we implement the  version of Leclaire's model, denoted as GG.

Biot Poro-Elastic Theory and Permeability
With our much better knowledge of likely hydrate morphology changes during hydrate formation from Sahoo et al. (2019), and our well-constrained ultrasonic data set, we use this model to gain insight into the most important physical phenomena controlling S-wave propagation in hydrate-bearing porous media: permeability and tortuosity. Once established and verified, the theory would be of considerable benefit to geological hydrate reservoir characterization using seismic and borehole sonic geophysical remote sensing methods.

Permeability
Fluid permeability in hydrate-bearing sediments is measured and defined as a bulk or effective permeability ( eff E k ) of the composite medium comprising hydrate, sediment and fluid phases (e.g., Johnson et al., 2011). However, Leclaire's model conceptually defines a separate permeability for the sand frame ( s E k ) and the hydrate frame ( h E k ).
The permeability for relative flow of pore-fluid in the sand frame, with no viscous friction from the hydrate matrix, s E k , is given by Leclaire related s E k to the intrinsic or absolute permeability of water saturated sand frame 0 s E k using the Kozeny-Carman relationship: where w E S is water saturation defined as the ratio of volume of water to volume of pore space. Similarly, h E k is defined as the permeability of the hydrate matrix (i.e., the patchwork of hydrate grains residing within the sand pores) for a water flow in the hydrate matrix with no viscous friction from the sand frame according to where  E is porosity of sample. In practice, it is difficult to measure 0 h E k because the hydrate matrix cannot exist separately from the sand frame. Nevertheless, several studies use Leclaire's model with Equation 4 to calculate k h using arbitrary values of k h0 that vary over an order of magnitude; for example, 10 −5 m 2 Leclaire et al., 1994), 5 × 10 −4 m 2 (Carcione & Tinivella, 2000), 10 −8 -10 −10 m 2 (Guerin & Goldberg, 2002), 10 −7 m 2 (Zhan & Matsushima, 2018). Guerin and Goldberg (2002) and Zhan and Matsushima (2018) showed that both Q p −1 and Q s −1 strongly depend on k h0 . Also, the value of k h0 depends on hydrate morphology and how hydrate morphology changes with hydrate saturation. Leclaire et al. (1994) also related k h and k s to the effective permeability in a partially frozen medium, which we can extend to a partially saturated hydrate-bearing sediment according to Several equations exist to calculate the effective permeability of hydrate-bearing sediments Masuda, 1999)  where arbitrary values of k h0 are needed.

Geometrical Aspect Ratio of Boundary of Different Phases
All versions of Leclaire et al.'s (1994) three-phase Biot model consider energy loss due to tortuosity, parameterized by the geometrical aspect ratio of the boundary between different phases. The GG model termed this as an inertial coupling coefficient but we will stick to the original geometrical definition of Leclaire et al. (1994) as this term seem to have a more physical meaning. The geometrical aspect ratios r sh,, r sf, r hf depend on the shape of the physical boundary that separates the phases; subscripts sh, sf, hf, refer to the sand/hydrate, sand/fluid and hydrate/fluid boundaries, respectively. The value of r can be calculated from the geometry or shape of the boundary separating the two-phases (Lamb, 1945). The geometry of a boundary of hydrate or sand in 3D, can be approximated by an ellipsoid: 8 of 19 Lamb (1945) calculated r for changes in the ratio a/b, as outlined in (Table S1 in Supporting Information S1). For a spheroid (a = b = c), r is 0.5, and for prolate ellipsoids (a > b, b = c), r decreases to zero as it becomes narrower (a/b → ∞). The lower the value of r, the greater is the surface area of the boundary, resulting in more attenuation. Leclaire et al. (1994) and Carcione and Tinivella (2000) used a spherical boundary with r hf = r sh = r sf = 0.5.  showed that maximum attenuation is obtained when r hf is zero for the hydrate/ fluid boundary. They tested the model of Leclaire et al. (1994) and a modification of the Leclaire et al. (1994) model by Carcione and Tinivella (2000) on borehole data from the Mallik 5L-38 gas hydrate well. They found that r hf = 0.5 matched well the sonic log P-and S-wave velocities, but gave very low attenuation (Q p −1 and Q s −1 ). However, they were able to obtain a good match for V p , V s , Q p −1 and Q s −1 by adding cementation, friction, squirt flow (specifically the BISQ model of ) and using r hf = 0 into the model (Figure 6 of . We next apply the Biot two-phase model and three-phase model (GG model) to our ultrasonic S-wave data to obtain consistent predictions for V s and Q s −1 based on our knowledge of different hydrate morphologies and their effect on model permeability and inertial coupling values. We explore a range of possible and plausible values. All the input parameters are given in Table S2 in Supporting Information S1.

S-Wave Velocity Versus Hydrate Saturation
We observe V s increases as hydrate forms (Figure 1a). Here, we present data from Cycle 3 and four in Figure 1a as ΔV s (change in V s relative to V s at zero hydrate saturation), which show similar results, although we would expect the hydrate to be more homogeneously distributed in Cycle 4 than in Cycle 3 (data presented in see above).
We see that ΔV s has no significant increase up to about 6% hydrate saturation (S h ), then it increases linearly to about S h = 12% where the curve is inferred to flatten (interpolation between available data points) before increasing rapidly at about S h = 20%, after which it increases slowly to the experimental maximum hydrate saturation of 25%. As noted by their Figure 9), the observed pattern of increasing V s with hydrate saturation for Cycle 3 cannot be explained by the transition from a pore-floating to a pore-bridging hydrate morphology that adequately describes the P-wave velocity results (according to the HBES effective medium model of Marin-Moreno et al., 2017).  were able to image changing hydrate morphology during an analogous experiment in sand using synchrotron X-ray CT, and concluded that a third hydrate morphology was likely controlling V s at higher hydrate saturations, namely an interpore hydrate framework. The weaker-than-expected effect of pore-bridging hydrate on V s was attributed to the presence of water films between the Berea sandstone grain framework and the interpore hydrate grains, consistent with three-phase Biot poro-elastic theory Leclaire et al., 1994).
Based on these observations, we now apply the Biot two-phase model for low hydrate saturations, where hydrate behaves as part of the pore fluid. We use the Biot three-phase model at higher saturations, when we treat the interpore hydrate framework morphology as a third phase.

S-wave Attenuation Versus Hydrate Saturation
In a series of resonant column laboratory experiments at seismic frequencies (200 Hz), Best et al. (2013) noted that, unlike P-wave attenuation, S-wave attenuation in methane hydrate-bearing sand (under similar excess water conditions to the ultrasonic experiments presented here) could not be explained by a squirt flow loss mechanism modeled theoretically (HEG model) for pore-floating and pore-bridging hydrate morphologies. There was, however, a better match of the S-wave attenuation to the three-phase Biot theoretical model. Figure 6h of Best et al. (2013) shows two attenuation peaks at about 10% and 35% hydrate saturations respectively; the three-phase Biot model of  was able to match the second peak quite convincingly using well constrained model input parameters. As noted by Best et al. (2013), the V s and Q s −1 observations suggest an underlying global viscous fluid flow loss mechanisms consistent with three-phase Biot theory originally developed for ice in sediments by Leclaire et al. (1994). Although the observations were made at disparate frequency ranges (seismic and ultrasonic), their similarities provoke further investigation into whether the same loss mechanism is operating in both cases.

of 19
Here, we present our novel ultrasonic S-wave attenuation observations in Figure 1b as changes in Q s −1 with respect to Q s −1 at zero hydrate saturation, denoted ΔQ s −1 . Similarly, the results of Best et al. (2013) show evidence for two attenuation peaks in ΔQ s −1 with increasing hydrate saturation. These peaks in ΔQ s −1 seem to coincide with changes in gradient of ΔV s. in Figure 1a.
There is a small peak in ΔQ s −1 at about S h = 6%, then missing data (due to night/weekend health and safety laboratory working restrictions), then a few points at about 20% and 23%-24% hydrate nevertheless indicating a much larger peak in ΔQ s −1 at about S h = 20%. These peaks are also seen in both cycles around similar hydrate saturations, and they seem to correspond to the hydrate morphology transitions from pore-floating to pore-bridging at about S h = 6%, and from pore-bridging to interpore hydrate framework at about S h = 20%. In the next section, we will use the Biot theoretical models discussed in Section 3 to verify if these peaks correspond to hydrate morphology transitions affecting permeability and tortuosity parametrized using geometrical aspect ratios or inertial coupling coefficients.

Ultrasonic Frequency Measurements
Formation of hydrate in the pore space results in a decrease of effective permeability k eff , and the amount of decrease depends on hydrate morphology (location of hydrate in the pores) and saturation (e.g., . Hydrate morphologies that block connectivity between pores decrease k eff more than hydrate morphologies that only constrict permeability within pores. That is, the former can switch off certain pore network pathways, while the latter can reduce flow through the entire network of pores more uniformly. We see two peaks in attenuation around S h = 6% and 20% in the ultrasonic data ( Figure 1). The small peak in ΔQ s −1 at about S h = 6% can be linked to a transition from pore floating to pore bridging morphology, while the second peak at S h = 20% corresponds to inter-pore hydrate framework formation. We first use the Guerin-Goldberg (GG) version of the Biot three-phase model to match the bigger attenuation peak at S h = 20% (Figure 2). We clearly see that if we apply the GG model without incorporating the effect of permeability variation with hydrate formation, attenuation is under-predicted by the model. GG models with permeability parameters N = 5 or M = 2 match the second peak shear wave attenuation ΔQ s −1 observed at S h = 20%, suggesting an optimum arrangement of hydrate in the pores for Biot-type attenuation (an interplay between permeability, fluid viscosity, inertia, and wave frequency; see Section 3.3). Physically, this corresponds to the frame suddenly becoming stiffer while permeability reduces only slightly, so it affects Biot fluid flow for S-waves most, as evident from the increase in shear modulus around this second Q s −1 peak at S h = 20% (Figure 1). The shear modulus does not show any sudden increase around the first Q s −1 peak, which supports our hypothesis that the second Q s −1 peak is due to a sudden stiffening of sample.
It is interesting to note that while Biot's three-phase modeled ΔQ s −1 is very sensitive to permeability changes, modeled ΔQ p −1 is much less dependent on permeability (Figure 1 and Figure S3 in Supporting Information S1). These modeling results suggest that the increase in experimentally measured ΔQ p −1 with S h is likely dominated by other mechanisms, and not so much by permeability. Sahoo et al. (2019) found the dominant mechanism for increase in ΔQ p −1 with S h is squirt flow in different aspect ratio pores created due to hydrate formation. Although the magnitude of change in ΔQ s −1 at the first peak around S h = 6% is small, we can still see that when Biot's two-phase model incorporating the Tokyo permeability method is applied with pore floating morphology, ΔQ s −1 increases and ΔV s remains almost constant, matching with data up to S h < 6% (inset of Figures 2c and 2d). When hydrate morphology changes to pore-bridging, ΔQ s −1 decreases and ΔV s increases as seen in the data for S h > 6%. This shows that the peak in ΔQ s −1 at S h = 6% is likely due to a transition from pore-floating to pore-bridging hydrate morphology.
If we do not account for permeability changes with hydrate formation, the Biot two-phase model calculates ΔQ s −1 -0 for both pore-floating and pore-bridging morphologies. We note that the Biot two-phase model incorporating the Tokyo permeability method with pore-bridging hydrate morphology still shows a lower velocity than the measurements. This could be due to some part of the system already starting to behave as a three-phase system and providing an extra shear strength/stiffness to the sample (i.e., a heterogeneous distribution of hydrate morphology). The Biot two-phase model is unable to match ΔQ s −1 for S h > 12%, as this corresponds to a threephase system described by the inter-pore hydrate framework morphology (Sahoo et al., 2019). Encouragingly, Dai and Seol (2014) also found N = 5 gives a good fit to several experimental and field data sets, even though N varies with hydrate saturation and morphology.

Seismic and Sonic Frequency Measurements
We extend our analysis to sonic logging and seismic frequencies used in field studies. We apply the three-phase Biot model (GG model) to the 200 Hz laboratory resonant column V s and Q s −1 results for methane hydrate-bearing sand of Best et al. (2013) and Priest et al. (2009). Figure 3 shows the resonant column data with GG model predictions incorporating different values of the Tokyo permeability method N parameter. The model shows a reasonable match to both ΔQ s −1 and ΔV s for S h > 18% when N = 3.37, including the inferred peak in ΔQ s −1 in the range 32% < S h < 44%. As expected, the N value of 3.37 for hydrate formation in unconsolidated sand with 40% porosity is less than that in Berea sandstone (N = 5) with 22% porosity. ) changes measured during methane hydrate formation in Berea sandstone with modified Biot-type three-phase (a and b) and two-phase (c and d) models incorporating different permeability calculation methods indicated in the legend. PF denotes pore-floating hydrate morphology and PB denotes pore-bridging hydrate morphology. The Y-axis shows the difference in V s and Q s −1 compared to their respective values at zero hydrate saturation. N and M denote the parameters used for relating permeability to hydrate saturation in the Tokyo (Masuda, 1999) and Dai and Seol (2014) methods respectively as shown in the inset figure. The Tokyo method with N = 5 is used in (c and d).
We also apply the permeability based three-phase Biot model to borehole sonic logging data from Mallik 5L-38 gas hydrate drilling project (Figure 4). The V s and Q s −1 data were taken from . We see a good match with this data set for N = 3.2 to 3.5. Our aim here is to show that our model gives reasonable results that fall within the range of logging data. It is interesting to observe that the model predicts a very low attenuation for N < 3.1. When N is increased to N = 3.2, modeled Q s −1 increases, and if N is increased further then Q s −1 decreases and the peak in Q s −1 shifts to lower S h ; however, V s is not affected much with changes in N. The Mallik borehole data have a wide scatter around their main trends, however our model Q s −1 and V s values are within the range of measured data. ) of methane hydrate bearing sand measured at 200 Hz (resonant column laboratory data from Best et al., 2013;Priest et al., 2009) with rock physics model modified from  incorporating effective permeability. The uncertainty in velocity is <4.5% (Priest et al., 2009). The parameter N for permeability change with hydrate formation is for the Tokyo method (Masuda, 1999) (see Section 2.2). In (a) the purple line (N = 4) is hidden below the green line (N = 5). ) data (shown as gray dots) from the Mallik 5L-38 well in methane hydrate bearing sand with the three-phase rock physics model modified from  incorporating effective permeability. The parameter N relates to permeability changes with hydrate formation calculated using the Tokyo method (Masuda, 1999) (see Section 2.2).

Inversion of permeability From V s and Q s −1
We observe that V s and Q s −1 are highly dependent to changes in permeability. Hence, we performed a data inversion for permeability while keeping other model parameters constant for both seismic (200 Hz) and ultrasonic datasets. For inversion, we used the permeability ratio  E k (also called permeability reduction due to hydrate formation) which is defined as ratio of permeability of hydrate bearing sediment for a given hydrate saturation to the permeability of the fully water saturated sediment (intrinsic or absolute permeability of sample). In the fluid flow literature, relative permeability refers to permeability of a porous solid to one fluid phase with respect to other fluid phases in a multi-phase system (e.g., oil and water in sediments). In hydrate literature, we define a decrease in permeability due to hydrate formation in the pore of the host sediment. Authors have used the terms permeability ratio, permeability reduction and relative permeability. Here we use permeability ratio as we find it less ambiguous than the alternatives above.
We calculated velocity (V m ) and attenuation (Q m ) for a given model (two-phase or three-phase Biot's model) by varying the relative permeability. Then we found which permeability minimized the difference between experimental and modeled values. We use an objective function of We minimize the objective function (Equation 9) to find the best fit inverted permeability. In previous sections, we used Cycle four of the ultrasonic data; here we test our inversion on a different data set, the Cycle 3 data.
We applied the inversion to the 200 Hz data from Best et al., (2013) to find relative permeability. Inversion with two-phase Biot model shows a sensible result for the first ΔQ s −1 peak for S h < 20% and matches the Tokyo permeability with N = 3.37 ( Figure 5). For the second ΔQ s −1 peak S h > 20% with the two-phase Biot model, inverted r E k is above 1, which is not feasible and shows that there is a need for other Biot models. Hence, we applied the three-phase GG model for this inversion of the second ΔQ s −1 peak; we see a good match for Tokyo permeability N = 3.1-3.37 for the whole data set. The inversion with the three-phase model also gives good results for the first ΔQ s −1 peak S h < 20%; this shows that the three-phase model also works well for a two-phase system when saturation of one phase is low. Figure 5. Inversion of permeability ratio from shear wave velocity and attenuation of methane hydrate bearing sand measured at 200 Hz Priest et al., 2009) using (a) two-phase and (b) three-phase Biot's type poro elastic models. Permeability ratio is the ratio of permeability of hydrate bearing sediment to absolute permeability of the host medium with zero hydrate saturation. The inverted permeability is also compared with the Tokyo method (Masuda, 1999) with parameter N which accounts for permeability changes with hydrate formation (see Section 2.2).
When we perform the inversion of ultrasonic data with the three-phase GG model ( Figure 6) we see that inverted permeability matches for Tokyo permeability N = 3 for S h < 20% and N = 5 for S h > 20%. There is a sharp decrease in permeability near the second ΔQ s −1 peak around S h = 20%. This agrees with our hypothesis that the ΔQ s −1 peak at S h = 20% is due to the hydrate morphology transition from pore-bridging to inter-pore framework. We also see that the apparent decrease in shear wave attenuation ΔQ s −1 after the peak for S h > 20% is linked to a slight increase in permeability. This small increase in permeability can be due to a decrease in tortuosity. There is a decrease in measured resistivity for S h > 20% which also supports our hypothesis of a decrease in tortuosity for S h > 20%.
This slight decrease in tortuosity could be explained physically by the emergence of easier or less tortuous fluid flow pathways. Although we have to be cautious in interchanging electrical tortuosity with hydraulic tortuosity as their definitions differ (Ghanbarian et al., 2013), they follow similar trends (Katagiri et al., 2017). Keeping this in mind, it is interesting to note that the peak in ΔQ s −1 is at slightly lower S h compared to the peak in electrical resistivity (Figure 1). The Leclaire et al. (1994) model accounts for tortuosity via the geometrical aspect ratio parameter, that depends on the geometry or surface roughness of the boundaries between different phases (Section 3.2). As hydrate forms, the boundary between the pore fluid and the host sandstone porous framework remains almost constant, but the geometry of the boundary between the pore fluid and the hydrate changes significantly. The geometry of the hydrate/pore-fluid boundary depends on hydrate morphology, which in turn affects the geometrical aspect ratio parameter r hf . As we have several separate hydrate grains in the system, we can think of r hf as representing an effective bulk average value. As discussed in Section 3.2, r hf is 0.5 for spheroids, and r hf for ellipsoids decreases to zero as they become narrower (Lamb, 1945). It is well known that the hydrate surface is not smooth but has several elongated inclusions (e.g., Kuhs et al., 2004) that suggest a r hf close to 0.  found a good match for r hf = 0 with borehole log data from the gas hydrate Mallik 5L-38 well. Hence, it is logical to expect a hydrate surface with several inclusions, giving r hf close to 0.

Permeability From Electrical Resistivity
Electrical conduction in our sample is mainly via the movement of ions in the pore fluid. This is controlled by the pore fluid flow (i.e., permeability) and capacity of pore fluid to conduct an electric charge. The pore fluid permeability (affecting the number and tortuosity of ionic fluid conduction pathways) is dependent on the hydrate saturation and morphology (e.g., Spangenberg, 2001). Hence, changes in electrical resistivity can be linked to changes in relative permeability by quantifying the interplay among these phenomena.
We calculated the permeability ratio, or relative permeability, of the brine phase from electrical resistivity using the method of  and Archies' equations. The permeability ratio  E k is equal to relative permeability ( rw E k ) of water in the hydrate-bearing sediment (e.g., . Relative permeability ( rw E k ) can be related to resistivity index I (K. where, w E S is water saturation and wir E S is irreducible water saturation. Resistivity index is defined as the ratio of partially water saturated rock to that of the fully water saturated rock. We calculated the water saturation from changes in pressure and temperature during hydrate formation (Sahoo, Chi, et al., 2018). We use wir E S = 0.23 Figure 6. Inversion of permeability ratio, due to hydrate formation from ultrasonic shear wave velocity and attenuation data using the three-phase Biot's type poro elastic models (black dots). Permeability ratio is the ratio of permeability of hydrate bearing sediment to absolute permeability of the host medium with zero hydrate saturation. Permeability ratio was also obtained independently from electrical resistivity measurements (blue diamonds and green star). Pore fluid salinity increases during hydrate formation, included for the green star results, excluded for the blue diamond results (see Section 8).
The inverted permeability ratio is also compared with the Tokyo method (Masuda, 1999) with parameter N which accounts for permeability change with hydrate formation (see Section 2.2).
from Attias et al. (2008); they determined this experimentally for Berea sandstone with a similar porosity. Archie related the resistivity index to where t E is the measured sample resistivity, Φ is porosity, m is the cementation coefficient and  w E is brine resistivity. We used m = 3.1, used earlier for this same data set.
When hydrate forms, dissolved ions in the pore fluid are excluded from the hydrate structure (e.g., Ussler & Paull, 1995). Hence, salt concentration   E C in brine increases during hydrate formation.
where 0 w E S is the initial water saturation and 0 E C is the initial salt concentration. This increase in E C can decrease the resistivity of brine  . w E However, we see an increase in the resistivity with hydrate formation, suggesting the increase in ion concentration has less effect on overall resistivity. We still account for the effect of change in salt concentration on resistivity. The resistivity of brine depends on dissolved ion concentration, mobility, ionic strength, and temperature (e.g., Rusydi, 2018;Sorensen & Glass, 1987;Williams, 1996). While in our experiment Na and Cl are the dominant ions, in seawater several ions are present. So, ion concentration is generally defined as total dissolved solids (TDS). TDS in mg/L can be empirically related to the resistivity of pore fluids in milli-Ohm-meters (mΩ m) using equation below (e.g., Rusydi, 2018;Sorensen & Glass, 1987): where F is an empirical factor that varies between 0.1 to 0.15 for sea water (e.g., Rusydi, 2018;Sorensen & Glass, 1987). However, Equation 12 may not be always linear as it depends on several factors such as concentration, mobility, ionic strength, temperature, and water type. Still, for NaCl or sea-water within our experimental conditions Equation 13 still works (e.g., Rusydi, 2018;Sorensen & Glass, 1987). Rusydi (2018) found a value of F = 0.13 for brine and Hervas et al. (2006) found F = 0.12 for NaCl solution at 25°C. Hydrates are generally formed at lower temperatures than 25°C and the resistivity of brine is also temperature-dependent. So we accounted for this temperature change using the widely used expression (Barron & Ashton, 2005;ISO, 1985;Sorensen & Glass, 1987): where T (°C) is temperature,  E is the temperature compensation factor and  25 E (Ωm) is the resistivity of brine at 25°C. Although the change of resistivity with temperature is not completely linear, for our experimental condition the non-linearity is very low and Equation 14 is valid (e.g., Sorensen & Glass, 1987 We calculated the relative permeability of water from electrical resistivity using Equations 10-14 to see if it matches the relative permeability obtained from inversion of ΔV s and ΔQ s −1 data ( Figure 6). We show results with and without considering the effect of increase in salinity due to hydrate formation and find this effect is minor. The results in Figure 6 agree quite well for the limited available resistivity measurements, although there is a slight mismatch; these could be associated with uncertainties in resistivity measurements and in the values of various parameters (  wir , , E S m ) used in Equations 10-14. Moreover, these equations only consider the conduction through pore fluid and not any surface conduction, for example, due to presence of clay. Our Berea sample has a very low clay fraction of 0.02, and so these equations should be valid.

Implications
Our results show that the change in effective permeability from pore-bridging to interpore hydrate framework morphologies seen in the laboratory by Sahoo et al. (2019), if replicated in natural systems, can influence significantly the magnitudes of V s , and especially Q s −1 , at seismic frequencies. Patterns of V s and Q s −1 observed on field seismic surveys could provide valuable information on the permeability of hydrate-bearing geological formations through interpretation using the GG model with modifications for effective permeability.
Currently, measuring S-wave attenuation in situ is challenging due to factors like generating and detecting shear waves at the seafloor, heterogeneous overburdens giving different propagation paths for each receiver, geometric spreading, source and receiver coupling, and so on. Although S-wave attenuation data for natural hydrate deposits are so far very limited, some S-wave studies have been made at Hydrate Ridge (Petersen et al., 2007), offshore India (Riedel et al., 2014), Cascadian Margins (Dash & Spence, 2011;Riedel et al., 2014) and offshore Svalbard, Norway (Westbrook et al., 2008).

Conclusions
We measured changes in shear wave velocity (V s ) and attenuation (Q s −1 ) and electrical resistivity during hydrate formation in Berea sandstone. Here, Berea sandstone is considered to be an inert porous medium with a pore structure that is similar to that of many natural formations where hydrate is found. For hydrate saturation S h < 5%, there is no significant change in V s and Q s −1 . As more hydrate forms, the V S increases uniformly to S h = 20%, then decreases slightly. The Q s −1 shows two peaks at hydrate saturations of 6% and 20% that correspond to respective changes in the rate of increase of V s .
While existing rock physics models use arbitrary values of permeability, we modified Biot-type rock physics models for two-and three-phase systems to incorporate the effect of hydrate saturation and morphology on the permeability of hydrate-bearing sediments. We applied these models to V s and Q s −1 data collected at ultrasonic, sonic and seismic frequency ranges during hydrate formation in sand and sandstone, and found a good match. The modeling results suggest that the dominant cause for increasing shear wave attenuation is the decrease in effective permeability associated with hydrate formation.
We can explain the first Q s −1 peak by the hydrate morphology transition from pore-floating to pore-bridging, and the second attenuation peak by the transition from pore-bridging to inter-pore hydrate framework. The second attenuation peak is accompanied by a sharp increase in shear modulus of the hydrate-bearing sediments that agrees with the formation of inter-pore hydrate framework morphology according to our models. We suggest that the observed decrease in ΔQ s −1 for S h > 20% is caused by the complex interplay between changing hydrate grain shape (i.e., increasing sphericity, thus affecting tortuosity), and the inertial and viscous effects according to a three-phase Biot-type fluid flow mechanism. This increase in tortuosity is also supported by a decrease in measured electrical resistivity. We were able to invert our measurements of V s and Q s −1 during hydrate formation in Berea sandstone for effective permeability, for hydrate saturation >10%, using the Tokyo permeability method with N of 5. For hydrate formation in sand, the effective permeability can be obtained using the Tokyo method with N = 3.2-3.5.
We inverted V s and Q s −1 data for permeability by using both three-and two-phase Biot-type rock physics models and found a very good match for 648 kHz ultrasonic measurements, and 200 Hz resonant column data. To the best of our knowledge, this is first study providing clear evidence linking V s and Q s −1 to the permeability of methane hydrate-bearing sands and sandstones typical of natural gas reservoir formations. This raises the possibility of inverting remote geophysical V s and Q s −1 measurements for in situ hydrate reservoir permeability estimates, and for remotely monitoring reservoir permeability changes. This capability would benefit hydrate reservoir production strategies, and the possible use of natural gas hydrate formations for injecting and sequestering CO 2 as hydrate.