Laboratory observations of frequency-dependent ultrasonic P-wave velocity and attenuation during methane hydrate formation in Berea sandstone


 Knowledge of the effect of methane hydrate saturation and morphology on elastic wave attenuation could help reduce ambiguity in seafloor hydrate content estimates. These are needed for seafloor resource and geohazard assessment, as well as to improve predictions of greenhouse gas fluxes into the water column. At low hydrate saturations, measuring attenuation can be particularly useful as the seismic velocity of hydrate-bearing sediments is relatively insensitive to hydrate content. Here, we present laboratory ultrasonic (448–782 kHz) measurements of P-wave velocity and attenuation for successive cycles of methane hydrate formation (maximum hydrate saturation of 26 per cent) in Berea sandstone. We observed systematic and repeatable changes in the velocity and attenuation frequency spectra with hydrate saturation. Attenuation generally increases with hydrate saturation, and with measurement frequency at hydrate saturations below 6 per cent. For hydrate saturations greater than 6 per cent, attenuation decreases with frequency. The results support earlier experimental observations of frequency-dependent attenuation peaks at specific hydrate saturations. We used an effective medium rock-physics model which considers attenuation from gas bubble resonance, inertial fluid flow and squirt flow from both fluid inclusions in hydrate and different aspect ratio pores created during hydrate formation. Using this model, we linked the measured attenuation spectral changes to a decrease in coexisting methane gas bubble radius, and creation of different aspect ratio pores during hydrate formation.


I N T RO D U C T I O N
Gas hydrates are naturally occurring, clathrate compounds of gas (predominantly methane) found in marine and permafrost environments (Sloan & Koh 2007). Seafloor gas hydrates may have an important role in climate change (Archer et al. 2009), carbon dioxide sequestration (Jung et al. 2010) and continental slope stability (Sultan et al. 2004) and are being considered as a viable alternative energy resource (Boswell & Collett 2011). Methane hydrates formed in sand offer the most commercially attractive hydrate reservoirs (Boswell & Collett 2011).
Natural hydrates commonly exist in several morphologies within host sediments: (i) hydrate forming cement between mineral grains, known as cementing hydrate; (ii) disseminated hydrate growing freely in the pore space away from grain contacts, known as porefloating or pore-filling hydrate; (iii) hydrate contacting neighbouring mineral grains, known as pore-bridging or load bearing or frame supporting hydrate (e.g. Holland et al. 2008). When the saturation of the pore-floating hydrate is high enough to fill a pore or bridge neighbouring sediment grains, then it becomes pore-bridging (Priest et al. 2009;Waite et al. 2009). We will use the terms pore-floating and pore-bridging hereafter. Furthermore, Sahoo et al. (2018a) recently provided evidence for (iv) an 'interpore hydrate framework' morphology that is created when hydrate from adjacent pores coalesces to interlock the host sediments. The non-cementing morphologies (ii-iv) are thought to dominate natural hydrate systems, and (ii) & (iii) have been sampled and/or inferred at locations such as Mallik, Mackenzie Delta (Uchida et al. 2000), the Nankai Trough (e.g. Fujii et al. 2015), Alaminos Canyon, Gulf of Mexico (Boswell et al. 2009) and Mount Elbert, Alaska North Slope (Stern et al. 2011). Morphology (iv) also probably occurs in situ, but has only just recently been identified in the laboratory (Sahoo et al. 2018a).
Gas hydrate deposits are generally quantified locally from elastic wave velocity, attenuation and electrical resistivity data (Ecker et al. 2000;Lee & Collett 2006;Weitemeyer et al. 2006;Westbrook et al. 2008;Goswami et al. 2015). The accuracy of hydrate estimates depends on our understanding of the effect of hydrate on the elastic wave velocity and attenuation (e.g. Ecker et al. 2000;Best et al. 2013), and electrical resistivity (e.g. Spangenberg & Kulenkampff 2006), of the host sediments. The effect of hydrate is not only dependent on its saturation but also on its morphology. The cementing morphology has a much greater effect on the mechanical, elastic and electrical properties of hydrate-bearing sediments than the non-cementing morphologies (Ecker et al. 1998;Waite et al. 2004;Priest et al. 2009;Best et al. 2013). For example, sediments hosting cementing hydrate show a significant increase in elastic velocities even for hydrate saturations as low as about 5 per cent (e.g. Priest et al. 2005;Dai et al. 2012), whereas sediments hosting non-cementing hydrate (both pore-floating and pore-bridging) do not show such significant changes at low saturations (e.g. Priest et al. 2009).
Attenuation of P wave can also be used for indirect estimates of hydrate saturation (e.g. Guerin & Goldberg 2002;Chand & Minshull 2004;Priest et al. 2006;Westbrook et al. 2008;Best et al. 2013). In contrast to velocity, attenuation shows more complex variations with hydrate saturation for both cementing and non-cementing hydrates (e.g. Best et al. 2013;Marín-Moreno et al. 2017). While some natural hydrate bearing samples suggest that the presence of hydrate reduces attenuation (Dvorkin et al. 2003;Westbrook et al. 2008;Dewangan et al. 2014), others show an increase in attenuation due to hydrates (e.g. Matsushima 2006). Laboratory studies by Priest et al. (2006) and Best et al. (2013) showed significant variation of seismic (200 Hz) P-and S-wave attenuation with hydrate saturation, including distinct attenuation peaks. Best et al. (2013) suggested that these peaks could be explained by local viscous fluid flow related to the microporous structure of hydrate containing gas and water inclusions (Schicks et al. 2006). Recent studies have shown coexisting gas and hydrate within the gas hydrate stability conditions (e.g. Milkov et al. 2004;Lee & Collett 2006;Sahoo et al. 2018b), and that the presence of gas can itself have a strong effect on attenuation (Anderson & Hampton 1980a, 1980bMarín-Moreno et al. 2017).
In general, attenuation in fluid-saturated porous media depends on measurement frequency, owing to the viscous interaction of pore fluids and solids (Biot 1956a(Biot , 1956b. Partial and patchy saturation of gas and hydrates adds further complexity. While it is often difficult to relate unambiguously attenuation measured from remote seismic data and sonic well logs to specific hydrate contents (e.g. because of spatial averaging effects; Lee & Collett 2006;Waite et al. 2009), they can be related with greater accuracy in controlled laboratory experiments, which could offer insights into attenuation mechanisms. It is clear that there is a need for understanding the complex effect of gas hydrate on attenuation of elastic waves in host sediments.
Here, we present results from laboratory experiments combining the highly accurate ultrasonic pulse-echo method for measuring velocity and attenuation with controlled methane hydrate formation in Berea sandstone. These results provide further insight into the different attenuation mechanisms that operate in hydrate bearing sediments, and how they are likely to influence different measurement frequency bands. This knowledge will be useful for inverting and interpreting in situ seismic measurements.

Laboratory experiment
We performed consecutive cycles of methane hydrate formation and dissociation in Berea sandstone using a stainless-steel high pressure cell (Fig. 1a). We measured P-wave velocity (V p ) and attenuation (Q p −1 ) using an ultrasonic (frequency 448-782 kHz) pulse-echo system and signal analysis based on the method outlined by Winkler & Plona (1982). The sample was sandwiched between two cylindrical Perspex buffer rods of the same diameter as the sample. The ultrasonic transducer transmits pulses through the buffer rod into the rock sample, and into the second buffer rod. The buffer rods act as delay lines, allowing the first reflected pulse from the proximal boundary of the sample to be separated from the driving voltage from the signal generator. The same transducer is used to receive the reflections from the Perspex-rock interfaces, from which the velocity and attenuation can be derived. To mitigate aliasing, the signal was sampled at a frequency of 2.5 GHz, much higher than the transducer's nominal frequency of 1 MHz.
The first step in the pulse-echo data analysis is to separate the proximal and distal reflections from the sample buffer rod interfaces using time-domain gating with Tukey windows (Fig. 1b). The waveforms are included in the Supporting Information (Fig. S1). These separated signals are then transformed to the frequency domain via fast Fourier transform (FFT). We then take a spectral ratio of these separated signals at each frequency. As time-domain gating is a convolution in the frequency domain, the results from timedomain gating are a function of both the gating envelope and the signal to be gated. Thus taking the spectral ratio of two different signals with the same gating envelope does not necessarily mitigate gating artefacts. In order to achieve satisfactory results, we used a Tukey window of a length sufficient to capture the majority of the proximal and distal reflections. These window lengths were kept constant throughout the entire experiment. Our method then diverges from the method of Winkler & 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 & Plonar (1982) approximate analytical solution which is strictly only applicable to elastic (zero loss) materials. The shortest time-domain gate was 9 ms for the reflection from the base of the sample, giving a frequency resolution of 111 kHz. Each spectrum contains four data points; while it is possible to artificially increase the frequency resolution by zero padding, this is essentially interpolation, and so we choose to show the actual raw data only. The finite size of ultrasonic transducers causes diffraction-like beam spreading. This effect causes phase and amplitude distortions which are a function of the transducer radius, frequency and path length (effective distance from the transducer). The effect of the phase distortion is an apparent frequency-dependent decrease in traveltime resulting in higher velocity. The effect of the amplitude distortion due to diffraction is an apparent frequency-dependent amplitude reduction, leading to higher values of attenuation than the true value. We correct for these artefacts using the method of Benson & Kiyohara (1974) and  Papadakis (1972). Using these methods, velocity can be measured to an accuracy of ±0.3 per cent and attenuation to ±0.2 dB cm −1 (Best 1992) for homogeneous samples.
We used a right circle cylindrical sample (4.97 cm diameter and 2.06 cm height) of Berea sandstone with a porosity of 22 per cent and a permeability of 448 mD as a stable, inert and well-characterized porous medium. The presence of microcracks, which close at high differential pressure (e.g. Prasad & Manghnani 1997), in Berea affects elastic wave velocity and attenuation. We ran four cycles of hydrate formation and dissociation, with a differential pressure of 10 MPa in Cycles 1 and 2, and 55 MPa in Cycles 3 and 4. In this study we only discuss the data from Cycles 3 and 4, as microcracks are closed at such high differential pressures (e.g. Prasad & Manghnani 1997), and hence any changes in attenuation can be attributed to hydrate formation and dissociation processes only; there is insignificant attenuation due to microcrack squirt flow at 55 MPa.
We followed the hydrate formation method of Waite et al. (2004). The sample was oven dried and placed in the ultrasonic rig. First we applied a vacuum to the sample to remove as much air as possible. Then, using a syringe pump, we injected brine to fill 83.5 per cent of the sample pore space, which allowed an excess water condition (Ellis 2008;Priest et al. 2009). We left the sample for 3 days so that the pore fluids could re-distribute and then injected methane gas. Finally, four cycles of hydrate formation and dissociation were performed by cooling and heating the system, in and out of the gas hydrate stability conditions. The confining pressure was always adjusted to maintain a constant differential pressure in a given cycle. The pore fluid pressure, sample temperature and ambient temperature were recorded throughout the experiment. These pressure and temperature measurements were used to calculate hydrate, gas and brine saturation continuously during hydrate formation and dissociation using a thermodynamic method (Sahoo et al. 2018b). The details of this experimental method are given in Sahoo et al. (2018a, b).

Rock-physics modelling
We interpret the P-wave velocity and attenuation using an effective medium rock-physics model, the Hydrate Bearing Effective Sediment (HBES) model (Marín-Moreno et al. 2017). The HBES model is based on the central idea that hydrate bearing sediment can be treated as an effective medium of sediment grains, solid hydrate and coexisting gas. The HBES model extends the model approach of Best et al. (2013) that considers that hydrate can have inclusions of gas and water, similar to clay assemblages in Leurer (1997) and Leurer & Brown (2008). The HBES model uses the approach of Ecker et al. (1998) and Helgerud et al. (1999) for different hydrate morphologies such as cementing, pore-floating and pore-bridging. When an elastic wave passes through a hydrate bearing sediment, the difference in elastic compliances of the porous host medium (e.g. sand grain framework) and the porous hydrate grains creates local fluid pressure gradients between the gas and water inclusions in hydrate and the sandstone frame pores, leading to viscous fluid flow or squirt flow (here denoted as sub-microsquirt flow) and associated wave energy loss (Fig. 2). The inclusions in hydrate can be occupied by both gas and water, which have different bulk modulus and viscosity leading to different attenuation. The attenuation due to gas and water inclusions are calculated independently and added based on their relative proportions. The formation of hydrate in the pore space would create pores of varying aspect ratio (denoted as type-2 pores). There are also viscous fluid flow losses due to these different aspect ratio pores that are created during hydrate formation between the hydrate grains and the sand frame pore walls (here denoted as microsquirt flow; Fig. 2). Gas bubble resonance effects are modelled according to Smeulders & van Dongen (1997). These attenuation mechanisms are embedded in the Biot-Stoll global fluid flow model (Biot 1956a(Biot , 1956b giving a frequency-dependent effective medium solution for P-wave velocity and attenuation in hydrate-bearing sediments and rocks, as a function of both hydrate saturation (S h ) and morphology.
To apply the model, we first develop a physical understanding of the evolution of the system with hydrate formation. Our system has methane gas bubbles and brine in the initial state (no hydrate). Here, we can expect gas bubbles to occupy the maximum available pore space, and the bubble size to reduce upon hydrate formation. The saturation of pore phases (gas, hydrate and brine) was calculated continuously from changes in pore pressure and temperature during hydrate formation (Sahoo et al. 2018b). This calculation showed that even at maximum hydrate saturation, gas was present in our system, hence we included the effect of gas bubble resonance in our model calculations. In a previous synchrotron imaging study of hydrate formation in sand using a similar hydrate formation method, we observed that initially gas bubbles would almost completely fill the pores, and as hydrate forms, the bubble size reduces (Sahoo et al. 2018a). The Berea sample pore size (longest dimension) varied from 11 to 73 μm, measured using synchrotron imaging at the Swiss Light Source (SLS), Switzerland (Sahoo et al. 2018b). So we choose to vary the bubble radius (r) from 30 to 9 μm. As more hydrate forms, the concentration of type-2 pores (C ϕ2 ) increases, and their aspect ratios (α φ2 ) decrease. This behaviour was observed by synchrotron imaging of methane hydrate formation in sand using an analogous hydrate formation method by Sahoo et al. (2018a). The C ϕ2 and α φ2 were chosen to fit both measured velocity and attenuation, and were based on the values used in Marín-Moreno et al. (2017). Marín-Moreno et al. (2017) showed that at ultrasonic frequencies, the squirt flow due to fluid inclusions in hydrate (sub-microsquirt flow) has a negligible effect. So we choose to use inclusion aspect ratio and concentration values of 1 × 10 −4 and 0.48 (Best et al. 2013) respectively for all hydrate saturations. The fluid in inclusions can be both gas and water, and we considered 50 per cent of the inclusions filled with gas and 50 per cent with water. For the Berea sandstone mineral bulk and shear moduli, we used the Voigt-Reuss-Hill average (e.g., Mavko et al. 2009) of the mineral moduli of the sample constituents measured using X-Ray diffraction (3.4 volume per cent k-feldspar and 1.7 volume per cent illite Han et al. 2015). These average moduli do not account for quartz cementation, so we followed the approach of Mavko et al. (1998) and increased the grain coordination number to match the initial measurement at zero hydrate saturation. A grain co-ordination number of n = 12.5 gave the best match (Supporting Information Fig. S2). The input parameters for the HBES model are given in Supporting Information Table S1.

Velocity and attenuation changes during hydrate formation
Hydrate formation leads to an increase in V p at 667 kHz (Fig. 3a). The gradient of this increase in changes at hydrate saturation (S h ) around 5 per cent and 20 per cent (Fig. 3a). Unlike V p , attenuation does not increase uniformly with hydrate saturation (Fig. 3). The attenuation shows peaks at S h of ∼5 per cent and ∼20 per cent that coincide with the gradient change for velocity (Fig. 3). The gradient changes in velocity were explained by changes in hydrate morphology using rock-physics modelling by Sahoo et al. (2018a). Hence, it is logical to relate the attenuation peaks also to changes in hydrate morphology. The first gradient change is due to a transition from pore floating morphology to pore bridging morphology. The second gradient change occurs when sufficient hydrate has grown to interlock the sand grains, that is, there is a transition from a porebridging to an interpore hydrate framework morphology. Similar attenuation peaks were observed at 200 Hz in methane hydrate bearing sand at S h of 13 per cent and 32 per cent (Best et al. 2013), despite the frequency dependence of attenuation mechanisms (e.g.

Marin-Moreno et al. 2017).
The relative change in magnitude of attenuation between zero S h (Q −1 p, Sh = 0 ) and maximum S h of 26.3 per cent (Q −1 p, Sh = max ) is while the relative change in magnitude of velocity between zero (V p, S h = 0 ) and maximum S h is (V p, Sh = max ) These calculations show that the relative change in magnitude of attenuation is higher than that of velocity; this result is consistent in both Cycles 3 and 4. This observation supports previous laboratory (e.g. Yun et al. 2005;Priest et al. 2009) and modelling studies (e.g. Ecker et al. 2000;Chand et al. 2006) which have shown that, for non-cementing hydrates, the P-wave velocity only increases significantly for higher S h (> 40 per cent). On the other hand, our results show a significant increase in attenuation even with low hydrate saturations, which agrees with other studies (e.g. Chand & Minshull 2004;Best et al. 2013). Hence, for accurate quantification of hydrate deposits with hydrate saturations below 40 per cent, using both attenuation and velocity could be more diagnostic than using velocity alone. However, measuring intrinsic attenuation of hydrate-bearing sediments in situ is often challenging due to several factors such as impedance mismatches in heterogeneous media or geometric spreading (e.g. Huang et al. 2009).

Attenuation frequency spectra during hydrate formation
For a given hydrate saturation, the Q p −1 spectra show broadly similar shapes and magnitudes during both hydrate formation Cycles 3 and 4 ( Fig. 4 and Supporting Information Fig. S4); this indicates systematic changes. Q p −1 is generally higher when hydrate is present, although the attenuation value at 562 kHz for an S h = 2.9 per cent in Cycle 3 (Fig. 4a) is lower than that for S h = 0 per cent. The Q p −1 spectra at S h = 0 show only a small change with frequencies above 550 kHz compared to when hydrate is present. The calculated low frequency cut-off for the pulse-echo system is around 400 kHz (e.g. Winkler & Plona 1982), so it is possible that the data at 448 kHz were affected by amplitude variations due to side-wall reflections (the calculated cut off frequency is approximate as it depends on a number of assumptions such as beam spreading angle). Alternately, the particular frequency dependence of gas bubble resonance effects may be playing a role here. Hydrate dissociation in Cycle 3 may lead to small changes in gas bubble size from that at the start of Cycle 3. This may result in small differences in attenuation between Cycles 3 and 4 for S h < 6 per cent (Figs 4a and b). As hydrate forms, the Q p −1 frequency spectrum also changes in shape. The Q p −1 spectrum can be classified broadly into two categories based on the slope. Category 1 is for S h below 6 per cent where Q p −1 increases with frequencies above 550 kHz for Cycle 3. In Cycle 4, the increase in Q p −1 with frequency is within the uncertainty of Q p −1 (Fig. 4 and Supporting Information Fig. S4).Category 2 is for S h > 6 per cent where Q p −1 decreases with frequency. These systematic changes in Q p −1 frequency spectra shape show that there are frequency-dependent effects within this bandwidth that change as hydrate and gas saturation change. As expected from the Kramers-Kronig relations of causality (Kronig 1926;Kramers 1927), for every change in shape of the Q p −1 frequency spectra, an associated change is observed in velocity (Supporting Information Fig. S3). However, the magnitude of the change in V p with frequency is much lower than that of Q p −1 (eqs 1 and 2). This phenomenon has been noted in other studies of velocity dispersion (e.g. Sams et al. 2002). We include the velocity data in the Supporting Information (Fig.  S3) for completeness and to demonstrate the stability of the results. We will now focus on the interpretation of the attenuation data.

Modelling
We use the HBES model (Marín-Moreno et al. 2017) to help interpret our experimental measurements. Model parameters were based on the physical understanding of the evolution of the system with hydrate formation as discussed in Section 2.2. As discussed in Section 3.2, we see broadly two categories of Q p −1 frequency spectra based on their slope.
Category 1 (Figs 5a and b): In general, we see that Q p −1 increases with frequency. The HBES model shows a good match for V p above about 550 kHz, and an excellent match for Q p −1 . The model result is quite sensitive to gas bubble radius, as shown for r = 20 and 22 μm. The aspect ratio (α φ2 ) and concentration (C ϕ2 ) of type-2 porosity are 0.14 and 0.01 respectively. Some of the hydrate may form close to the grain contacts, as seen in synchrotron imaging (Sahoo et al. 2018a). Therefore, based on velocity and attenuation matching, we attributed 0.84 of hydrate (0.84 × S h ) to the pore floating morphology, and 0.16 to the cementing morphology (0.16 × S h ). Category 2 (Figs 5c and d): Here, the Q p −1 magnitude has increased from that of Category 1, and Q p −1 decreases with frequency. This change in slope could be explained by a change in the viscous loss relaxation frequency. In Category 1, an attenuation peak at a higher frequency than our measurement band could explain why we see attenuation increasing with frequency. In category 2, an attenuation peak at a lower frequency than our measurement band could explain the observed attenuation decreasing with frequency. This behaviour is consistent with the rock-physics model; the attenuation peak's magnitude increases, the peak moves towards a lower frequency when the aspect ratio of type-2 pores decreases, and the peak's magnitude also increases when the concentration of type-2 pores increases (Marín-Moreno et al. 2017). We used a bubble radius of 9-17 μm, which is lower than that of Category 1. The C ϕ2 was increased to 0.025-0.027 and α φ2 was decreased to 0.08. With increasing S h , the amount of hydrate growing close to grain contacts was increased. Based on velocity and attenuation matching, we used 0.6 × S h for hydrate in the pore floating morphology, and 0.4 × S h in the cementing morphology. We show results for a range of HBES model parameters to indicate their sensitivity (Figs 5c and d). It is possible to obtain a good match to the measured data for both V p and Q p −1 using reasonable parameter values. We can see that as more microsquirt flow occurs (C ϕ2 increases), V p decreases and Q p −1 increases (Fig. 5), allowing an optimum match to be obtained for both V p and Q p −1 . If there is no microsquirt flow, the model calculates very low Q p −1 and high V p . This specific C ϕ2 would vary if other model parameters were changed. However, this duality between velocity and attenuation provides a useful way to constrain the choice of parameters to match the C ϕ2 to both V p and Q p −1 . It was also evident that bubble size affects the shape of the frequency spectra for both S h categories, and this behaviour may be used to better constrain possible bubble size (which might represent a dominant or average bubble size). This observation confirms that broad-band measurements can give us better control over the choice of model input parameters, like aspect ratio and concentration (e.g. Chand et al. 2006;Best et al. 2013;Marín-Moreno et al. 2017). Dotted lines consider attenuation due to gas bubble resonance. Dashed lines consider attenuation due to viscous fluid flow between different aspect ratio pores generated due to hydrate formation (microsquirt flow) and gas bubble resonance. All model runs include the attenuation due to fluid inclusions in hydrate. The accuracy of attenuation is ±0.2 dB cm -1 (Best 1992). The gas bubble radius is denoted by r and aspect ratio for microsquirt flow is denoted by α φ 2 . The inset diagrams show possible pore phase distributions: brown for sand grains, black for gas, white for hydrate and blue for water. Figure 7. HBES rock-physics model predictions of Q p −1 from seismic to ultrasonic frequencies using the parameters calibrated by the ultrasonic Berea sandstone data at hydrate saturation of 24 per cent for Cycle 3 with differential pressure of 55 MPa. The purple line is very close to the horizontal axis. Microsquirt flow is due to different aspect ratio pores created during hydrate formation. Sub-microsquirt flow is due to fluid inclusion in hydrate. The black dot represents the attenuation measured at 200 Hz at hydrate saturations of 13 and 32 per cent (Best et al. 2013).
To understand the most dominant loss mechanisms at different S h , we ran the model both with and without microsquirt flow at 667 kHz (Fig. 6). Although the gas bubble radius and aspect ratio of type-2 pores would vary with S h , we kept them constant in each model run. This allows a conceptual understanding of the system. The results show that, for low S h , only gas bubble resonance can explain the elevated attenuation in this observation frequency band. Still, even for low S h , some Q p −1 (at S h = 4.6 per cent) are higher than the model results without microsquirt flow. This result indicates that even at low S h , attenuation due to microsquirt flow may contribute to the overall attenuation. At high S h , the attenuation due to microsquirt flow has a dominant effect compared to that of gas bubble resonance. The conceptual diagrams in Fig. 6 show the expected distribution of the system at high and low S h . Our measurements are in the ultrasonic frequency range, well above the frequencies used in seismic exploration and sonic well Downloaded from https://academic.oup.com/gji/article-abstract/219/1/713/5533330 by Hartley Library user on 21 August 2019 logging. However, having constrained (or calibrated) the HBES model input parameters using the ultrasonic data, we can extend our model to make predictions of V p and Q p −1 behaviour at seismic and sonic frequencies (Fig. 7). We chose the best fit Category 2 parameters including the inclusion aspect ratio and concentration of Best et al. (2013). The model gives similar Q p −1 magnitudes to seismic frequency (200 Hz) measurements of Best et al. (2013), who give Q p −1 at both S h = 13 per cent and S h = 32 per cent of 0.1 ± 0.02; those data are for hydrate formation in sand under excess water condition (Fig. 7). Our model shows that fluid inclusions in hydrate (for our choice of parameters) do not have much effect on Q p −1 in the ultrasonic frequency band, while they are expected to have a strong effect at seismic or sonic frequencies. Our model results are very sensitive to the chosen gas bubble size, especially the Q p −1 frequency spectrum shape and peak frequency. Without gas bubble resonance the magnitude of Q p −1 decreases significantly (Fig. 7). However, even in the absence of gas bubble resonance (yellow line, Fig. 7) Q p −1 in the frequency range of 10 3 -10 4 Hz is higher than with gas bubble resonance (blue line, Fig. 7). The yellow line represents a model run without gas bubble resonance, even though gas inclusions may be present. This result suggests that, in our model different attenuation mechanisms are not necessarily independent of each other. Depending on frequency, aspect ratio of hydrate inclusions and gas bubble radius, the total attenuation calculated for the partially water saturated sample can be higher or lower than the water saturated sample. This is because the presence of gas (in hydrate inclusions and the pores) affects the effective bulk modulus and viscosity of the pore fluid, and hence affects attenuation caused by all three mechanisms: sub-microsquirt flow, microsquirt flow and bubble resonance.
We used model input parameter estimates based on matching both velocity and attenuation data, exploring values within the range used by Marín-Moreno et al. (2017) and Best et al. (2013). Another approach would be to use the model to invert for concentration and aspect ratio of type-2 porosity. Such an inversion would be needed to explain the changes in gradient and peaks observed in V p and Q p −1 for the whole range of hydrate saturations. While rock-physics models can help us to understand intrinsic attenuation mechanisms, field measurements of intrinsic attenuation are still a challenge. Geometric spreading and scattering due to impedance mismatches in spatially heterogeneous media can significantly affect field measurements (e.g. Huang et al. 2009). For example, acoustic energy from borehole sources can reflect off the borehole wall due to the high stiffness of hydrate bearing sediments, and may not propagate through the hydrate bearing sediment at all

C O N C L U S I O N S
We measured the ultrasonic frequency spectra of P-wave velocity and attenuation during laboratory methane hydrate formation, and observed distinct changes in the shape of the frequency spectra with hydrate formation. We interpret two separate trends for hydrate saturations below and above about 6 per cent; attenuation increases with frequency for hydrate saturations below 6 per cent and decreases with frequency for higher saturations. Our measurements also support earlier experimental observations of attenuation peaks at specific hydrate saturations that depend on measurement frequency. Using a theoretical rock-physics model, we were able to provide a plausible explanation for the measured P-wave attenuation and velocity changes through a combination of different hydrate formation effects; these include reduction of coexisting gas bubble size and creation of different aspect ratio pores with hydrate formation. Our results show that measurements at ultrasonic frequencies can be used to better constrain various parameters used in rock-physics models, thus providing a means to constrain in situ seismic and sonic data inversions. Also, coexisting gas with hydrate is likely to enhance seismic and sonic wave attenuation, and must be considered when using attenuation to constrain in situ hydrate saturation estimates from elastic wave data.

A C K N O W L E D G E M E N T S
We acknowledge funding from the UK Natural Environment Research Council (grant nos. NE/J020753/1 and NE/N016041/1). TAM was supported by a Royal Society Wolfson Research Merit award. Data used in this paper are available at the National Geoscience Data Centre, UK and can be searched using DOI of the manuscript.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Typical examples of P-wave signals measured during hydrate formation in Berea sandstone. Only some of the wave forms are shown. All the data are included in the Supporting Information and also available from the National Geoscience Data Centre, UK and can be searched using the DOI of the manuscript. S h denotes hydrate saturation. Figure S2. Comparison of results from the rock-physics model of Marín-Moreno et al. (2017) to measured V p for Cycle 3 of the hydrate formation in Berea sandstone at zero hydrate saturation. Figure S3. Variation of P-wave velocity with frequency during hydrate formation in Berea sandstone. The differential pressure is 55 MPa. (a and c) Cycle 3 and (b and d) Cycle 4. Category 1 (S h = 0 per cent) and Category 2 (S h < 6 per cent) are shown in (a) and (b) while Category 3 (S h > 6 per cent) is shown in (c) and (d). Saturations of hydrate S h and methane gas S g are indicated in per cent. The accuracy of velocity measurement is up to ± 0.3 per cent (Best 1992). Figure S4. Uncertainty in P-wave attenuation with frequency for various hydrate saturations during hydrate formation at a differential pressure of 55 MPa in Berea sandstone. (a) Cycle 3 and (b) Cycle 4. Saturations of hydrate S h and methane gas S g are indicated in per cent. Table S1. Values used in the rock-physics model.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.