Estimating the composition of gas hydrate using 3 D seismic data from Penghu Canyon , offshore Taiwan

Direct measurements of gas composition by drilling at a few hundred meters below seafloor can be costly, and a remote sensing method may be preferable. The hydrate occurrence is seismically shown by a bottom-simulating reflection (BSR) which is generally indicative of the base of the hydrate stability zone. With a good temperature profile from the seafloor to the depth of the BSR, a near-correct hydrate phase diagram can be calculated, which can be directly related to the hydrate composition. However, in the areas with high topographic anomalies of seafloor, the temperature profile is usually poorly defined, with scattered data. Here we used a remote method to reduce such scattering. We derived gas composition of hydrate in stability zone and reduced the scattering by considering depth-dependent geothermal conductivity and topographic corrections. Using 3D seismic data at the Penghu canyon, offshore SW Taiwan, we corrected for topographic focusing through 3D numerical thermal modeling. A temperature profile was fitted with a depth-dependent geothermal gradient, considering the increasing thermal conductivity with depth. Using a pore-water salinity of 2%, we constructed a gas hydrate phase model composed of 99% methane and 1% ethane to derive a temperature depth profile consistent with the seafloor temperature from in-situ measurements, and geochemical analyses of the pore fluids. The high methane content suggests predominantly biogenic source. The derived regional geothermal gradient is 40°C km-1. This method can be applied to other comparable marine environment to better constrain the composition of gas hydrate from BSR in a seismic data, in absence of direct sampling. Article history: Received 29 November 2016 Revised 5 September 2017 Accepted 11 September 2017


InTroDuCTIon
A bottom-simulating reflection (BSR), is commonly associated with the presence of gas hydrate under the seafloor.The BSR generally indicates the boundary between free gasbearing sediments below and hydrate-bearing sediments above (Shipley et al. 1979;Yamano et al. 1982;Singh et al. 1993).BSRs can be observed in seismic reflection data due to the acoustic impedance contrast between sediments bearing hydrate and free gas (Kvenvolden 1995).The hydro-static pressure at the BSR can be estimated from the water depth and the BSR sub-bottom depth.Using the gas hydrate phase boundary, the temperature at the BSR can then be estimated (Yamano et al. 1982).
For a known seabed temperature and a linear geotherm at the site of each BSR depth measurement, a regional geothermal gradient pattern can be derived more efficiently than that can be done by in-situ measurements.Using the seafloor temperature, depth, and temperature at the BSR, a geothermal gradient can be determined.The depth of the BSR can be interpreted from 2D (e.g., Yamano et al. 1982;Chi et al. 1998;Chi and Reed 2008), and occasionally 3D Terr. Atmos. Ocean. Sci., Vol. 29, No. 2, 105-115, April 2018 seismic reflection datasets, dependent on the availability of data (e.g., Martin et al. 2004;Hornbach et al. 2012;Zander et al. 2017).This interpreted BSR depth can then be used to derive the temperature at depth of the BSR using a given hydrate phase boundary (e.g., Yamano et al. 1982;Chi et al. 1998;Chi and Reed 2008).
To accurately estimate temperatures at the depth of the BSR, we need to have an accurate phase boundary, which is dependent on the composition of gas forming the hydrate, and salinity of the pore fluid.In most cases, gas hydrate is formed by pure methane but other gases may also be present (Kvenvolden 1995).The variation of gas composition perturbs the system by changing the phase boundary, and in some cases it even changes the structure of the hydrate (Lu et al. 2007).Another factor that will affect the phase boundary is the salinity of pore fluids at the depth of the BSR, for which the exact alterations are still under debate.Some studies have used the salinity of sea water in the pore fluid (e.g., Grevemeyer and Villinger 2001;Shyu et al. 2006), while others have used the salinity of fresh water (e.g., Hyndman et al. 1992), due to possible repeating formation and dissociation of hydrate at depth, i.e., the formation of the hydrate will extract fresh water into the lattices, and such water will be released when hydrate dissociates.Information on the salinity of the pore fluid at the BSR depth can be important for marine electric resistivity studies under the seabed.Capillary effects can also change the depth of the phase boundary (Henry et al. 1999;Daigle and Dugan 2011;Liu and Flemings 2011), although such effects are thought to be minimal.Minshull and Keddie (2010) proposed an innovative approach of calculating the geotherm at a gas hydrate BSR using high resolution 3D seismic data considering several possible hydrate compositions.Several subsequent studies have discussed similar approaches (e.g., Daigle and Dugan 2011;Hornbach et al. 2012;Serié et al. 2017).If the water depth and the BSR depth are known, pressure at the BSR depth can be estimated.In addition, for a known gas hydrate phase boundary, the temperature at the depth of the BSR can be calculated.The pressure at the phase boundary increases with increasing water depth.Hence, for different sub-bottom depths of the BSR, we can derive a temperature-depth profile, which is a profile showing the variation of temperature at different sub-bottom depths.An accurate phase boundary is required to correctly extrapolate seafloor temperature.When the temperature profile is linear or follows a known nonlinear (exponential or cubic) function then the regression of BSR temperature as a function of BSR sub-bottom depth can be extrapolated to derive seafloor temperature, which is the temperature at zero BSR sub-bottom depth.The correct phase boundary must give an accurate estimate of seafloor temperature based on the regression function and from in-situ seafloor measurements.A small variation in the chemical composition of the gas and pore fluid salinity will shift the phase boundary by several degrees at the seafloor (Kvenvolden 1995;Henry et al. 1999;Minshull and Keddie 2010).This relationship provides an opportunity to constrain the correct phase boundary.This method can be used when dense seismic data collected over a small region with a bathymetric relief are available.For a large areal extent of seismic data, the geothermal gradient might vary with water depth for other reasons, such as heat transfer due to the temperature field disturbance from displacement of faulted blocks (e.g., Ganguly et al. 2000).However, if the BSR is densely sampled at a small site with large bathymetric relief, the regional tectonic processes, which are assumed to have longer spatial wavelengths, should not vary significantly over the area.But such bathymetric relief can cause large lateral temperature changes in short wavelengths due to topographic focusing of heat flow from depths (e.g., Kinoshita et al. 2011), and might need correction as we wish to document in this work.Topographic focusing of heat flow occurs due to the fact that some heat from depth will transfer through the shortest path to the seafloor, causing higher heat flow in the submarine canyons; similarly defocusing effect will conversely lead to lower heat flow in ridges (e.g., Birch 1950;Turcotte and Schubert 1982;Shankar et al. 2010).In this study, the effect of a poorly fitting curve on the phase boundary may be separated from that of larger-scale regional thermal structures due to the small footprint of the dataset used in this study.Some difficulties may arise when applying this method.One issue is that widespread variation of different BSR temperatures at a given sub-bottom depth makes the regression trend of a temperature profile less well constrained.In this study, we demonstrate that the scattering in temperature can be reduced by applying a correction to account for the effect of topography on the temperature field.
In this study, we first mapped the BSR across a 3D seismic dataset.We then performed a topographic correction based on our 3D modeling method using the depths of the BSR from seismic data and measured seafloor temperature to derive a 1D temperature-depth profile in the area, and thus, the gas composition of hydrate at depths of the BSR.The temperature-depth profile is constrained by insitu seafloor temperature measurements.We then derived the best-fitting regional geothermal gradient for this area.We wish to document the improved results from the 3D bathymetric corrections and the non-linear temperature profiles used in this work.This method can be applied to comparable marine environment (e.g., other passive margins) to indirectly estimate gas hydrate composition from an observable BSR in a 3D seismic data set.

GEoloGICAl sETTInG
Taiwan is located along the boundary between the Eurasian plate and the Philippine Sea plate where oceanic lithosphere of the South China Sea is subducting eastward beneath the Philippine Sea plate (Bowin et al. 1978) at a rate of about 8 cm yr -1 (Yu et al. 1997), forming the Manila trench, the Taiwan accretionary prism, the Luzon Trough (forearc basin), and the Luzon Arc.The subduction changes into arc-continent collision where the thicker and more buoyant Chinese passive margin enters into the convergent zone.The Penghu canyon is located 80 km offshore SW Taiwan on the passive margin that is the foreland basin of Taiwan (Fig. 1).Along the passive margin, there are several similar canyons and other erosional features running towards the south.The BSR is wide-spread in both active and passive margins of this region (e.g., Chi et al. 1998;Chi and Reed 2008;Lin et al. 2009).

DATAsETs AnD METhoDs
A pseudo-3D seismic dataset (Fig. 1) was collected at the upper Penghu channel area, offshore SW Taiwan by National Taiwan University in 2013 on the R/V Ocean Researcher I (OR1).A 475 in 3 air-gun array was used as source.A single 1037 m long streamer was used, composed of 84 channels, each with a channel spacing of 12.5 m, and a shot interval of 25 m.In total, 78 in-line profiles were collected with a line spacing of 50 m.The area of this pseudo-3D cube is 17.5 km by 3.85 km.The seismic data were processed using ProMAX software.The seismic processing workflow included a band-pass filter, normal move out correction, 2D/3D geometry, 2D/3D stacking, 2D/3D migration, and auto gain control (Han et al. 2017).The BSR is widespread in this region (Lin et al. 2009(Lin et al. , 2014)), including the area of seismic data coverage (Fig. 1b) for this work.Figure 2 shows an example of an in-line profile and the location of the BSR.
We used the BSR sub-bottom depth and water depth to calculate pressure at depth of the BSR assuming a hydrostatic gradient (Yamano et al. 1982).The BSR sub-bottom depth was calculated from two-way travel time using a velocity function derived from velocity analysis of the seismic data: where t = BSR(t) -SF(t), BSR_mbsf is depth to the BSR (in mbsf), "BSR(t)" is the BSR depth (one-way travel time, s), and "SF(t)" is the seafloor depth (one-way travel time, s).Using pore pressure, temperature at depth of the BSR was calculated using a phase boundary for varying gas compositions and an assumed ratio of pore water salinity (Fig. 3), based on Tohidi et al. (1995).These phase boundaries are usually nonlinear and thus we do not show these thermodynamic equations.For each composition, we used the program Heriot-Watt Hydrate (HWHYD), developed at Heriot-Watt University (Tohidi et al. 1995;Webber et al. 2007), to calculate a range of pressures and temperatures at the phase boundary.We then selected only those phase boundaries that match the pressure at depth of the BSR.We  used a porewater salinity of 2% based on samples from the Shenhu Area, northern South China Sea (Wu et al. 2011), which are located in a similar tectonic setting to the west of our study area.
We corrected for the effect of topography on geothermal gradients in sediments inside the 3D seismic cube in the shallow crust (Minshull and Keddie 2010).The seafloor geothermal gradient can be lower under ridges but high under valleys due to the deflection of the heat flux from depths to the valleys away from the ridges.In other words, it is easier for heat from depth to reach to the submarine valleys due to the shorter distance, compared to the distance to the ridges.Such variation is not related to the chemical composition of the hydrate and thus must be corrected, even though such correction has not been done before for this type of studies.
Here we use the high-resolution bathymetry data with a grid size of 50 m (S.-K.Hsu, personal communication) acquired using a deep-towed multibeam echosounder system.We used a 3D finite element method (Braun 2003) to calculate the temperature field in the 3D seismic cube (Fig. 4) using the following boundary conditions: (1) a constant seafloor temperature of 4.13°C (Shyu et al. 2006), based on the in-situ thermal probe data and (2) four insulated side boundaries.A regional geothermal gradient of ~40°C km -1 was used (Shyu et al. 2006) after many iterations of trial and error testing of different geothermal gradients during the modeling, to fit the BSR depth in the 3D seismic cube under the continental shelf.This area was used as the bathymetry is relatively flat, and therefore topographic focusing should have minimal influence here.As we only want to derive the first-order effect, we assume a homogeneous crust and the topographic effect is in a steady state.A 3D heat transfer equation for an isotropic and homogeneous medium can be defined as: where T is temperature, t is time, and a is thermal diffusivity, which can be defined as k c p t , where t is density, c p is specific heat capacity, and k is thermal conductivity.Due to the steady-state assumption, time dependent variations are zero, therefore: This shows that for this particular type of modeling the temperature field is controlled by the boundary conditions, but not related to the thermal diffusivity, and thus thermal conductivity, used in our model.

rEsulTs
We calculated the topographic effect across the extent of the BSR within the 3D seismic cube.Firstly, we acquired temperature and geothermal gradient fields at the BSR locations (T BSR_num_topo and GG BSR_num_topo ) from the numerical thermal model.The temperature field is usually different from the temperature derived from depths of the BSR (T BSR_no_topo ), which was forward-modeled based on a 99% methane hydrate phase diagram using water depths, subbottom depths of the BSR, and a regional geothermal gradient of 40°C km -1 .We assumed that the discrepancy from the seismic data and from the 3D numerical modeling was primarily the result of the topographic effect.We calculated the topographic effect coefficient (TEC), which is the ratio The BSR sub-bottom depth is greater under the ridge, but shallower under the flanks, due to the deflection of heat source from depth away from the ridge to the flanks.Such topographic effects can be corrected using the information from Fig. 4a, allowing us to derive a consistent geothermal gradient, i.e., the same BSR sub-bottom depth for the whole BSR, as shown by the corrected line.between the regional geothermal gradient of 40°C km -1 and GG BSR_num_topo .In other words, TEC represents the predicted effect of topography on temperature field at the BSR depth based on our numerical modeling.Subsequently, we applied a correction to the observed data by multiplying the TEC with each particular BSR-derived geothermal gradient to correct for the topographic effects on the seismically-derived temperatures at the BSR depths.After 3D topographic correction, the scatter of BSR-derived temperatures at various sub-bottom depths is reduced, supporting the theory that this correction method can be used to better constrain the results.However, Hornbach et al. ( 2012) have proposed that these variations of BSR sub-bottom depths can also be caused by lateral changes of hydrate concentrations above the BSR.Here we assume there is negligible variation in hydrate concentration in this region because we do not observe clear velocity pull-up in this seismic dataset.
After topographic correction, the temperature profile was fitted using a least square regression to derive a geothermal gradient (Fig. 5).The results show that the geothermal gradient varies with depth.We then tried to derive extrapolated seafloor temperature data using several different methods.At first, we used a linear fit to the data where the BSR sub-bottom depth is between 100 to 170 mbsf.However, the extrapolation of resulting seafloor temperatures by linear fitting is much higher than the observed value (Fig. 5).Next, we use a cubic function rather than a linear function; however, the extrapolation result for this also gave higher seafloor temperatures than in-situ measurements (Fig. 5).In the third estimation, we considered decreasing porosity, and thus, increasing thermal conductivity with depth.This results in a higher geothermal gradient in shallower sediments (Fig. 5).Once we considered the topographic correction and depth-dependent thermal conductivity, we found the bestfitting trend was represented by a quadratic function.0.00010429 0.070726 4.1468 .
where T is temperature in °C and d is the sub-bottom depth in meters.The depth-dependent temperature profile trend in the best-fitting model is similar to that used by/described in Chi and Reed (2008).Chi and Reed (2008) derived the geothermal gradient variation with depth due to increasing thermal conductivity, this was obtained using depth-dependent porosity profiles from several nearby wells (Lin et al. 2003).This provided a more accurate seafloor temperature than the linear regression method as the estimated temperature agrees with the measured temperatures from multiple in-situ thermal probe measurements, indicating a temperature of 4.13°C (Shyu et al. 2006).
Next, we discuss how we tested the robustness of our Fig. 5. BSR sub-bottom depths vs. temperatures calculated at the depth of the BSR for a gas composition of 1.2% ethane and 98.8% methane, with a pore-water salinity of 2%.Various curves were used for fitting to derive our preferred model.Here we show the regression results for the linear, quadratic, and cubic fits: T = -0.029849d+ 7.8694, T = -0.00010429d 2 + 0.070726d + 4.1468; T = -0.00000057493d 3 -0.00022748d 2 + 0.010619d + 7.4959, respectively.The root-mean square deviations are 0.81, 0.85, and 0.85°C, respectively.The linear fit has smaller deviation but the extrapolated seafloor temperature is too high.The best fit is generally consistent with the one based on Chi and Reed (2008)'s data, which was derived from a regional dataset of over 1000 BSR locations, taking porosity decrease with depth into consideration.See the text for details.
best model.We tested a range of different compositions with the resulting temperature profiles (Fig. 6).These calculations consistently showed results that indicate the presence of a high concentration of methane with up to 1% ethane.The presence of hydrogen sulphide increased the calculated seafloor temperature significantly while an increase in salinity decreased the calculated seafloor temperature (Fig. 6a).
The effect of changes in salinity and gas composition had a non-uniform effect on the geothermal gradient (Fig. 6b), for example, a general decreasing trend in the geothermal gradient with an increase in percentage of other gases.Carbon dioxide has a relatively small effect on the inferred geother-mal gradient; therefore, we have discounted its influence on the seafloor temperature in this case.Unfortunately, this test also shows that we cannot exclude the possibility of having a higher percentage of ethane if the pore fluid has a high salinity.
We conclude that the mean geothermal gradient for the BSRs at a depth of 100 -170 mbsf in the study area is 40.1°Ckm -1 (Fig. 7).Previous studies have reported gradients around 35 -55°C km -1 (Shyu et al. 2006;Chi and Reed 2008) for locations near that of the present study.
Using a molecular composition of pure methane results in a significantly lower seafloor temperature than the extrapolated value.Based on our model we therefore infer that at this site, the composition of the gas hydrate is 99% methane and 1% ethane (Figs. 3 and 5).
There are ongoing debates as to whether the gas composition near this site is biogenic or thermogenic.A systematic study of pore fluids in the sediments recovered from many gravity cores found that the gas is composed of pure methane (Chuang et al. 2006).However, Wu et al. (2011) found that the methane content in volume percentage is 96.10 -99.91% with small amount of ethane and propane in this region.Zhu et al. (2009) and Chen et al. (2013) documented biogenic methane with some possible percentage of thermogenic methane at shallow depths for a site to the west of our study area.All those results were from core samples from the seafloor, down to 20 m below the seafloor (mbsf).These results didn't have any date from higher depth as typical piston and gravity coring systems cannot sample such depths.Our study can supplement these results as we use the seismic data to derive the gas composition at higher depth of about 100 -300 mbsf.Our results of 99% methane extends these finding of very high methane concentration from upto 20 mbsf (e.g., Zhu et al. 2009;Chen et al. 2013) to about 100 -300 mbsf.

DIsCussIon
In this study, we use a model based on an adaptation of the method described by Minshull and Keddie (2010) to constrain estimates of gas hydrate composition.As a result, we infer that at this site, the composition of the gas hydrate is 99% methane and 1% ethane.This is consistent with previous work offshore Taiwan by Chuang et al. (2006) in a largescale geochemical study, at Formosa Ridge on the same passive margin to the west (Feng and Chen 2015), and with a slightly higher methane concentration from the study of core samples from Shenhu Area, in the northern South China Sea (Wu et al. 2011), located in similar tectonic setting to the west of our study area.If the pore fluid salinity differs from 2% then the composition of gas hydrate may vary.However, 2% salinity was used based on core samples from the Shenhu Area, northern South China Sea (Wu et al. 2011), which are located in a similar tectonic setting to the west of our study area and the results of our geochemical modelling are consistent with previous studies (e.g., Feng and Chen 2015).
Unlike previous examples, in this study we compensated for the effect of 3D bathymetry on the temperature at the BSR based on a static 3D finite element model.Through this correction, we reduced the scattering of the data in the temperature profile.We also investigated the theory that if the variation in conductivity with depth is known, more accurate estimates can be obtained using curve fitting methods based on the conductivity profile.We found such a depthdependent temperature profile is needed to better fit the in-situ seafloor measurements in this setting.We derived a geothermal gradient of 40°C km -1 , which is consistent with the regional geothermal gradient of 40°C km -1 used in the 3D finite element thermal modelling in this study, and is around the average value for the BSR-derived geothermal gradients in this region.These results suggest that topographic affect may be significant in this region in terms of seafloor geothermal gradient.In addition, because thermal conductivity is treated as depth-dependent, the derived geothermal gradient provides an average value that can be multiplied by the average shallow crust thermal conductivity of this region (1.3 W m -1 °C-1 ), to obtain an average heat flow value of 53 mW m -2 .Heat probe temperature measurements and geochemical analysis of core samples are consistent with these results.
For this particular site, there is evidence of mass wasting (Han et al. 2012(Han et al. , 2017)).The scale of this mass wasting features is about 1300 -2100 m in length, 300 -700 m in width, and 30 -150 m in thickness (Han et al. 2012).If such mass movements are active and wide-spread, temperature field within the shallow crust might not be in steady-state, depending on the ages of the mass movement events, and such processes might affect our results.However, the scale of regional temperature field model (3.85 km by 17.5 km) is much larger than these local mass wasting features.We have not considered these effects in our model as they do not cover the entire 3D seismic cube.Even though we have not been able to consider effects of such mass wasting, they may have local effects.
We propose that this method can be applied to other regions as suggested by Minshull and Keddie (2010), where seafloor temperatures are relatively constant in terms of temporal and spatial changes and where there is good velocity control in the depth interval between seafloor and the BSR.We recommend that topographic correction and depth-dependent thermal conductivity should be considered where possible.Through conducting sensitivity tests, we determined that there are limitations to the accuracy of gas composition using this estimation method.Higher impurities in the gas composition may be possible if capillary effects are significant (Fig. 6).Another issue which was not considered in this study is that we have used only hydrate structure I for calculating the phase boundary.
If a linear geothermal gradient, constant gas composition, and constant pore fluid salinity are assumed, then there is a large scatter as observed in the temperature profile calculated from the BSR depths in Fig. 5.The uncertainty from the linear regression approach is so high that it indicates that the geotherm is not linear, but changes with depth, so a higher order function was used (Fig. 5).A cubic or an exponential function is usually used to fit the depth-dependent porositydepth data (e.g., Sclater and Christie 1980), and thus the temperature profile for the shallower sub-bottom depths.The results from this study allow us to reduce the uncertainty in estimations of geothermal gradients arising from uncertain phase boundary conditions.The overall uncertainty is larger because of the uncertainties in the interval velocity between the seabed and the BSR, although these are unlikely to be greater than 5% in the BSR temperature as shown in sensitivity tests done for this region (Chi and Reed 2008).Here we have used an average 1D velocity model derived from detailed velocity analyses of a 1 km long streamer, which should provide relatively good controls on the shallow crust velocity.With the hydrate molecular composition derived from this study, we have derived an updated geothermal gradient map from this seismic dataset (Fig. 7).

ConClusIons
We have used available 3D seismic data, in-situ thermal probe measurements, 3D finite element thermal modeling, and depth-dependent thermal conductivity data from offshore southwest Taiwan to estimate the chemical composition of gas hydrate at Penghu Canyon.Based on this data, we demonstrate that the hydrate composition here is similar to the results of in-situ geochemical analyses of the pore fluids recovered from gravity core samples.An advantage of this improved method is that it provides a technique that can be used to derive the molecular composition of hydrate-bound gases at depths of a few hundred meters below seafloor, where typical piston and gravity coring systems cannot sample.A similar method could be applied to many other areas in comparable settings (e.g., other passive margins) where available 3D seismic data containing signs of a BSR is available, allowing researchers to better constrain the molecular composition of natural hydrates in the marine environment.
Fig. 1.(a) Location of the study area offshore Taiwan.The location of the 3D cube is marked by a rectangle.The line inside the rectangle shows the location of the seismic profile in Fig. 2. EP: Eurasian Plate; PSP: Philippine Sea Plate.(b) BSR sub-bottom depth.(c) Temperature at the BSR.

Fig. 2 .Fig. 3 .
Fig. 2. Example of an in-line seismic profile from the 3D seismic volume, the position of the BSR is indicated by the white arrows.
Fig. 4. (a) Temperature field of the study area modelled using the 3D finite element method.Based on this result we can estimate the effect of topography on temperatures derived at the sub-bottom depths of the BSR at different parts of the volume.(b) We then use this information to correct the seismic-derived BSR temperature before using the Minshull and Keddie method.This schematic figure shows the blue line as the location of a BSR.The BSR sub-bottom depth is greater under the ridge, but shallower under the flanks, due to the deflection of heat source from depth away from the ridge to the flanks.Such topographic effects can be corrected using the information from Fig.4a, allowing us to derive a consistent geothermal gradient, i.e., the same BSR sub-bottom depth for the whole BSR, as shown by the corrected line.

Fig. 6 .
Fig. 6.Sensitivity tests on the change in geothermal gradient and seafloor temperature with varying compositions of gas hydrates.For varying chemical compositions, the salinity is fixed at 2%.(a) Seafloor temperature (°C).(b) Geothermal gradient (°C km -1 ).