Gas hydrates are a potential energy resource and hazard for drilling and infrastructure, yet estimates of global volume vary by over three orders of magnitude. Hydrates are electrically resistive compared to water saturated sediment and so electromagnetic methods provide an additional tool to seismic surveys and drilling for determining hydrate saturations. A marine electromagnetic survey was carried out at Hydrate Ridge, Oregon, USA, with the aim of testing the use of controlled source electromagnetic (CSEM) and magnetotelluric (MT) methods to map gas hydrate and free gas below the gas hydrate stability zone. A 2-D CSEM inversion supports the scenario deduced from previous seismic and drilling results, which indicate two mechanisms of hydrate emplacement: a transport-dominated and reaction-dominated regime. A prominent resistive region of 2.5–4 Ωm at a depth of about 130 mbsf, near the seismic bottom simulating reflector (BSR), suggests that 27 to 46 per cent of the bulk volume is filled with hydrate, depending on whether Archie′s Law or the Hashin-Strikman bounds are used. This is representative of a reaction-dominated regime for hydrate emplacement, and where a significant low velocity zone exists based on a seismic tomography inversion, suggests large quantities of free gas below the BSR. Electrical resistivity logging while drilling (LWD) data show general agreement with the CSEM inversion model except for a CSEM-derived resistive region at seismic horizon A, known to transport free gas into the gas hydrate stability zone. Inversion of MT data collected simultaneously during the CSEM survey provides a complimentary low-resolution image of the shallow sediments and shows folding in the accretionary complex sediments similar to that imaged by a tomographic seismic velocity model.
Hydrate Ridge is located on the accretionary complex of the Cascadia subduction zone where the Juan de Fuca Plate subducts obliquely (N69°E) beneath the North American Plate at a rate of 42 mmyr−1 (MacKay et al. 1992; DeMets et al. 1990) (Fig. 1A). The subducting plate′s thick sediment cover, 3–4 km of sandy and silty turbidites, is accreted to North America by offscraping at the deformation front or by underplating beneath the accretionary complex tens of kilometres east of the deformation front, creating an extensive fold and thrust belt on the continental slope (Tréhu et al. 2006a). One of the resulting ridges is a 25 km long by 15 km wide north–south trending feature called Hydrate Ridge, located approximately 80 km offshore Newport, Oregon, USA. Hydrate Ridge is located where the dominant direction of thrusting at the deformation front undergoes a transition from landward vergence to the north to seaward vergence to the south (MacKay et al. 1992). As the underthrust sediments dewater, there is upward migration of fluids, which is likely responsible for the occurrence of gas hydrate here and at other active accretionary margins (Tréhu et al. 1999; Peacock 1990).
Natural gas hydrate, a type of clathrate, is an ice-like solid that consists of a gas molecule, typically methane, encaged by a water lattice (Sloan 1990). Methane hydrates are found worldwide in marine and permafrost regions where the correct thermobaric conditions exist and sufficient water and gas molecules are available (Sloan 1990; Kvenvolden 2003). The quantity and distribution of gas hydrate in sediments is important because of its potential as an energy resource (Moridis & Sloan 2007) and as a trigger for slope instability (Mienert et al. 2005; Nixon & Grozic 2007; Paull et al. 2007; Sultan et al. 2004; Smith et al. 2004; Field & Barber 1993), which may threaten seafloor infrastructure (Kvenvolden 2000; Hovland & Gudmestad 2001). As more deep drilling and production operations are carried out within the thermodynamic stability conditions for hydrate (Dawe & Thomas 2007) the consequences of drilling into hydrate sediments will become a bigger threat, since drilling and production fluids can cause hydrate to dissociate and cause wells to blow out (Ostergaard et al. 2000).
Seismic data alone are often insufficient for accurately resolving the amount of gas hydrate in sediments. One seismic signature often associated with gas hydrate occurrence is a bottom simulating reflector (BSR), which typically marks the phase change of solid hydrate above and free gas below the BSR (Shipley et al. 1979). However, the BSR may not indicate the existence of hydrate, as was observed on DSDP Leg 84 site 496 and site 596 (Sloan 1990, p. 424; Sloan & Koh 2007, p. 575). In fact, it requires very little gas to form a strong seismic reflector (Domenico 1977). Other types of seismic signatures have been noted at Blake Ridge by Hornback et al. (2003) and Gorman et al. (2002), such as a fossil BSR, seismic blanking and seismic bright spots. While seismic methods are often able to detect the lower stratigraphic bound of hydrate, the diffuse upper bound is not well imaged and there is often no seismic reflectivity signature from within the hydrate region.
Hydrate is electrically resistive compared to the surrounding water saturated sediments (Collett & Ladd 2000), which provides a target for marine electromagnetic (EM) methods. Marine EM methods can be used to image the bulk resistivity structure of the subsurface and are able to augment seismic data to provide valuable information about gas hydrate distribution in the marine environment (Edwards 1997; Yuan & Edwards 2000).
At Hydrate Ridge evidence for gas hydrate comes from the BSR present over much of the area (Tréhu et al. 1999), recovered samples of massive hydrate (Bohrmann et al. 1998), as well as logs and cores from Ocean Drilling Program (ODP) Leg 204 (Tréhu et al. 2006a). Hydrate emplacement is lithologically controlled and due to two main mechanisms: a focused high-flux regime and a distributed low-flux regime (Tréhu et al. 2006c). The most representative example of the focused high-flux regime is the gas-charged seismic horizon A, which is a conduit transporting thermogenic or altered biogenic gas from great depth within the accretionary prism to the summit of Hydrate Ridge (Tréhu et al. 2004a; Claypool et al. 2006). Hydrate formation due to the presence of methane within the sediments from in situ microbial methane production leads to diffuse fluid flow and dispersed hydrate throughout the sediment, and is evidenced by the pervasive BSR over much of the Cascadia accretionary complex (Tréhu et al. 1999; Claypool et al. 2006). No gas hydrate is present in the upper 30 metres below seafloor (mbsf) because the methane content of the pore water is below saturation (Tréhu et al. 2006a), except in regions of vigorous focused fluid flow. The ODP Leg 204 drilling and 3-D seismic data provide information about the distribution of hydrate, which may be used to assess electrical resistivity models derived from the controlled source electromagnetic (CSEM) and magnetotelluric (MT) data collected in the study presented here.
Application of csem and mt techniques
The application of marine CSEM methods to hydrate detection was first considered by Edwards (1997). He modelled the transient electric dipole–dipole method as a means of estimating hydrate volume and argued for the usefulness of EM methods in augmenting drilling and seismic techniques. Field studies conducted at the Northern Cascadia margin off the west coast of British Columbia, Canada, demonstrated the merits of CSEM by showing the existence of hydrate when no BSR is present (Yuan & Edwards 2000), and the existence of hydrate or free gas in seismic blanking zones thought to represent hydrate-bearing pipes (Schwalenberg et al. 2005). The number of EM surveys to image gas hydrate is slowly increasing: see (Goto et al. 2008; Dunbar 2008; Evans 2007; Ellis et al. 2008; Schwalenberg et al. 2009; Zach & Brauti 2009; Darnet et al. 2007).
The CSEM method used in this paper is a frequency-domain technique whereby a horizontal electric dipole transmitter is towed on or close to the seafloor and receivers record the transmitted electric and magnetic fields at various frequencies and ranges (Constable & Cox 1996). The method has recently been used commercially for hydrocarbon exploration (e.g. Ellingsrud et al. 2002; Eidesmo et al. 2002; Hesthammer & Boulaenko 2005; Constable & Srnka 2007). The purpose of this study was to see if the equipment and method developed to image many kilometres into the crust could be adapted to image a shallow gas hydrate reservoir, in the top hundreds of metres of the seafloor. Fig. 2 shows the general CSEM layout, whereby autonomous seafloor receivers spaced a few metres or kilometres apart are typically arranged in lines or grids. A man-made EM source field is generated by a deep-towed transmitter that is towed in a pattern around the seafloor receivers. Electric fields recorded by receivers are larger over resistive seafloor structures such as basalt, salt, carbonates, hydrocarbon reservoirs or gas hydrates (Constable 2006). The same seafloor receivers can also be used to collected MT data.
Prior to the advent of a broad-band marine MT instrument (Constable et al. 1998), only CSEM techniques could be used to study the oceanic crust, as MT studies were limited to shortest periods of several hundred seconds at best (Filloux 1987), which images only mantle structure. [There are a few exceptions. As a result of energetic source fields at high latitudes, Heinson et al. (2000) obtained MT responses to 40 s periods in water depths around 3000 m on the Reykjanes Ridge, and Jegen & Edwards (1998) used vertical gradient sounding (VGS) in 2600–2700 m water at the Juan de Fuca Ridge to obtain 10 s data. The VGS method uses magnetometer data only and makes some assumptions about 1-D structure to generate a MT response; it is probably better considered a magnetic transfer function (e.g. Key & Constable 2011).] By extending the high frequency limit to about 0.1 Hz, the broad-band MT instrument allows imaging at crustal depths, as shallow as a few hundred metres and overlapping with CSEM measurements to some extent. The CSEM method is still more sensitive to shallow and resistive parts of the seafloor, while MT sounding is sensitive to deeper and relatively conductive features.
The survey at Hydrate Ridge provided an opportunity for the simultaneous collection of MT and CSEM data. Although the focus of the experiment was to use the CSEM data to image gas hydrates, MT data provides information about the large scale regional features of the accretionary complex sediments at Hydrate Ridge. In this paper, we expand on the details of the modelling, collection and analysis of the CSEM data presented previously in Weitemeyer et al. (2006b) and Weitemeyer et al. (2010). A geologic interpretation is presented by combining other data sets available such as resistivity well logs from ODP Leg 204, seismic stratigraphy from Chevallier et al. (2006), and seismic velocity from a tomographic inversion of first arrival times recorded on an array of ocean bottom seismometers (Arsenault et al. 2001). We present the first formal treatment of the MT data analysis and inversion. MT data is almost always recorded by the seafloor receivers in CSEM surveys, even though it is less often processed and interpreted. Comparisons are made between MT and CSEM inversions.
1-D CSEM experimental design study
Prior to the field work we carried out an experimental design study using 1D models based on ODP Leg 204 resistivity logs. Well logs showed that the background resistivity at Hydrate Ridge is ≈1 Ωm; and that sediments containing hydrate vary between 2 and 6 Ωm; 2 Ωm was conservatively used as the resistivity for hydrate-containing sediment in the model studies. According to ODP Leg 204 Initial Reports, the hydrate distribution north of the southern summit begins at 45 mbsf, which was used as the starting depth for the hydrate layer in the model studies. The thickness of the hydrate layer comes from the seismic BSR depth, about 150 mbsf. An average water depth of 910 m was used. The transmitter tow altitude is 90 m for the 1-D model studies.
There are two geometric modes in CSEM data, to excite purely radial modes the transmitter is towed directly over the top of the receivers, and to excite azimuthal modes the transmitter is towed at a distance perpendicular to the receiver. Using the 1-D CSEM codes of Flosadóttir & Constable (1996) and Key (2009) it was found that the azimuthal CSEM mode has very little sensitivity to a hydrate layer and that the radial mode is most sensitive to such a layer. A suite of frequencies (0.1–300 Hz) and ranges (0–4000 m) were explored to determine the largest signal from the hydrate layer. Fig. 3 is a shaded anomaly plot of frequency versus range for the radial (Fig. 3 left-hand side) and azimuthal (Fig. 3 right-hand side) modes. The contours are electric field amplitude and the colour scale is the ratio of the hydrate layer response to a background sediment of 1 Ωm.
The largest normalized radial response from the hydrate model occurs at ranges and frequencies to the right of the thick white line at 10−15 VA−1 m−2, which is the instrument system noise floor, and are not measurable. The azimuthal mode has a significant normalized response at high frequencies that is below the noise floor of our instruments and so we only consider the radial mode. Despite this, a large hydrate signal is detectable at high frequencies (>10 Hz) and short ranges (<2000 m), above the noise floor. However, the electric fields attenuate very quickly at these high frequencies, which means that navigation of the transmitter becomes very important and the range window of detection is narrow. A 5 Hz square wave was chosen based on a compromise between the larger hydrate signal at higher frequencies and the stricter requirements of transmitter navigation at higher frequencies.
Other 1-D models of a thin hydrate layer and a thick hydrate/free gas layer were also examined to show that we can distinguish the existence and thickness of a hydrate/free gas layer (Weitemeyer et al. 2006a). However, Weinberger & Brown (2006) have shown that hydrate at Hydrate Ridge is distributed along fractures and permeable sand horizons, and is not distributed in a simple hydrate layer as depicted in this model study. Nonetheless, these 1-D model studies are useful for determining the CSEM frequencies, ranges and geometries we needed to detect hydrate.
The Hydrate Ridge Experiment conducted in August 2004 was an opportunistic use of shiptime aboard the R.V. New Horizon, during a transit from San Diego, California, to Newport, Oregon, for another project. Although we had funding for only 3 days on station for the survey, a marine EM experiment to image shallow gas hydrates was successfully conducted. The experimental layout, shown in Fig. 1(B), consisted of a single east–west line of 25 ocean bottom electromagnetic (OBEM) receivers spaced about 600 m apart and a CSEM tow line (CSEM tow 1) along this line to generate radial data. [An azimuthal tow was made 2 km parallel to the line of receivers, but navigation uncertainties resulted in these data being of little use (Weitemeyer 2008)]. The receivers and CSEM tow 1 are co-located with the ODP Leg 204 drill sites 1244, 1245, 1246 and 1252 (Shipboard Scientific Party 2003b) and seismic line 230 from the high resolution precruise 3-D seismic site survey (Tréhu & Bangs 2001), to provide ground truth for the EM results. In addition to the CSEM data, about half the receivers had magnetometer sensors to record Earth′s natural time varying MT signal during the experiment, which allowed us to image the accretionary complex sediments using this method.
The CSEM transmitter used for this study is a horizontal EM source similar to that described in Constable & Cox (1996), nominally capable of 200 A transmission—a moderate current output considering that SIO now has a transmitter capable of 500 A, and industry commonly transmits 1000 A when collecting CSEM data. Nevertheless, this is a sufficient current for imaging hydrate in the top 100s of metres, compared to the deeper targets (1000s of metres) for which CSEM technology has been commonly used. In fact a higher output will increase the saturation of the OBEM receivers at close ranges, where the electric fields are most sensitive to the hydrate. During CSEM tow 1 the transmitter, called Scripps Undersea Electromagnetic Source Instrument (SUESI) was flown approximately 100 m above the seafloor and transmitted a 100 A, 5 Hz square wave. The electric dipole moment of the transmitter is dependent on the separation of the two copper pipe electrodes (90 m) and the output current (100 A). The Fourier series representation of a square wave gives a coefficient of 4/π that multiplies the peak current at the fundamental frequency, and so the dipole moment for CSEM tow 1 was 11.5 kAm—incorrectly reported as 1.15 kAm in Weitemeyer et al. (2006b).
The transmitter was deep-towed at an average speed of 1.5 knots (46 metres per minute) at the end of an armored coaxial 17 mm (0.680 inch) cable that is used both to power the transmitter and for telemetry between the transmitter and the shipboard control console. A stand-alone AC power source takes shipboard 60 Hz three phase power and transforms it to 2000 V, 400 Hz single phase power. The 400 Hz frequency is generated under control from a GPS time base (Zyfer 565-210) to ensure that transmitter phase does not drift, a common problem in earlier academic sources. Power is transmitted down the tow cable with bi-directional FSK telemetry overlain at 70 kHz and 118 kHz. At the transmitter the 2000 V power is transformed down to about 100 V, and internal control circuitry, a set of rectifiers and bipolar transistors are used to generate a square wave with a lower frequency envelope. A display system of the environmental parameters of the transmitter (internal temperatures, current and navigational sensors) is monitored on board the ship.
SUESI is equipped with a Paroscientific Inc. depth sensor, a Valeport Limited CTDV (conductivity, temperature, depth and sound velocity) metre, a Kongsberg Simrad 1007 altimeter, and an acoustic transponder. These sensors, along with acoustic ranges from the ship to the transmitter, aid in locating the position of the transmitter as a function of time. The Valeport also measures an accurate conductivity, sound velocity and temperature profile of the seawater as the transmitter is lowered to its tow depth. The sound velocity profile is used for acoustic navigation of the receivers, and the seawater conductivity profile is used when inverting the CSEM data. The altimeter readout is continually monitored and digitally recorded so that the depth of the transmitter can be adjusted in real time to ensure a tow height of 100 m above the seafloor. These altitude measurements are used later in the modelling of the data.
The seafloor receiver is similar to the one described by Constable et al. (1998). Two types of receiver configurations were used: a vertical electric (VE) receiver (odd numbered sites) and a MT receiver (even numbered sites). All receivers consisted of two horizontal perpendicular sets of Ag–AgCl (silver–silver chloride) electrodes on a 10 m dipole. The VE receiver had additional Ag–AgCl electrodes along a vertical 1.5 m dipole. The MT instrument instead had the addition of two horizontal and orthogonal induction coil magnetometers, allowing the collection of MT data at every other site. Although the coils recorded to >10 Hz, the ocean acts like a low pass filter on the natural MT signal, only allowing frequencies below about 0.1 Hz to be detected. We did not combine the VE and MT sensors on one instrument because at the time we were concerned that the VE sensor would cause motion of the instrument and corrupt the MT data. The receivers record both the natural time varying MT field of the earth and the manmade source from the transmitter, as well as any ocean-generated environmental noise.
MT data and regional structure
MT data can provide important constraints on fluid distribution and mantle temperature associated with subduction zones when collected in long profiles (hundreds of kilometres) on both the seaward and landward side of the subduction zone (e.g. Evans et al. 2002). The MT profile at Hydrate Ridge is relatively short (20 km), and while sensitive to deep mantle structure does not resolve lateral structure beyond the limits of the profile in any detail. However, we are able to gain some insight into the shallow electrical conductivity structure of the accretionary complex sediments and obtain some evidence for the electrical conductivity structure of the mantle and crust from these data. The MT data were processed in two groups (s02–s12 and s16–s24) using the multi-station MT processing code of Egbert (1997); of the 12 MT receivers 10 had good MT responses (a median standard error of 0.2 Ωm in apparent resistivity and 1.8° in phase) at 10–1000 s period, using only 46 hrs of continuous data.
Impedance (Z) polar diagrams and Swift skews for periods between 10 and 1000 s are shown (Fig. 4) to provide an indication of the dimensionality of the data. Circular Zxy and small Zxx indicate 1-D structures, while the more peanut shaped Zxy indicates 2-D structures. When the Swift skew is above 0.2, 3-D effects need to be considered [see Simpson & Bahr (2005)]. These data are generally consistent with 2-D interpretation.
These MT data were inverted using a 2-D MT OCCAM inversion program of deGroot Hedlin & Constable (1990). A starting model of 10 Ωm was used and in 11 iterations a RMS misfit of 1.3 was achieved using a 10 per cent error floor for apparent resistivity and phase data. Fig. 5 shows the MT responses and the TE/TM mode fits of a smooth OCCAM inversion for all sites.
Fig. 6 shows the MT derived resistivities from a combined TE and TM mode inversion of the MT data with an overlain seismic interpretation by Gerdom et al. (2000). The agreement of the 15 km MT line with the seismic interpretation is surprising for such a short profile. The sediment thickness agrees, and the increase in resistivity of the lower crust/upper mantle is seen. This image suggests that with longer period data and receivers placed from the toe of the accretionary complex to the coast it will be possible to image the subduction zone well. For the purpose of the current experiment, the MT data provides the regional background resistivity structure while the CSEM data provides a detailed image of the shallow sediments in which gas hydrate is found. However, a second inversion using only the the top decade of periods (10–100 s) and a finer mesh was run with the aim at getting a better image of the shallow structures in the MT data, less effected by the coastlines and bathymetry. This result is shown in Fig. 9 and achieved a RMS of 1.12 from a starting model of 10 Ωm.
CSEM data processing
The CSEM data (d) were processed by sectioning the raw time-series (t) recorded by the receiver into 120 s stack frames and fitting sinusoids of an angular frequency ωi= 2πfi, where i is the index of the frequency (f):
The average tow speed of 1.5 knots means that the 2 min stack frames are averaging over a range of about 92 m between data points; this is approximately a movement of one dipole length (90 m). Only the 5 Hz and 15 Hz data have been analysed because the higher frequencies had a poor signal-to-noise ratio due to a combination of instrument performance, the 1/f falloff, and rapid attenuation with range (Weitemeyer 2008). Our discussions in this paper will focus only on the fundamental frequency of 5 Hz.
The data loggers store several hours of data in RAM and then write these data to disk. The write operation creates electronic noise lasting about a minute, but the processing routine ignores the coefficients at these times. Digital counts are converted into volts, using the least count (count V−1) of the analogue to digital converter (1165084.3 counts V−1 for MkII and 2516582.2 counts V−1 for Mk III instruments), and then the data are normalized by the receiver′s antenna length, the amplifier gain (1000000), and the amplifier transfer function at 5 Hz, to get electric field in volts per metre. The data are finally normalized by the transmitter dipole moment.
Examples of calibrated amplitude and phase data versus distance along transmitter UTM easting are shown for two sites in Fig. 7, for a VE and MT instrument. The closest approach of the transmitter to the receiver occurs at just about 330.1 km easting for site 9 and 330.8 km easting at site 10. Unfortunately, the electric field data are saturated at about 10−10 V A−1 m−2 at source–receiver ranges of 750 m and less because the electric field amplifier gains were set to 1000000, a hold-over from an earlier generation of 16 bit instruments. Subsequent work shows that the 24 bit instrument gain can be lowered without compromising long range data. The noise floor for the horizontal electric fields is about 10−15 V A−1 m−2 and for the VE field 10−14 V A−1 m−2, which is proportional to dipole length. The magnetic field sensors rarely saturate, as the gains were set lower than for the electric field amplifiers, and have a noise floor of about 10−17 T A−1 m−1. The noise floor for any given receiver is most evident in the scatter in the phase data. The transmitter navigational parameters of antenna position (x, y, z), antenna azimuth, antenna altitude and antenna dip as a function of time are merged with the receiver amplitude and phase data versus transmission time (see -2">Section 7.1.2). The 5 Hz data reached a noise floor of 10−15 V A−1 m−2 at about a 2.5 km range.
Initial interpretation using pseudosection analysis
The quality of marine CSEM data is dependent on accurate navigational information for the transmitter and the receiver positions and orientations. Because the relatively high CSEM frequencies used here attenuate rapidly with source–receiver separation, the navigational data for the transmitter and receivers had to be more accurate than for earlier academic crustal-scale CSEM experiments (e.g. Cox et al. 1986). Accurate navigational data were meant to be collected using a rented, commercial short baseline (SBL) acoustic navigation system, but unfortunately this system failed. Instead, long baseline (LBL) acoustic navigation data were collected using a backup system, recording ranges between ship, receivers and the transmitter.
Preliminary receiver positions were estimated using a Marquardt inversion (a non-linear least squares method) of LBL acoustic traveltimes collected during the CSEM transmission and during instrument recoveries, using ray-tracing and a variable seawater sound speed versus depth for the forward calculations. Due to the short time available for this experiment, there was no time for dedicated acoustic receiver navigation, and we estimate that receiver positions obtained this way are accurate to about 50 m compared to 3–5 m for fully surveyed instruments using LBL.
Recording compasses were mounted inside the logger pressure cases to measure instrument orientations, but these proved to be subject to field distortion from nearby magnetometers and batteries. This issue has since been resolved by making the compass external to the instrument, but in this case receiver orientations were estimated by using the geometry of polarization ellipses of the electric or magnetic fields (Behrens 2005). Forward modelling of the x and y components of electric and magnetic fields was used to resolve 180° ambiguities [see chapter 4 in Weitemeyer (2008) for more details].
The LBL acoustic navigation, direct acoustic pings, relays from receivers and layback of the transmitter from the ship based on transmitter depth and direct ranging to the transmitter provided limited accuracy for the cross-tow set and no indication of transmitter dip because we had not yet included a recording depth gauge on the antenna. However, the cross-tow position and dip of the transmitter can be estimated by modelling the electromagnetic (EM) radiation pattern of the horizontal electric dipole source. Close range electric and magnetic data collected by seafloor receivers (in this case <1.5 km in source–receiver offset) are less sensitive to seafloor resistivity than long range (>1.5 km) data, and can be used to refine the geometry of the transmitter and receivers. A Marquardt inversion was developed to solve for navigational parameters, including transmitter position, rotation, dip and receiver positions. The inversion program uses a 1-D dipole forward modelling code, Dipole1D (Key 2009), and requires an initial model of half-space seafloor resistivity and the geometry of the transmitter and receivers. The program updates the model parameters until convergence is reached between the synthetic EM responses and the observed EM data. We call this technique ‘total field navigation’ and applied it to the Hydrate Ridge data.
Because multiple receivers observe a single transmitter the inversion is well constrained. The model requires at least four free parameters to achieve good fits to the data: transmitter rotation, dip and x and y positions. In our case a model which found the transmitter x, y, dip, rotation and 22 out of 25 receiver (x, y) positions, gave an RMS misfit of about 2.79, with a 10 per cent error associated with each data point. We deem this a reasonable misfit and the model of the transmitter is smooth, with reasonable values for the rotation and dip of the transmitter antenna. While this method is useful to improve navigation in CSEM surveys it also demonstrates the sensitivity of each of the EM components to the geometry of the transmitter, especially at ranges shorter than 750 m.
Navigational errors in transmitter and receiver positions obtained this way are less than 10 m. The transmitter dip and azimuth are known to within 0.6° and 2°, respectively, which is fairly good for this generation of CSEM data (ca. 2004). Table 1 provides some insight into how these navigation errors influence interpretation, by giving the percentage change in electric field amplitude when transmitter geometry is perturbed by these errors. For comparison we include the effects of 1 per cent, 5 per cent and 10 per cent changes in resistivity of a 1 Ωm half-space. Resistivity estimates derived from the transmitter geometry will thus have an error equivalent to between a 1 per cent and 5 per cent error at all ranges considered. Error in the radial position has the largest effect, and the dip the second largest effect. The effect of our transmitter geometry error is thus equivalent to between 1 and 5 per cent error on estimated apparent resistivities. We assigned a 5 per cent error floor to the data amplitudes for the 2-D CSEM inversion shown below, which is large enough to capture the navigational errors.
To obtain an image of the subsurface structure and heterogeneity across the CSEM profile, a pseudosection technique (a method used extensively in land DC resistivity and IP surveys) was used. We converted the major axis of the polarization ellipse Pmax electric field amplitudes into apparent resistivities using models of half-space resistivities for each receiver–transmitter geometry. Models were computed using the 1-D layered code of Key (2009) given the water depth at each receiver, transmitter height, the UTM easting and northing positions of the transmitter and receiver, and the transmitter azimuth and dip. The pseudosection provides a way to look at all of the CSEM data collected at every site in one single image. The midpoint between the source and receiver is plotted at a depth given by a 45° projection from the source and receiver. Since EM induction is not a purely geometric phenomenon, this image is not a depth section.
Shaded apparent resistivity pseudosections are shown in Fig. 8. The values of apparent resistivity versus range varied between 0.3 and 3 Ωm. Reciprocity between the transmitter and receivers creates a two-fold redundancy in the data, with separate pseudosections from east-side and west-side transmissions (or, in-tow and out-tow in recent industry terminology).
The pseudosection projection technique causes the west pseudosection to have a striping pattern to the west and vice versa for the east pseudosection. This pattern is most obvious for the shallow conductor under site 6. However, the east-side and west-side pseudosections (top and middle panels) are sufficiently similar to justify taking an average of the two (bottom panel). The combination of the east and west images reduces the striping pattern except under site 6, where a classic pseudosection ‘pant-leg’ feature occurs where a surface conductor has been mapped into depth because of the data projection technique.
All three pseudosections display a more conductive basin under sites 18–24 that increases in resistivity with depth, likely a result of a decrease in porosity due to compaction. The pseudosections also show a resistive anomaly where an anticline is evident in the seismic data under sites 16 and 17. Sites 1 to 14 generally have a shallow surface conductor likely due to high porosity (Fig. 11), but below this are resistive features that may be associated with hydrate and free gas. Particularly, the ridge itself is resistive below sites 9–13, and there is a prominent resistor under sites 2–8, particularly evident in the west pseudosection.
A collaboration with EMI Schlumberger allowed access to a 2.5D finite difference CSEM inversion code called 2.5D DeepEM inversion (Gao et al. 2008) to interpret Hydrate Ridge CSEM survey. The major axis of the polarization ellipse (Pmax), and a relative coordinate system between transmitter and receiver was sufficient for pseudosection projections. However, the 2-D inversion code requires an absolute coordinate system (x, y, z) and so we use the radial electric field data for inversion. A full development for the selection of an appropriate finite difference grid to represent the bathymetry profile of the CSEM transect and the final inversion are presented in Weitemeyer et al. (2010). The inversion program was given 59 transmitter positions spaced about 240 m apart and 25 seafloor receivers spaced about 600 m apart. The observed in-line imaginary and real electric field data were assigned a noise floor of 5 per cent of the maximum datum amplitude. The 2.5D inversion achieved an RMS misfit of 4.73 in 22 iterations from a starting RMS of 12.03 (see Weitemeyer et al. 2010, for plots of data and fits).
The 2.5 DeepEM inversion model is shown in Fig. 9(A). The inversion provides a depth scale unattainable from the pseudosections, giving both lateral and vertical resistivity structure. The extent of the resistor to the east of site 4 was not obvious in the pseudosections because the conductive pant leg dominated the image. The pant leg has been collapsed to a surface conductor in the inversion, confirming that it was indeed an artefact of the pseudosection projection technique. The resistive feature associated with the seismic anticline under site 16 is still present, and the conductive basin to the east also remains. The inversion result includes a shallow conductive basin below sites 18–25, and a shallow resistor at about the depth of the BSR to the west below sites 1 to 7. Deeper in the inversion at about 1600 m there is evidence of folding in the accretionary prism. However, the MT data provide a much better image of the deep accretionary margin sediments.
Comparison of CSEM with MT results
Fig. 9(B) shows a close-up of the MT model with the colour scale saturated to match the CSEM model′s resistivity scale. The CSEM and MT inversions provide complimentary images, derived independently and with inherently different sensitivities. MT data are preferentially sensitive to conductive features and have sensitivity to depths of tens of kilometres, while the CSEM method is preferentially sensitive to resistive features shallower than a few kilometres. The resolution of the MT model is fundamentally lower than the CSEM model at the depths of overlap. However, MT data provide resolution at depths smaller than the shortest skin depth (here 2 km for 16 s data in 1 Ωm), if only because of near surface galvanic effects, and it appears there is an overlap of at least a few hundred metres between the MT and CSEM data. The estimated depth of investigation for a 5 Hz CSEM transmission, based on a 1-D Fréchet kernel analysis at a source–receiver range of 3 km, is about 1.1 km. The two methods also have different along strike sensitivities; CSEM data are not very sensitive to offline structures, whereas MT data are very sensitive to offline structure because of the ubiquity of the source field (Constable 2010). The two inversions are, however, consistent, despite different sensitivities of the two methods and the vastly different depths of penetration. The conductive basin from sites 16 to 24 is present in both data sets. The MT model characteristically highlights the conductors more than the resistors, but both models show the presence of folding associated with the accretionary complex, which was less apparent in the CSEM pseudosection projection technique. An anticline is obvious under site 16 in all images. The MT model includes a dipping conductor at around −2000 m that is less pronounced in the CSEM inversion, possibly because of a lack of depth penetration. Under site 20, the MT model places a conductive body where the CSEM places a resistor bounded by a conductor above and below at the same depth. The MT model has a shallow conductor to the west (sites 2–8), whereas the CSEM model has a shallow resistor around the BSR depth in this region. Both of these differences may reflect the limits in the resolution of the MT data and could be associated with shallow non-2-D structure. However, in both the CSEM and MT models we see a deep resistor likely associated with the compacted accretionary complex sediments. Both the MT model and CSEM model place a surface conductor just below site 6. It suggests that our use of a 2-D approximation for both the MT and CSEM is capturing the main features of the geologic and bathymetric structure. (Both the MT and the CSEM inversions have taken into account the bathymetry across the profile.) The conductance (conductivity of a layer multiplied by the vertical thickness) is equivalent for the surface layers and is comparable for both inversions; for example, below site 8, the CSEM conductance is 125 S for a 125 m thick conductor of 1 Sm−1 and the MT conductance is 125 S for a 62.5 m thick conductor of 2 S m−1. There is an increase in resistivity with depth as a result of compaction and lithologic changes. The MT data provides no detailed information about the shallow sediments where one finds gas hydrate, but rather provides information about the background resistivity structure of the accretionary complex sediments.
Comparison of CSEM with Logging While Drilling (LWD)
There are three ODP Leg 204 LWD boreholes along the CSEM profile with well logs that we can compare to the CSEM inversion. CSEM soundings provide a bulk resistivity measurement at a scale of tens of metres to a kilometre, compared with the detailed centimetre scale resolution of the well logs, and so we expect some differences between the two very different samplings of seafloor resistivity. For example, the logging while drilling (LWD) deep focussed resistivity measurement has penetration depths of 13 cm and vertical resolution of 5–8 cm (Shipboard Scientific Party 2003b), and the resistivity at the bit (RAB) has a lateral depth of investigation of 30 cm (Shipboard Scientific Party 2003b). The LWD resistivity data are used for comparison with the CSEM inversion resistivities rather than wireline logging because LWD takes place during excavation of the hole (or shortly after), providing a measurement of resistivity before conductive drilling fluids invade the formation (Schlumberger Oilfield Glossary, ) or hydrates destabilize. The LWD deep resistivity log also samples the greatest volume. The CSEM resistivity values in the inverted model closest to each of the ODP Leg 204 well logs (1245, 1244, 1246) are used to make a comparison between the logging data, shown in Fig. 10. The RAB image maps the electrical resistivity around the borehole wall and is shown here as unwrapped borehole images with a colour scale provided by ODP Leg 204 Scientific Party.
The CSEM inversion provides an average resistivity value compared to the many small changes in resistivity observed in LWD, and will never be able to provide the centimetre detail offered by LWD, yet there is broad general agreement between the CSEM and the LWD resistivities at sites 1246 and 1244. On the other hand site 1245 has a very obvious difference between the LWD and CSEM at the central portion of the log. The upper 100 m and last 75 m in depth appear to agree well with the logged resistivities, but in the region between 100 m and about 300 m depth there are differences of nearly a factor of two. The CSEM inversions give a relatively large resistivity value, about 2.75 Ωm at the location of seismic horizon A, whereas the wells show an average resistivity of about 1.5Ωm and the RAB shows a high electrical resistivity in this region. The differences between the CSEM and LWD could be because CSEM is sensitive to the cummulative effect of all the resistive layers observed in the RAB, and is thus seeing the bulk effect of all these thin resistors. The CSEM inversion will also smooth the effect of any thin layers. Considering that horizon A is known to be a fluid conduit carrying quantities of free gas to the summit (Tréhu et al. 2004a), it is also possible that it is being seen as a resistor by the CSEM data but not by the well logs. Another consideration is electrical anisotropy from the inherent heterogeneity of gas hydrate distributions. CSEM methods measure mostly vertical resistivity but the well logs are mostly sensitive to horizontal conductivity. Seismic anisotropy has been attributed to hydrate veins and distribution and fabric of gas hydrates, with a vertical velocity higher than the horizontal velocity due to preferred alignment of hydrate veins (Kumar et al. 2006). In addition, Weinberger & Brown (2006) examined the electrical resistivity well logs from Hydrate Ridge and found that the hydrate distribution at ODP sites 1244–1246 is within well ordered hydrofractures dipping 20–70° and strike parallel to the trend of the ridge, and within permeable sand horizons that intersect the fractures. Since CSEM data are most sensitive to vertical resistivity in the overflight direction (Newman et al. 2010), and there is vertical resistivity variation due to the presence of horizon A, a type of macroanisotropy exists such that the bedding or layering is thinner than the resolution of the CSEM measurement. Therefore CSEM may ‘see’ a homogeneous but macroscopically anisotropic formation, whereas well logs do not as they sample a smaller sediment volume. A similar effect was seen in VSP data collected during ODP Leg 204, where a velocity–thickness product satisfied the traveltimes equally well at horizon A: either a 15 m thick region with velocities of 1100 ms−1 or a 4 m thick region with velocities of 595 ms−1 surrounded by 1600 ms−1 sediments (Tréhu et al. 2006b). From the direct sampling in well logs we know horizon A to be about 4 m thick, and so the CSEM inversion may also be trading off the resistivity thickness product in a similar way to the VSP data.
Site 1252 has no LWD measurements because it was an add-on at the end of the leg and all LWD were collected at the beginning of the leg. There are wireline logging induction measurements, but any hydrate that was present would have largely been disturbed or dissociated. The well log resistivity measurements are in general agreement with the CSEM inversion results but CSEM results are slightly more resistive throughout the sequence, consistent with a sensitivity to in situ hydrate.
Comparison of CSEM with 3-D Seismic Reflection Line 230
The colocation of the CSEM tow 1 with seismic reflection line 230 (Tréhu & Bangs 2001) allows us to make a comparison between the CSEM inversion and the seismic reflection, seismic stratigraphy and tomographic velocity data (Figs 11A, B and 12). Fig. 11(A) shows an overlay of the CSEM inversion with seismic line 230 and Fig. 11(B) the overlay of seismic sequence stratigraphy from Chevallier et al. (2006). The seismic horizons and geologic features discussed below are labelled according to Fig. 11(B). Seismic sequence stratigraphy of Chevallier et al. (2006) identifies geologic structures: Fold F, faults F1 and F2, anticline B and a dome; and two main types of sedimentary layers: (1) those that originated from the accretion of abyssal plan sediments of the Astoria Fan form the core of the accretionary complex; and (2) the subsequent deposition of the overlying slope basin sediments formed from the formation and evolution of the accretionary wedge fold-thrust belt system. These two main sedimentary packages are distinguished by the gradual change in resistivity with depth across the profile from the more resistive sediments at depth to more conductive shallow sediments. The top of the accretionary complex is marked with a thick grey line that consists of a major angular unconformity K and the fault F2 (Chevallier et al. 2006). The resistivity increase in the accretionary complex sediments is due to lithification and decrease in sediment porosity of these deep sea fan sediments. Anticline B is cored by a resistive anomaly from the accretionary complex sediments, and may also contain free gas. Fold F is cored by a splay fault F1 and is composed of younger deep sea fan sediments with a nannofossil-rich clay stone interspersed with turbidites and wood fragments that were rapidly deposited (references within Chevallier et al. 2006) and may be more conductive due to the presence of clays and higher porosity due to its rapid deposition. In sediments above the accretionary complex a resistor to the west of the profile is resolved by the inversion at about the depth of the seismic BSR, which typically marks the phase change from solid hydrate above and free gas below. In this region of Hydrate Ridge, the methane in the hydrates is biogenic and a concentration of gas at the BSR is created from this microbial methane production (Shipboard Scientific Party 2003a). The chaotic seismic region between sites 2 and 4 (Fig. 11A) was interpreted as having high free gas or gas hydrate saturations in an inversion by Zhang & McMechan (2003). The resistor in the same region is above the BSR and thus probably associated with high hydrate concentrations. There is also evidence of a resistive region (below sites 4–8) likely associated with free gas below the BSR. The shallow resistor between sites 6 and 7 may correspond to seismic horizons Y, a regional geologic unconformity (Chevallier et al. 2006). Seismic horizon A is a gas-charged fluid conduit taking methane gas to the southern summit (out of the page) (Tréhu et al. 2004a), which also shows up as a resistor in the CSEM inversion. Seismic horizons B and B' are largely faulted volcanic ash-lined conduits carrying free gas into the gas hydrate stability zone, which then freezes into hydrate (Tréhu et al. 2004b). These have a subtle resistivity contrast below ODP Site 1244 and not much of one below ODP Site 1246, suggesting that the bulk volume of hydrate decreases upwards and is too small for the CSEM to detect. A conductive region exists within the hydrate stability zone at the summit of this profile, suggesting lower hydrate concentrations and/or the presence of brines. A more conductive anomaly occurs at the surface below site 6 which may be associated with a recent formation of hydrate which causes the expulsion of salt. An EM survey over hydrates in the Gulf of Mexico also found conductive regions attributed to high salinity and temperatures, due to the rapid formation of hydrate (Ellis et al. 2008). The shallow sedimentary units have various resistivities which result from hydrate occurrence and changes in lithology that consist of turbidites and silty clays as the depocenters are filled. Two debris flows DBF1 and DBF2 occur on either side of anticline B and are characteristically conductive, probably due to the unconsolidated nature of sediments within debris flows.
Comparison of CSEM with tomographic seismic velocity inversion
Fig. 12 shows a vertical slice extracted from a 3-D velocity model derived from first arrivals recorded at 20 ocean bottom seismometers deployed during a seismic experiment to collect 3-D reflectivity data in preparation for ODP Leg 204 (Tréhu & Bangs 2001; Arsenault et al. 2001). Velocities within a 20×25×5 km volume were determined using the FAST algorithm (Zelt & Barton 1998). Inversion grid spacing was 0.5 km horizontally and 0.2 km vertically. After eight iterations, the RMS misfit of 14000 first arrival traveltime picks decreased from 203 ms to 24 ms. One of the most striking features of the model is a region of very low velocities northwest of the summit of south Hydrate Ridge. These very low velocities are required to fit the data and are interpreted to indicate the presence of free gas beneath the BSR in this region, which includes the gas-charged Horizon A, which was imaged by the 3-D seismic reflection data and sampled at three drill sites (1245, 1247 and 1250) during ODP Leg 204 (Tréhu et al. 2006a). The low velocity region coincides with the resistive region below sites 3–7, but because P-wave velocity is smoothed in the tomographic model, which was constrained by a minimum bound of 1.5 kms−1 on the velocity beneath the seafloor, the boundaries of the low velocity and high resistivity regions do not match in detail. Unconstrained inversions resulted in a narrow zone of lower velocity. However, folds in the underlying accretionary complex, as indicated by the 2.1 kms−1 boundary and the 2 Ωm boundary in Fig. 12, match surprisingly well between the velocity and resistivity models.
We follow the methodology presented in the Shipboard Scientific Party (2003b) and Collett & Ladd (2000) to estimate hydrate concentrations from both the deep resistivity well log data and CSEM resistivity-depth profiles at three Sites: 1244, 1245 and 1246. We begin by defining Archie′s LawArchie 1942).
Although Archie′s Law was developed to estimate water saturations in gas–oil–water–matrix systems (Archie 1942), it has been used to obtain quantitative hydrate concentrations assuming that hydrate fills the remaining pore space, as demonstrated in Collett & Ladd (2000), ODP Leg 204 Initial Reports (Shipboard Scientific Party 2003b), and (Collett 1998). A hydrate saturation, Sh, is calculated using
Archie′s Law requires knowledge of the saturation exponent, n, which is dependent on pore shape, connectivity, constrictivity of the pore network and distribution of the conducting phase (Spangenberg 2001). We used the coefficients reported in Tréhu et al. (2003) (a= 1, m= 2.8, n= 1.9386) to compute Sw. We also require ϕ, Rw and Rt-depth profiles to compute a hydrate saturation. The ϕ-depth profile is computed from the bulk density, ρb, measured by the azimuthal density–neutron tool, the grain density, ρm, measured from the cores and a fixed density for water (ρw= 1.05 gcm−3), as was done for the Shipboard Scientific Party (2003b)
R w is computed by converting interstitial pore water salinity at the top of the core into a resistivity () as a function of depth at a surface temperature, T1, using the algorithm of Fofonoff (1985; Fofonoff & Millard 1983) and Arps Law is then used to extrapolate the surface resistivity to depth using the temperature gradient (Arps 1953; Collett & Ladd 2000)
Using eq. (6) we get the water saturation and then eq. (7) gives the hydrate saturations. Fig. 13 shows the hydrate saturations for the CSEM and deep resistivity LWD data—note that the hydrate saturations are only valid within the gas hydrate stability zone, above the BSR. Below this region one could consider it to be a free gas saturation. Our Sh saturations differ somewhat from those reported in ODP Leg 204, because in the ODP Leg 204 Initial Reports computations were not made when the differential caliper log was of poor quality (>1 inch), whereas we placed no constraint on the log quality. The fit for the Site 1244 is excellent, including the fit for the accretionary complex rocks, which are quite different from the slope basin fill. The accretionary complex rocks are less dense, less resistive and have caliper logs of >1 inch (Tréhu et al. 2003). The fit at site 1245 is not so good below the BSR, but this may be the result of free gas beneath the BSR and within horizon A, as seismic velocities are also quite low.
The resistive region below sites 2 to 4 is thought to contain hydrate above the BSR and free gas below. There is no ground truth from well logs to prove this, and the seismic inversion indicates a LVZ suggestive of free gas. Despite this we venture to speculate on the hydrate saturation for this resistive zone (Rt= 2.5–4 Ωm). We again use Archie′s Law with ϕ= 55 per cent, Rw= 0.25 Ωm, a= 1, m= 2.8, n= 1.9386, based on the porosity and connate water resistivities found in site 1245 at about this depth (≈130 mbsf). From this we calculated that about 49 per cent of the pore space is filled with hydrate, which is about 27 per cent of the bulk volume. Archie′s Law is valid if hydrate is disseminated through the pore space, as hydrate will act to reduce the porosity (Hyndman et al. 1999). However, gas hydrate occurrence will not necessarily be modelled correctly by simple mixing rules (Lee & Collett 2001), and depending on the geometric distribution of hydrate, Archie′s Law may not be a representative model, especially if hydrate is found in veins and fractures, as was seen by Weinberger & Brown (2006) and Cook et al. (2008). The extremal bounds for effective conductivity σ are the Hashin Shtrikman bounds (HS-bounds) (Schmeling 1986; Hashin & Shtrikman 1963)Hashin & Shtrikman 1963). In terms of hydrate, the HS lower bound may represent a low concentration of granular disseminated hydrate distributed in isolated spheres within the conductive sediment. In clay-rich sediments hydrate may occur in veins or fractures and be better represented by the HS upper bound—where resistive material occurs in sheets impeding current flow through the matrix of fluid. The volume fraction of fluid can be related to the porosity and water saturation via
Discussion and conclusions
The use of CSEM data to augment seismic data and resistivity well logs provides an additional attribute to seismic velocities for estimating the volume of hydrate and provides resistivity information beyond a single vertical profile at a single well location. Our study at Hydrate Ridge represented an opportunistic use of ship time with only 3 days of data collection on station. For these reasons it should be considered a pilot study aimed to demonstrate that CSEM is a valuable tool in discriminating hydrate and indicate where improvements to experimental methods should be made (e.g. navigation and data quality). The original treatment of the data as 1-D apparent resistivity pseudosections was an effective way to determine the nature of the electrical heterogeneity across the tow line, but 2-D inversion is the only way to resolve the depth extent, particularly the relationship between resistive regions and the seismic BSR.
Marine EM methods have produced an electrical resistivity image of the subsurface consistent with models for hydrate emplacement at Hydrate Ridge. Both the focused high flux regime, as represented by the resistive region at horizon A, and a distributed low flux regime, as represented by the resistor along the length of the BSR to the west, are present in the models. The agreement with an independently obtained seismic velocity model is excellent and highlights the effect of free gas to depress p-wave velocity and increase resistivity.
While the EM data collected at Hydrate Ridge were meant to target hydrate, the simultaneous collection of MT data allowed us to obtain an image of the deeper folding of the accretionary complex associated with the subduction zone. Although the footprint is too small to image the subducting plate, these data suggest the potential for a much larger scale survey to successfully image the subducting slab. Our study demonstrates how marine EM can be used to understand the relationship between hydrate, gas and fluid flow in the accretionary complex.
Funds for ship time were provided by ExxonMobil and GERD, Japan. We thank the Captain and Crew of the RV New Horizon and the members of the SIO marine EM Laboratory. David Alumbaugh and Guozhong Gao of EMI-Schumberger are thanked for providing access to the 2.5D Deep Pixel inversion code and help in inverting the data. Kerry Key is thanked for assistance with MT processing. K.W. was supported by the Department of Energy grant DE-NT0005668. We thank editor Martin Unsworth and reviewers Nigel Edwards and Marion Jegen for their suggestions, which greatly improved the manuscript.