## Summary

The presence of gas hydrate in marine sediments alters their physical properties. In some circumstances, gas hydrate may cement sediment grains together and dramatically increase the seismic *P*- and *S*-wave velocities of the composite medium. Hydrate may also form a load-bearing structure within the sediment microstructure, but with different seismic wave attenuation characteristics, changing the attenuation behaviour of the composite. Here we introduce an inversion algorithm based on effective medium modelling to infer hydrate saturations from velocity and attenuation measurements on hydrate-bearing sediments. The velocity increase is modelled as extra binding developed by gas hydrate that strengthens the sediment microstructure. The attenuation increase is modelled through a difference in fluid flow properties caused by different permeabilities in the sediment and hydrate microstructures. We relate velocity and attenuation increases in hydrate-bearing sediments to their hydrate content, using an effective medium inversion algorithm based on the self-consistent approximation (SCA), differential effective medium (DEM) theory, and Biot and squirt flow mechanisms of fluid flow. The inversion algorithm is able to convert observations in compressional and shear wave velocities and attenuations to hydrate saturation in the sediment pore space. We applied our algorithm to a data set from the Mallik 2L–38 well, Mackenzie delta, Canada, and to data from laboratory measurements on gas-rich and water-saturated sand samples. Predictions using our algorithm match the borehole data and water-saturated laboratory data if the proportion of hydrate contributing to the load-bearing structure increases with hydrate saturation. The predictions match the gas-rich laboratory data if that proportion decreases with hydrate saturation. We attribute this difference to differences in hydrate formation mechanisms between the two environments.

## Introduction

Gas hydrates are found in marine sediments and permafrost regions within the gas hydrate stability zone and where sufficient methane is present. In most places, their presence is inferred from the presence of a bottom simulating reflector (BSR) (Shipley 1979) and velocity increase in the sediments relative to the normal velocity (Singh 1993). Recent studies have shown that the attenuation of seismic waves also increases with hydrate saturation ranging from sonic frequencies (Pratt 2003) to log frequencies (Guerin & Goldberg 2002; Matsushima 2005; Murray 2005). Since no direct measurements are available in most locations, the compressional and/or shear wave velocity distribution in marine sediments is often used to derive quantitative estimates of gas hydrate in the pore space (e.g. Lee 1996; Ecker 2000; Jakobsen 2000; Chand 2004). Attenuations in hydrate-bearing sediments are generally higher than those of water-filled, normally compacted marine sediments, so attenuation measurements can also be used to quantify the amount of gas hydrate (Chand & Minshull 2004). It has been inferred that the velocity and attenuation ratios of *P* and *S* waves give a further constraint on hydrate saturations (Matsushima 2005). Such anomalies can be used to infer the concentration of hydrate if the relationship between velocity/attenuation and gas hydrate content is known.

In a recent study comparing different methods relating gas hydrate saturation and velocity/attenuation, it was found that the self-consistent approximation/differential effective medium theory with a hydrate-sediment bi-connected morphology agrees well with measured velocities of hydrate-bearing sediments (Chand 2004). The attenuation behaviour is well explained by a theory based on Biot and squirt flow mechanisms (Chand & Minshull 2004). We modified these theories to develop a new inversion algorithm for predicting gas hydrate saturation using measurements of velocity and attenuation. We tested our inversion algorithm on laboratory measurements of hydrate-bearing sand samples (Priest 2005; Waite 2004) and from the Mallik 2L-38 well (Collett 1999; Guerin & Goldberg 2002). Using the inversion algorithm and measured parameters, we are able to predict hydrate saturations comparable to those determined independently by non-seismic methods.

## Forward Model

Jakobsen (2000) relate the seismic properties of clay-rich hydrate-bearing sediments to their porosity, mineralogy, microstructure, clay-particle anisotropy and hydrate saturation, following approaches developed by Hornby (1994) and by Sheng (1990). The building blocks are considered to be transversely isotropic due to the alignment of clay platelets, and the theory is based on a combination of a self-consistent approximation (Willis 1977; hereafter called ‘the SCA’, though there are other self-consistent approximations in the literature), a differential effective medium theory (DEM) (Nishizawa 1982) and a method of smoothing for crystalline aggregates (Bonilla & Keller 1985). The SCA/DEM model structure consists of inclusions embedded in a homogenous matrix; the effective elastostatic parameters can then be expressed in terms of a tensor relating the applied stress to the average strain in each inclusion (Hudson 1991). The model starts with a self-consistent method in which the elastic tensor is determined approximately using a single inclusion in a host having the elastic properties of the yet to be determined effective medium (Hornby 1994). Individual components are then replaced with an effective material. Starting with such a bi-connected medium created with self-consistent approximation, which is valid only for porosities below 60 per cent, the stiffness values at higher porosities and at other sediment compositions are calculated by successive operations of removing an infinitesimal subvolume of the host material, and replacing it with another component, using differential effective medium theory (Nishizawa 1982).

Since the DEM preserves the connectivity of the phases (Sheng 1990), this procedure produces sediment that is bi-connected at all porosities. The material is considered as made up of blocks of a fully aligned composite medium, but arranged in different orientations. A limitation of Sheng's (1990) approach is that it will not take into account the weaker bonding and greater compliance likely to exist at the edges of the individual domains of preferential particle alignment (Bennett 1991). The effective stiffness is calculated using the Voigt-Reuss-Hill average (Hill 1952), followed by a ‘method of smoothing’ (Bonilla & Keller 1985), which can be used to calculate the effects of interactions between neighbouring grains. This averaging and smoothing also minimizes error due to effects at edges when the orientation of individual components is different. The SCA/DEM theory does not consider pressure effects other than through those resulting from changes in porosity.

Hydrate can be added as a part of the load-bearing framework (called here the ‘cementation model’), so that grains of sediment are replaced by grains of hydrate, or as an inclusion. Since the bi-connected material created using the SCA at certain porosity involves porosity filled with water, replacement of water later using the DEM will not create the connectivity between hydrate and the matrix implied by the cementation model. In such an approach, the substitution of water by hydrate results in an inclusion model in which hydrate behaves like other unconnected grains in the composite. Jakobsen (2000) assumed that complete connectivity between the matrix and hydrate remains if hydrate is included in the starting SCA medium. However, at 100 per cent hydrate saturation, the velocities predicted using a matrix-hydrate-starting SCA model are the same as those predicted using a matrix-water-starting SCA model. This result is similar to that obtained using some other cementation models (Ecker 1998), where the complete cementation occurs at very low hydrate saturation and velocity then increases gently with increasing hydrate saturation to reach a value similar to that of a three-phase Hertz-Mindlin-Hashin-Shtrikman (HMHS) model (Ecker 1998). Recent studies have shown that such models are only valid for a limited range of hydrate saturations (Waite 2004). Therefore, we used a different approach by creating an SCA model starting with hydrate as part of the matrix as cement rather than as a part of the pore fluid as inclusions. In this case, the overall strength of the SCA bi-connected medium increases significantly with increasing hydrate content. In this new approach to modelling cementation, the connectivity of hydrate with the sediment matrix varies with hydrate saturation. The approach is implemented by creating a bi-connectivity between the sediment and the hydrate using SCA.

We start with the creation of an SCA medium consisting of the matrix (normally clay) and the pore fluid (normally water) at sediment porosity, if the porosity is less than 40 per cent. For porosities above 40 per cent the bi-connected medium is created with 40 per cent porosity. In either case, the porosity is adjusted using DEM to create unconnected space for the remaining (non-matrix) mineral grains. Then these grains are added into the unconnected pore space left out for them using DEM. This will give the base values for the seismic properties (*V*_{pbase}, *V*_{sbase}, *Q*_{pbase} and *Q*_{sbase}) in the absence of hydrate. We used 40 per cent porosity as the limiting value, since Dvorkin & Nur (1996) and Krief (1990) have predicted that most sediments lose their connectivity beyond this porosity and because loose-sand samples often exhibit ∼40 per cent porosity at atmospheric pressure (Priest 2005). Other predictions, such as those using the SCA, assume connectivity up to 60 per cent porosity, as may occur where large amounts of clay are present.

Hydrate in nature may occur as part of the load-bearing matrix, like clay, and/or as an inclusion like other larger grains. Therefore, we model a portion of hydrate saturation as load-bearing cement and the remaining as pore-filling inclusions (Fig. 1). Our algorithm first creates an SCA medium of matrix and pore fluid at a porosity, which is reduced from the true value by an amount equal to the part of hydrate saturation that is assumed to act as bi-connected with the matrix. Then the porosity is readjusted to leave vacant spaces for other grains to be included. The bi-connected, load-bearing part of the hydrate is added by replacing an equal volume of the matrix by hydrate using DEM, thus keeping the connectivity of the constituents in the matrix. The inclusion part of the hydrate is added as a replacement for pore fluid (water or gas). The remaining grains are added into the unconnected pore space as described earlier. Using this approach, we can simulate partial cementation by hydrate at low saturations and incomplete cementation at 100 per cent saturation. We use standard values for the physical parameters of the constituent mineral grains and the pore fluid (Table 1).

The model predictions for sediment containing equal amounts of quartz and clay for the pure cementation and pure inclusion models are shown in Fig. 1. Here, we assume sediment grain connectivity is a function of porosity. The SCA model estimates changes with porosity up to 40 per cent porosity and thereafter the porosity is established using DEM. Adding hydrate increases the connectivity and results in cementation and reduction of the porosity. Since some or all of the hydrate is added as a part of the bi-connected medium, the cementation effect of hydrate is also scaled according to hydrate saturation. In this modified, partially cemented model the velocities increase gradually from zero hydrate saturation and reach the fully cemented value at 100 per cent hydrate saturation. The *P*- and *S*-wave velocities increase from 1.54 km s^{−1} and 0.336 km s^{−1} to 4.45 km s^{−1} and 2.6 km s^{−1}, respectively, with increasing hydrate saturation and decreasing porosity for the cementation model (Figs 1c and d). For the inclusion model, the velocity increases are 1.54–3.56 km s^{−1} and 0.336–2.03 km s^{−1} for *P* and *S* waves, respectively.

The attenuation of acoustic waves in a medium can be due to many different mechanisms. These include friction, scattering, viscous shear relaxation, macroscopic fluid flow and intercrack squirt flow. At seismic frequencies the frictional and scattering effects can be neglected since the strain amplitudes are very low and wavelengths are large. It has been suggested that the major factor influencing attenuation at seismic frequencies is fluid flow and the viscous effects related to it (Klimentos & McCann 1990). Hence, the influence of fluids on attenuation can be represented by a combined poro-elastic wave equation including both Biot (Biot 1956a, b) and squirt flow (BISQ) (Mavko & Nur 1975, 1979) mechanisms (Parra 2000). The model assumes that the Biot flow takes place in the direction of wave propagation while squirt flow takes place perpendicular to it. The BISQ model relates the dynamic poro-elastic behaviour of a saturated rock to traditional poro-elastic constants such as porosity, permeability, fluid compressibility and viscosity, and the characteristic squirt-flow length. Chand & Minshull (2004) used the BISQ model in a multistructural framework model to establish the attenuation properties of hydrate-bearing sediments.

To model both the wave velocity and attenuative properties of the hydrate-bearing sediments, we combined the modified SCA/DEM theory (Jakobsen 2000; Chand 2004) and the BISQ attenuation mechanisms of fluid flow for hydrate-bearing sediments (Chand & Minshull 2004). The model predicts that, *P* and *S* wave velocities and attenuations increase with hydrate saturation. Attenuation also increases with frequency until it reaches a peak, which is related to the hydrate permeability, and then decreases (Chand & Minshull 2004). The peak attenuation frequency depends on the competing effects of hydrate saturation. With increasing hydrate saturation, the matrix strengthens and the attenuation decreases. However, increasing hydrate saturation also increases attenuation by increasing the proportion of fluid flow in the hydrate matrix and between hydrate and the sediment grains (Chand & Minshull 2004). Observations suggest that the latter effect dominates. Hence, at low hydrate saturations the frequency of maximum attenuation is determined by the properties of the sediment microstructure, while at higher hydrate saturation, it is determined by the properties of the hydrate microstructure. Chand & Minshull (2004) suggested a model in which the hydrate microstructure has its own pore-fluid flow properties that are different from those of the sediment matrix.

## Inversion for Hydrate Saturation

A formal inversion approach should generate both an estimate of the hydrate saturation and uncertainties of that estimate. Input parameters include the sediment composition with depth, the physical properties of the hydrate, water and the grains constituting the sediment, and the porosity. If the sediment composition and porosity are available as functions of depth, the prediction of hydrate saturation will be more accurate. The largest uncertainty comes from the uncertainty in porosity values used. In a field situation, the porosity often has to be estimated from a nearby borehole or sediment core. Seabed porosity may be determined from cores, and the porosity variation with depth is assumed to match that of the borehole. Comparing the different data sets of Collett (1999), Guerin (1999) and Jarrard (1995), we found such estimates introduce a maximum uncertainty of 5–10 per cent in the porosity values. Another factor influencing the estimates is the mineralogy of the sediment. Since elastic moduli of minerals vary widely, errors in mineral composition can introduce another 5–10 per cent uncertainty. Uncertainties in parameters relating to characteristics of fluid flow, such as permeability, viscosity, etc. result in uncertainties in the attenuation parameters of the composite. Also, the physical properties of the hydrate depend on the gases involved in its formation: structure I hydrates are the most common, but structure II hydrates form in the presence of higher order hydrocarbons.

Our inversion algorithm takes account of both the above uncertainties in the constants used and the uncertainties in the measured values of velocity and attenuation (*Vp*, *Vs*, *Qp* and *Qs*, hereafter denoted ‘inversion parameters’). The calculated seismic properties are compared with experimentally determined values. As the hydrate saturation is adjusted to improve the fit between the calculated and observed values, the amount of hydrate cementing the grains is also adjusted. A further complication arises in the case of laboratory measurements on pure sand samples, since sand grains are not cemented in the absence of hydrate. In this case, to obtain the correct seismic properties in the absence of hydrate, we must replace part of the connected quartz modelled by the SCA with unconnected quartz as inclusions using DEM.

The inversion process is iterated until the normalized misfit between the calculated and observed values is minimized. The minimization is done using Brent's algorithm (Brent 1973). Once the misfits are normalized, we follow a joint minimization approach giving equal weight to each of the inversion parameters. We followed the approach used by McKenzie & O'Nions (1991), modified from Parker (1977), to incorporate the uncertainties in the inversion. The function we minimized is given as

where,*R*and

^{o}_{i}*R*are the observed and calculated inversion parameters obtained and N is the number of inversion parameters (up to 4). To account for the contribution of uncertainties in the fixed parameters (porosity, etc.) to the uncertainty in the computed seismic properties, the standard deviation, σ for each inversion parameters is estimated as follows. where

^{c}_{i}*M*is the number of fixed parameters contributing to the uncertainty and (

*C*

^{j}

_{0i}−

*C*

_{0i}) represents the perturbation in the

*i*th inversion parameter resulting from changing the

*j*th fixed parameter from its standard value (Table 1) to its limiting value, and ɛ

_{i}is the uncertainty in the measurement of the

*i*th inversion parameter. These perturbations vary little with hydrate saturation, so for computational efficiency they were computed at zero hydrate saturation. Each misfit is weighted using the observed parameter values, so that the minimization of each parameter is given equal importance. The uncertainty in hydrate saturation is estimated by perturbing the inversion parameters with their standard deviation and then calculating the perturbation induced in hydrate saturation. Thus, the inversion algorithm gives an estimate of hydrate saturation and its total uncertainty.

## Application to Gas-Rich Laboratory Measurements

To test and tune our inversion algorithm, we used two laboratory data sets for which the hydrate saturations were determined independently without the use of effective medium approximation. We compared the measurements from these experiments with the predictions from our inversion algorithm for different amounts of hydrate cementation. In both cases, the laboratory hydrate was made in a gas-rich environment where hydrate saturation is limited by the availability of water.

### (a) Gas Hydrate Resonant Column (GHRC)

Structure I methane hydrate-bearing sand samples were created in the laboratory by Priest (2005, 2006), by melting fine-grained ice particles in the presence of methane gas, a method pioneered by Stern (1996). The elastic moduli of these composite samples were measured using the Gas Hydrate Resonant Column (GHRC) apparatus developed at University of Southampton (Priest 2005, 2006). The shear modulus, and hence Vs, was determined from the resonant frequency of vibration of the sample and attached drive mechanism. The new instrument can carry out flexural excitation as well, allowing the estimation of the bulk modulus, and hence the compressional wave velocity from the flexural and torsional resonance frequencies. Though the GHRC can be used to measure both compressional and shear wave attenuation (Priest 2006), the measurements are not made on fully saturated samples. Since fluid motion is a primary factor controlling attenuation, we do not use the laboratory attenuation measurements in the inversion for hydrate saturation.

It is estimated that measurements using the present apparatus produce an error of 0.9 per cent on values estimated for torsional mode vibrations and 4.5–6.3 per cent for values estimated based on flexural mode vibrations, depending the resonant frequency of vibration (Priest 2005). We used these values in our inversion to estimate the uncertainty in our model predictions. Also, since the resonant column works on the principle of free vibration, the resonant frequencies for torsional mode will be different from the flexural mode (Priest 2005). Hence, not only were the measurements for compressional and shear waves made at different frequencies but the resonance frequencies change due to changes in the sample strength with the addition of gas hydrate. Since the measurements were based on flexure of the sample, the pore space gas would have caused little effect on the moduli measurements.

The samples used for the study consist of fine-grained Leighton Buzzard sand (average grain size 100 μm). Hydrate formation in the pore space was controlled by mixing known quantities of ice (180–250 μm) particles with sand (Priest 2005). The torsional and flexural measurements were made at low strains (<10^{−5}) to minimize the attenuation due to friction. The sample was essentially a dry sand-hydrate composite since the amount of water/ice in the mixture was used to control the amount hydrate formed. Measurements at different pressures and porosities show that velocities increase with increasing pressure and decreasing porosity. The velocity dependence on pressure for hydrate-bearing samples is weak, with shear velocity changing by 2–4 per cent between pressures of 500 and 1500 kPa. For samples with zero hydrate saturation, the shear wave velocity increases by almost 28 per cent for the same pressure increase (Priest 2005). The *P*- and *S*-wave velocities increase from 1.639 km s^{−1} and 0.292 km s^{−1} at 500 kPa (46.5 per cent porosity) to 1.675 km s^{−1} and 0.374 km s^{−1} at 1500 kPa (46.2 per cent porosity), respectively, for loose sand. For dense sand the increase is from 1.654 km s^{−1} and 0.339 km s^{−1} (41.6 per cent porosity) to 1.697 km s^{−1} and 0.434 km s^{−1} (41.3 per cent porosity), respectively. The velocity difference between the loose and dense sand at 500 kPa pressure can be entirely attributed to the difference in porosity. However, there is only a very small change of porosity with the application of pressure (2 × 10^{−4} per cent kPa^{−1}), and the increase of velocity with pressure is entirely due to an increase in grain contact.

Our inversion approach requires that the velocities of the samples without hydrate are modelled first to determine the base properties of the host sediment medium. If we use a quartz–air bi-connected starting model and then use DEM to reach the desired porosity, velocities are over-predicted, suggesting that the sand is only partly connected. To represent this effect in our model, we assume that only some of the quartz grains participate in the starting SCA quartz–air bi-connected medium. To match the velocities at zero hydrate content correctly, we use a quartz–air starting SCA medium at the sediment porosity but replace some of the connected quartz with quartz inclusions using DEM.

We can represent the velocity increase with pressure as due to higher connectivity between quartz grains. The increased connectivity can be modelled with more sand grains participating in the starting SCA bi-connected medium. The observed increase of velocity with increasing pressure can be modelled with 2 per cent quartz participating in the starting SCA bi-connected medium at 500 kPa and 4 per cent quartz participating in the starting SCA bi-connected medium at 2000 kPa. To model the acoustic velocity dependence on pressure as purely a porosity effect would require an unrealistic 16 per cent porosity change, from 40 to 34 per cent, to produce the velocities observed at 1500 kPa. Measured velocities in hydrate-bearing sand samples are shown in Fig. 2. Velocity is observed to increase steadily from zero hydrate saturation, with a disturbed zone at 3–5 per cent hydrate saturation where it shows a chaotic behaviour. This chaotic behaviour may be partly attributed to porosity changes resulting from a difference in packing density (Fig. 2), but could also result from a heterogeneous distribution of grain contacts (Priest 2005). The sharp increase of velocities for saturations up to ∼5 per cent may be attributed to the fact that the hydrate forms initially at the grain contacts where the water is concentrated, thereby increasing the connectivity and cementation between the grains. This hypothesis is supported by the reduction in the pressure dependency of velocity (fig 5 of Priest 2005) and the reduction in grain contact compliance (Clayton 2005), as hydrate saturation increases above 5 per cent. At higher water/ice saturations, water may remain at the initial ice locations partially wetting the grains and partially occupying the pore space while leaving some pore space totally dry.

Hence, at hydrate saturations of less than ∼8 per cent, hydrate is considered as part of the frame (cement). However, at higher saturations, if we assume that all of the hydrate is cementing the grains together, the hydrate saturations predicted by our inversion algorithm are lower than those observed (Fig. 3). To match the observations, the proportion of hydrate cementing the grains decreases to 85, 55 and 38 per cent at 10, 20 and 35 per cent hydrate saturations, respectively (Fig. 3). Conversely, if we assume that hydrate forms as pore space inclusions, the amount of hydrate present in the pore space is overestimated (Fig. 3). Hence a microstructure in which hydrate partially cements the grains and partially forms a framework within the pore space explains the observed phenomena.

### (b) Gas Hydrate And Sediment Test Laboratory Instrument (GHASTLI)

Physical properties measured with the US Geological Survey's GHASTLI facility are used here to study unconsolidated, partially water-saturated Ottawa sand samples containing an interconnected structure I methane gas hydrate phase. A description of GHASTLI and a detailed description of the sample construction and experimental methods used in this work are given by Winters (2004). The experiment was conducted using 1 MHz compressional wave piezo-electric transducers. Details of the measurements are given by Waite (2004). Here we model the *P*-wave velocities for three samples, GH083, GH084, GH085, to infer the hydrate microstructure. The samples were quartz sand with an average grain size of 370 μm (Table 2).

The partially saturated Ottawa sand pack (GH084) without hydrate did not pass any compressional wave, but measurement on a water saturated sample (33.8 per cent porosity) gave a compressional wave velocity of 1.88 km s^{−1}. Similar measurements by Winters (2004) on Ottawa sands with porosities 30.4 and 33.5 per cent gave velocities of 1.89 km s^{−1} and 1.861 km s^{−1}, respectively. The net change in porosity due to the application of pressure was negligible (Waite 2004). In our modelling, we assume that the gas phase remains connected within the pore space, because the measurements were made in a gas-rich environment and the pore pressure was maintained at a constant value of 12 MPa.

Winters (2004) and Waite (2004) modelled the observed velocities using theories proposed by Ecker (1998). The rock-physics models proposed by Ecker (1998) and Helgerud (1999) are based on theories of Dvorkin & Nur (1996) for different modes of occupancy of secondary material in the pore space of a random pack of spherical grains. The models describe extreme scenarios where hydrate forms only at grain contacts, only by enveloping grains, or only in the pore space. Another accompanying model uses the HMHS averaging model to describe the average properties of a partially structured composite. In HMHS model, the Hashin-Shtrikman bounds are used to extend the properties of the composite away from the totally bi-connected medium created using Hertz-Mindlin theory at 36 per cent porosity. Winters (2004) proposed that hydrate-bearing natural sediment properties are best described using the HMHS model, but frozen samples and hydrate formed in the laboratory show evidence of cementation. Measurements and modelling by Waite (2004) showed the grain-enveloping cementation model agreed with measured results within 3 per cent, whereas the grain-contact cementation model over predicted the observations by >20 per cent, and the pore space model under predicted the observations by 60 per cent.

We modelled the measurements of *P*-wave velocity of these hydrate-bearing sand samples containing residual gas using the SCA/DEM cemented approach (Fig. 4). There is ±0.5 per cent uncertainty in the velocity measurements, resulting from a combination of uncertainties in sample length and measured time of flight for the acoustic signal through the sample. The gas is considered as a continuous phase, and all initial water is assumed to have been converted into hydrate. The physical properties of the gas are estimated using the expressions of Mavko (1998) using the sample temperature of 279 K and pore pressure of 12 MPa used during the measurements. The details of the physical parameters are given in Table 1. The properties of the hydrate-bearing samples are given in Table 2. Again, we compared *P*-wave velocities from our forward model with observed velocities, and our inferred hydrate saturations with the known experimental values, to estimate the amount of cementation occurring in gas-rich systems. In the GHRC experiments, the percentage of hydrate participating in cementation appears to decrease with hydrate saturation: 100, 60 and 48 per cent for 20, 40 and 70 per cent hydrate saturations, respectively (Fig. 5). Relative to the GHRC experiments a greater fraction of hydrate was cement in the GHASTLI experiments. This difference likely results from the larger grain size used in GHASTLI experiment (250 to 500 μm) than in the GHRC experiment (100 μm), which decreases the number and increases the size of the pore spaces.

## Application to Data Set from Mallik 2L-38 Well

We tested our inversion algorithm by predicting hydrate saturations using various combinations of input inversion parameters from the Mallik 2L-38 well. We chose these data because, as shown in Fig. 6, all the inversion parameters (*P*- and *S*-wave velocities and their attenuations) are available, and there is reliable information regarding the sediment composition and porosity (Collett 1999; Guerin & Goldberg 2002). Our results (Figs 7a–c) may be compared with hydrate saturations predicted from resistivity data (Collett 1999). We modelled the hydrate saturation in three different ways, using the inversion parameters individually and in combination to predict the hydrate saturation. In the first model, cementation increases linearly with increasing hydrate saturation until all components became totally interconnected at 100 per cent hydrate saturation (Fig. 7a). In the second model we assumed that all of the hydrate acts as cement (Fig. 7b). The predictions of this second model agree well with those inferred from resistivity at high hydrate saturations, but poorly at lower saturations. In the third model, all of the hydrate acts as inclusions. The predictions of this model (Fig. 7c) are well above the resistivity-derived values. In contrast to our results from the laboratory data sets, we conclude that in nature, the amount of hydrate cementing the grains increases with hydrate saturation.

There is a large misfit between observed and predicted attenuations (Fig. 6c); this discrepancy may be due to limitations in our attenuation model, or alternatively, to variations in unconstrained parameters such as the permeability. Here we used only a single permeability value to model the whole sediment column (Table 1), since there are no available data to justify the use of variable permeability. The use of more than one inversion parameter improves the results considerably, though caution is required in the use of attenuations as inversion parameters since there are more unknowns involved. However, it is clear that attenuation can be used as a qualitative guide to indicate the presence of hydrate. The results shown here differ from those of Chand (2004) because that study used the approach of Jakobsen (2000), in which cementation occurs only through replacement of water in the pore space and does not affect the matrix.

## Application to Water-Saturated Laboratory Measurements

In a recent study, Yun (2005) measured *P*- and *S*-wave velocities on structure II hydrate-bearing sand samples (100 μm) created using tetrohydrofuran (THF) in solution. They observed low velocities that were well matched by a non-cementing effective medium model for all but the highest saturations. We analysed their observations using our SCA/DEM algorithm and compared them with results from the Mallik well and laboratory analyses (Fig. 8). From our models using log data from the Mallik 2L-38 well and laboratory samples, we infer greater cementation than Yun 's (2005) infer from their observations. The maximum amount of cementation we infer from their data is ∼5 per cent. Their result may be due to the formation of hydrate from dissolved THF in uniform and large pore sizes, which will result in the formation of hydrate entirely at the pore centre. In nature, the pore sizes will vary and, as we inferred from the experiments, there will be greater cementation when there is a larger variation in pore and grain sizes.

## Discussion and Conclusions

We have developed an algorithm for predicting the physical properties of sediment containing hydrates, and have compared our calculated results with well-log measurements and also with laboratory-formed sand samples containing hydrate. Acoustic velocities increase dramatically in sediment saturated with hydrate. Our model treats only hydrate in the pore space; we do not consider the case where hydrate fills cracks and fissures, as sometimes observed in nature, and this case merits further research. The attenuation also changes dramatically in the case of well log measurements, but only to a limited extent in the laboratory samples since they are essentially dry. We modelled the hydrate saturation using theories describing behaviour of composites and fluid flow when excited by a seismic wave. To match laboratory data from gas-rich samples, we used a model in which hydrate starts cementing from the start of hydrate formation in the pore space, and the cementation decreases with increasing hydrate saturation. However, using observations from the Mallik 2L-38 well, we infer that the cementation increases linearly with hydrate saturation. We infer a similar linear increase of cementation, though smaller, with hydrate saturation for water-saturated samples of structure II hydrate (Yun 2005). We attribute this difference in cementation to the formation environment. In the gas-rich laboratory measurements, the hydrate formation was controlled by locations of residual water. At low water saturations (for creation of low hydrate saturations) the water seems to concentrate along the grain boundaries due to surface tension rather than in the centre of the pore spaces. In such cases, the hydrate is initially entirely cementing, but later forms as inclusions in centre of the pore space. In water-rich environments, hydrate formation may be more uniformly.

If a linear increase of cementation with saturation is widespread in nature, our algorithm may be generally applicable. In our model, the hydrate is porous and results in increase of attenuation of the composite material. Our approach is further justified by recent observations of cementation and microporosity in both laboratory and natural samples using SEM imaging (Stern 2004, 2005). These images show that hydrate starts forming as meso-porous substance, but with the passage of time it becomes more crystalline and microporous, with grain sizes of the order of several tens of microns. This fine scale structure of hydrate may explain its formation as a structure even in sediments containing large amounts of clay. A similar study of both laboratory and reconstituted natural samples showed a large increase in velocities for both sand samples and samples containing clay (Winters 2005).

From our study we reached the following conclusions:

An observed increase in velocity with pressure for hydrate-free laboratory sands can be explained through a combination of decreasing porosity and increasing connectivity between grains; the latter effect may be modelled through appropriate adjustments to the SCA/DEM approach.

The observed increase in velocity due to formation of hydrate in pore space can be explained only through hydrate forming as a partially cemented phase that increases the total connectivity and thus strengthens the composite.

Our analysis of laboratory and field data sets suggests that the hydrate microstructure is mainly determined by the location of the water in the pore space and by the distribution of pore and grain sizes in the sediment matrix.

### Acknowledgments

This work was funded by European Commission under the contract EVK3-CT-2000-00043 (HYDRATECH). JAP and AIB were partly funded by Natural Environment Research Council. We thank G. K. Westbrook for many constructive discussions during the course of this work, and M. W. Lee, D. Twichell and two anonymous reviewers for constructive comments.

## References

*Microstructure of Fine-Grained Sediments: From Mud to Shale*

*Algorithm for Minimization without derivatives*

*Practical handbook of physical properties of rocks*

*Scientific Results from JAPEX/JNOC/GSC Mallik 2L-38 Gas Hydrate Research Well, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada Bulletin 544*

*in situ*elastic properties of gas hydrate-bearing sediments on the Blake Ridge

*Proceedings of the Ocean Drilling Program, Scientific Results 146*

*Scientific Results from JAPEX/JNOC/GSC Mallik 2L-38 Gas Hydrate Research Well, Mackenzie Delta, Northwest Territories, Canada, Geological Survey of Canada Bulletin 544*

*P*and

*S*waves (full waveform sonic)

*The Rock Physics Handbook—Tools for Seismic Analysis in Porous Media*

*Proc. Vth Int. Conf. Gas Hydrates*

*Proceedings of the Vth International Conference on gas hydrates*

*Proceedings of the Vth International Conference on Gas Hydrates*