Gas hydrate and free gas estimation from seismic analysis offshore Chiloé island ( Chile ) *

In this study one seismic section offshore Chiloé Island was analyzed to better define the seismic character of the hydrate-bearing sediments. The velocity analysis was used to estimate the gas-phase concentration and relate it to the geological features. The velocity model allowed us to recognize two important layers that characterize hydrateand free gas-bearing sediments above and below the BSR respectively: one located above the BSR, characterized by high velocity (1,800-2,200 m/s) and a second one, below the BSR, characterized by low velocity (1,600-1,700 m/s). A weak ref lector at about 100 m below the BSR marks the base of the second layer. AVO analysis and offset stack sections confirming that the reflector interpreted as BGR is related to free gas presence in the pore space. The velocity field is affected by lateral variation, showing maximum (above the BSR) and minimum (below the BSR) values in the sector. Here, the highest gas hydrate and free gas concentrations were calculated, obtaining 9.5% and 0.5% of total volume respectively. A variable BSR depth (from 300 to 600 mbsf) can be justified supposing a variable geothermal gradient (from 25 to 45 °C/km).


Introduction
Gas hydrate occurrence in the natural environment has a global significance because of the potential energy resource represented by the large amount of hydrocarbon trapped in the hydrate phase (Milkov, 2004).The methane plays a role in global climate change (Kennett et al., 2003), and its accumulations are potential hazards which may cause damage to drilling and seabed installations (Hovland and Gudmestad, 2001).
During the last decades many efforts have been made to identify and quantify gas hydrate presence.Geophysical studies along the Chilean margin through seismic ref lection data analysis have given us the opportunity to recognize the Bottom Simulating Reflector (BSR) and to map its presence (Bangs et al., 1993;Brown et al., 1996;Díaz-Naveas, 1999;Morales, 2003;Vargas-Cordero et al., 2011;Villar-Muñoz et al., 2014).This ref lector, representing the acoustic limit between a gas hydrate layer above the BSR and free gas layer below the BSR, is extended along the Chilean margin between 33° to 47° S. In seismic sections, the BSR presence allows defining the boundary between the high seismic velocity (i.e., the gas hydrate-bearing sediments) and the low seismic velocity because of the free gas presence (e.g., Shipley et al., 1979;Hyndman and Spence 1992;Berndt et al., 2004).In literature, few studies related to the hydrate quantification along the Chilean margin are available (Rodrigo et al., 2009;Vargas-Cordero et al., 2010).Anyway, the results of these studies have allowed estimating the gas hydrate and free gas amounts in limited areas of Valdivia, Itata and Coyhaique offshore.
Both compressional and shear wave velocities provide information about the presence of gas hydrate and free gas in marine sediments (Tinivella and Accaino, 2000;Tinivella, 2002;Bünz et al., 2004).To interpret seismic data in terms of hydrate content, it is necessary to establish a relation between hydrate saturation in sediments and their velocity.Lee et al. (1996) developed a weighted equation model, in which hydrate is considered as part of the pore f luid.However, this model is not based on physical principles and requires empirical calibration.The velocities predicted from the cementation theory proposed by Dvorkin and Nur (1996) are much higher than those normally observed in nature (Ecker et al., 1998).Helgerud et al. (1999) introduced an effective medium theory that considers the hydrate as part of the rock frame and successfully applied their approach to the P-wave velocity data from Site 995 of ODP Leg 164 in the Blake Ridge area.Lee (2002) proposed a method based on Biot-Gassmann theory to relate the elastic properties of sediments to those of the matrix and the pore fluid.Finally, Chand et al. (2004) compared various models and found modest variations in the pattern of velocity variation with hydrate saturation if the same f luid and matrix properties are assumed.Here, it was adopted the method proposed by Tinivella (2002), based on the Biot-Gerstmann-Smith equation, successfully applied to direct measurements (ODP Leg 164; Tinivella and Lodolo, 2000) giving results in agreement with the three-phase Biot theory (Tinivella, 1999;Carcione and Tinivella, 2000).
Offshore Chilean margin and offshore Antarctic Peninsula important gas hydrate reservoirs have been detected in the last years (see for example, Vargas-Cordero et al., 2010, 2011and Loreto et al., 2011).These regions are very sensitive to climate change as pointed out by recent studies (i.e., Marín-Moreno et al., 2015), confirming that it is necessary to increase our knowledge about the quantity of gas trapped in marine sediments in order to include it in the global carbon cycle.In this context, we focused our attention to the southwest part of Chiloé Island across the continental slope between 43° and 44° S (Fig. 1).This area is characterized by subduction of Nazca Plate below the South American continental plate along the Perú-Chile trench at a rate of about 66 km Ma-1 (Angermann et al., 1999;Kendrick et al., 2003).Seismic line SO161-40, analyzed in this study, is located in the area where the last glaciations have modeled the Chilean southern margin (Rabbasa and Clapperton, 1990) and the tectonic evolution of the fore arc basin included an inversion from tectonic erosion to accretion (Melnick et al., 2006).These aspects are important to consider in the seismic analysis across the BSR, because tectonic processes and glaciations determine seismic velocity changes that must be distinguished from those induced by the distribution of hydrate and free gas in the sediment pore space.
The main goal of this research is to estimate the gas hydrate and free gas amounts along the seismic line SO161-40 in order to understand the potential energy source of this area, as well as the environmental influence of the gas-phase presence.
In fact, a quantification of gas presence in the marine sediments is important not only from the energetic point of view (even if the exploitation of hydrate in marine environment is still in the early stages of development).On the other hand, it is a priority to understand if the gas released by hydrate dissociation due to pressure-temperature change can dissolve in the sea water, consume by the ecosystem or potentially reach the atmosphere increasing the green-house effect.

Material and Methods
The adopted procedure includes: a) detection of gas hydrate and free gas presence by performing an advanced seismic processing; and b) estimation of gas hydrate and free gas amounts by using velocity.Then, the geothermal gradient was estimated from sea water and BSR depths, knowing the sea bottom temperature.This procedure has been partially tested along the Chilean margin (Vargas-Cordero et al., 2010).In addition, in this study Amplitude Versus Offset (AVO) analysis was included (Yilmaz, 2001).

Seismic data
In this study seismic data acquired during the RV Sonne cruise SO161 (January-February 2001) in the frame of the project "Subduction Processes off Chile (SPOC)" were analyzed.The acquisition parameters of Line SO161-40 are the following: an array of 20 air guns, tuned to suppress the bubble pulse with a total volume of 54.1 l a shot interval of 20 s, producing a shot spacing of approximately 50 m; and a streamer length of 3 km, consisting of 132 channels.From channel 1 to 24 the inter-trace was 12.5 m, while from channel 25 to 132 was 25 m.

Advanced processing and inversion modeling
To detect and quantify the gas hydrate and free gas presence a procedure already tested was adopted (Tinivella et al., 2009;Vargas-Cordero et al., 2010).Once recognized the BSR in the post-stack section, a part (18 km) of the line, where the BSR was strongest and continuous, was selected to perform an advanced processing.The target includes: a) Kirchhoff Prestack time migration (PreSTM) and b) a detailed velocity model using the layer stripping approach (i.e., Yilmaz, 2001).The layers are modeled in depth with an iterative approach by using a Kirchhoff pre-stack depth migration (PreSDM) (see details in Tinivella et al., 2009).
The initial velocity is assumed constant and equal to 1,480 m/s (water seismic velocity), assuming horizontal and vertical spacing of the grid equal to 25 and 10 m, respectively.After six iterations the sea-bottom reflection within CIGs was flat suggesting that a correct migration velocity was used.Then, the velocity analysis of a second layer (top: the seafloor; bottom: a horizon between the seafloor and the BSR, hereafter called horizon 1) was performed.After 30 iterations, the velocity for the second layer was fixed.A third layer (top: horizon 1; bottom: BSR) was updated after 45 iterations.Finally, the velocity in the free gas zone below the BSR and above the so-called Base of free Gas Reflector (BGR) was updated by 20 iterations.
After inserting a velocity gradient below the BGR, the final velocity model was smoothed because the PreSDM is more stable for smoothed models (Liu and Bleistein, 1995).In order to obtain an accurate image and attenuate the stretching effects, a stack of CIGs was produced considering a maximum offset of 2,500 m.Trace mixing and band-pass filter were applied to the final stacked section.

BSR-derived geothermal gradient
The geothermal gradient (dT/dZ) has been calculated with the following formula: where T BSR and T SEA are the temperatures at BSR and at sea floor, respectively, while Z BSR and Z SEA are the depths of the BSR and the seaf loor respectively (i.e., Grevemeyer et al., 2003).The BSR and the seaf loor depths were obtained from the PreSDM section, while the seaf loor temperature (equal to 2.2 °C) was taken from CTD measurements off Chile (Grevemeyer and Villinger, 2001).The estimate of the temperature at BSR is based on the dissociation temperature-pressure function of gas hydrates, considering a salinity of 35‰ (Sloan, 1998;Dickens and Quinby-Hunt, 1994).This study considered that only methane is present, which is the most restrictive case (i.e., shallower BSR).In fact, even if well data pointed out the presence of 1% of ethane (Froelich et al., 1995), Vargas-Cordero et al. (2010) suggested that the effect of the ethane on theoretical BSR depth is very slight and can be neglected.

Estimate of gas hydrate and free gas concentrations.
The method to estimate gas hydrate and free gas concentrations consist of comparing velocity anomalies with theoretical velocity curves in absence of gasphases.Because of the lack of direct measurements, a simplified theory, instead of more refined theory (i.e., Tinivella and Carcione, 2001;Chand et al., 2004), was used.In this case, the main error is due to the sediment property assumptions and not to the theoretical model adopted to determine the theoretical velocity versus hydrate/free gas concentration.Modified Hamilton's curves (Hamilton, 1976(Hamilton, , 1979) ) were used in order to reproduce the velocity field in absence of gas or full water saturation (Tinivella, 1999).
A qualitative estimate of concentrations can be obtained by comparing the theoretical velocity for full-water saturation to the seismic velocity, evaluated by the PreSDM velocity analysis.Positive anomalies indicate the presence of gas hydrates, whereas negative anomalies indicate the presence of free gas.In order to perform a quantitative estimation, the gas-phase concentrations are increased in the theoretical formula starting from 0. The process is stopped when the theoretical velocity fits the velocity obtained by PreSDM.The method can model two main distributions to calculate the concentrations of free gas in the pore space: a) uniform distribution (gas and water in pore space), and b) patchy distribution (all gas in patches without water).In the present study, a uniform distribution of free gas in pore space was considered, due to the low velocity observed in the free gas layer (less than 1,600 m/s; Tinivella, 2002).Missing direct measurements, the Hamilton trend for the porosity/density versus depth (Hamilton, 1979) was considered, based on a 50% porosity measured in a ODP well of the same region (Grevemeyer and Villinger, 2001).

Amplitude versus offset curves and amplitudes change with offset
The data processing was focused on the signalto-noise ratio increase and vertical resolution enhancement.A "preserving-amplitude" processing was adopted in order to allow a successive AVO analysis (Yilmaz, 2001).In fact, AVO analysis was performed in order to have useful information about the petro-physical properties of the shallow layers, such as the fluid content (i.e., Giustiniani et al., 2009).
The AVO effects are evident in the stacked sections obtained by using different offset ranges: lower than 500 m (hereafter called near-offset stacked section); from 500 to 1,700 m (hereafter called medium-offset stacked section); and higher than 1,700 m (hereafter called far-offset stacked section).This choice was performed evaluating the angle of incidence versus offset by using the velocity model obtained by seismic data analysis as described in the previous paragraphs.In fact, the near offset corresponds to a range of incidence angle from 0° to 10°, the medium offset from 10° to 30° and the far offset means incidence angles greater than 30°.

PreSTM section
The PreSTM section (Fig. 2) allows recognizing the following features: 1. Chaotic and low amplitude reflectors of the upper continental slope (from CDP 100 to 300).
2. Landward sub horizontal and discontinuous reflectors characterized by normal faults in the upper continental slope (CDP 300 to 600).Below these sediments, a high amplitude reflector with evident folding about 1 s from the seafloor is recognized.Note that a weak and discontinuous BSR is detectable.

Velocity model and PreSDM sections
The final velocity model shows the layer above the BSR characterized by high velocity (1,800-2,200 m/s), with respect to the normal compacted marine sediments (Fig. 3), which can be associated to gas hydrate presence.A layer below the BSR is characterized by low velocity (1,600-1,700 m/s) and can be associated to free gas presence (Fig. 3).A strong lateral velocity variation is evident along the section.In particular, a low velocity across the BSR is reached upward of upper continental slope (km 12-13), while a high velocity above the BSR is present downward of it (km 4-6).
The BSR is present from 4 to 16 km.It presents a gap where the area is affected by intense deformations, characterized by faults and fractures that cut sediments going up to the seafloor.Moreover, a variable BSR depth was recognized reaching maximum value in correspondence to the boundary between the lower and upper continental slope (Fig. 3; 600 mbsf -meters below sea floor), whereas the minimum value was identified in correspondence to the uppermost part of continental slope (300 mbsf).The PreSDM section reveals possible gas hydrate and free gas thickness of about 300 m and 100 m, respectively.The amplitude of the BSR is variable along the section and inversely proportional to the velocity, as expected (Fig. 3).
The BGR is present along the whole seismic section indicating the base of a layer of low velocity, which can be associated to the free gas presence.As evidenced in the blow-up in figure 3, the BGR presents a normal polarity with respect to the BSR, as expected.Note the loss of ref lectivity below the BSR that could be due to the free gas presence.

BSR-derived geothermal gradient
A variable geothermal gradient along the seismic section is recognized, with a maximum value equal to 45 °C/km, located in correspondence to the uppermost part of the continental slope, and a minimum value equal to 25 °C/km westwards of it (Fig. 4).Note that the highest value is in proximity of the border line and can be affected by a higher uncertainty than the values estimated in the central part of the seismic line.The theoretical BSR depth is evaluated considering that the hydrate contains only methane, as discussed in the previous paragraphs.

Estimate of gas hydrate and free gas concentration and AVO curves
The results of the gas-phase estimate indicate high variability along the section.The highest gas hydrate concentration is located between 4 and 10 km (about 9.5% of total volume; Fig. 4), while the highest free gas concentrations (0.5% of total volume) are detected where highest gas hydrates concentrations are present (Fig. 4).An average of gas hydrate and free gas concentrations is obtained, 6% and 0.3%, respectively.The velocity field and the related gas hydrate concentration are well characterized by the analysis described above.On the contrary, the free gas layer is difficult to analyze, because the intrinsic characteristics of the free-gas bearing sediments (i.e., the drop velocity for very small content of free gas).The AVO analysis is widely used in literature as a tool in gas-phase identification; so, in order to detect AVO effects due to gas presence, we obtained three stacked sections by using selected offset range: near, medium and far offset (Fig. 5).Then, we evaluated the theoretical PP reflection coefficient at the BGR, following the method described in Carcione and Tinivella (2000).The adopted values are reported in Table 1 and the result is shown in figure 6.As can be observed, the PP ref lection coefficient is quite constant for near offset and increases slightly for medium offset.On the contrary, the PP ref lection coefficient increases rapidly at far offset.The same trend is observed in the offset-stacked sections, as reported in figure 5, confirming that the ref lector interpreted as BGR is effectively related to the free gas presence in the pore space.For example, this effect is evident between CDP 3,200 and 3,400, as shown in figure 5.

Discussion and conclusions
The seismic velocity analysis indicates strong lateral variations in the seismic section.Highest and lowest velocities are observed only where the BSR is present (Fig. 3).Note that the reflector selected below the BSR can only be locally associated to the BGR.As evidenced from seismic profiles (Figs. 2 and 3), faults and fractures show small slips affecting shallowest sediments.Moreover, small deformations within sediments reach the seaf loor, thus we infer the presence of active faults.The occurrence of active faults is related mainly to lateral variation and secondly to depth variation of the BSR.Generally, along active faults intense fluids up-welling can be guided and accumulated within a trap, resulting in temperature increase of sediments.In fact, a variable BSR depth (from 300 to 600 mbsf) can be justified supposing a variable geothermal gradient (from 25 to about 45 °C/km; Fig. 4).The estimated geothermal gradient values are in agreement with values reported in literature (Bangs et al., 1993;Brown et al., 1996;Villar-Muñoz et al., 2014).The maximum depth of the BSR, reached in the down side of upper continental slope (600 mbsf), is associated to the decrease of geothermal gradient from 45 to 25 °C/km.This decrease could be justified by a combination between reduced f luid circulation and increased sedimentation rates, or with morphology effect, as can happen along the boundary between the lower and upper continental slope, which is characterized by a dip change (i.e.,  VP: compressional-wave (P) velocity.VS: shear-wave (S) velocity.ρ: density.QP: quality factor of P-wave velocity.QS: quality factor of S-wave velocity.BGR: base of the free gas layer.The values of velocities and densities are in agreement with the theoretical model adopted in this study.The values of QP and Qs are in agreement with Tinivella and Accaino (2000), because the lithological and tectonic characteristics in the South Shetland margin are similar to our studying area.Ruppel, 1997;Ganguly et al., 2000).Clearly, the estimation of the geothermal gradient can be used for other purposes, such as regional studies including land region.The highest velocities can be related to two main factors: 1) the gas hydrate presence and 2) the change in petro-physical properties because of the different marine sediments compaction, which can be explained assuming high dewatering (Loreto et al., 2007) due to fault presence at the border of the upper slope (Fig. 2).Moreover, continental glaciations effects (Rabassa and Clapperton, 1990) may cause the porosity reduction and the consequent velocity increase as reported by Accaino et al. (2005).Similar conclusions may be drawn analyzing the estimate of the gas-phase concentrations, derived from the velocity model.In fact, the high velocity above the BSR means high gas hydrate concentrations and the low velocity below the BSR observed in uppermost side of continental slope means high free gas concentration.It is important to underline that, in absence of direct measurements; the estimate is affected by inaccuracy: the main reliable result is the distribution of both hydrate and free gas in the pore space along the margin.Nevertheless, this estimate is consistent with the values obtained by other authors analyzing different data in the same area (Bangs et al., 1993;Froelich et al., 1995;Brown et al., 1996;Vargas-Cordero et al., 2010).For example, Brown et al. (1996) estimated that the gas hydrate concentration ranges from 5 to 20% of the pore volume, while the free gas concentration ranges from 1 to 3% of volume.In any case, this procedure can be used to determine variability at regional scale, as confirmed by the order of magnitude of our estimate.Indeed, regional and local studies (e.g., Rodrigo et al., 2009;Vargas-Cordero et al., 2009, 2010, 2011) suggest that gas hydrates are present in accretionary complexes in subduction zones, reaching 15-20% of total volume in frontal accretion, and 3-5% of total volume in basal accretion systems.
This study suggests that gas hydrate can play an important role in this part of the Chilean margin.In fact, the high local concentrations of both hydrate and free gas, as suggested by several studies, could be considered as a potential future energy resource.
FIG. 2. Pre-stack time migrated (top) and line drawing (bottom) sections, in which the main features are indicated.

FIG. 3 .
FIG. 3. Top: final velocity model superimposed to the pre-stack depth migrated section.Bottom: zoom of the selected parts, where the BSR and BGR are clearly detected.

FIG. 4 .
FIG. 4. Top: BSR-derived geothermal gradient.Bottom: Gas phase's concentration model superimposed to the pre-stack depth migrated section.Positive values are referred to the gas hydrate amount and negative values are referred to free gas amount.