# Model-Based Identification of the Base of the Gas Hydrate Stability Zone in Multichannel Reflection Seismic Data, Offshore Costa Rica

Authors: Romina Gehrmann; Christian Müller; Peter Schikowsky; Thomas Henke; Michael Schnabel; Christian Bönnemann
Journal: International Journal of Geophysics
 (2009)

---

## Abstract

Along the pacific margin offshore Costa Rica the Bottom Simulating Reflector (BSR) shows a patchy occurrence in 2-D seismic reflection profiles. The reason for this can be either lack of free gas beneath parts of the gas hydrate stability zone (GHSZ) or poor seismic imaging. We compare far to near offset stacked common midpoint sections to reduce imaging ambiguity utilizing the amplitude variation with offset effect and thus successfully distinguish BSRs from regular sediment reflections. In combination with 1-D modeling of the base of the GHSZ we disqualify or qualify reflections in the predicted depth range as BSR. Additionally we calculate the heat flow and compare it with an analytical solution to detect thermal anomalies, for example, at the frontal prism. The higher confidence in BSR depths based on the far offset stacks and heat flow calculations allows further analyses on gas hydrate concentration estimates and tectonic evolution of the margin.


---

## Body

## 1. Introduction

The Bottom Simulating Reflector (BSR) that has been recognized on continental margins worldwide (e.g., [1–6]) is caused by a negative impedance contrast at the base of the gas hydrate stability zone (GHSZ). Either high amounts of gas hydrate above and/or a few percent of free gas below the GHSZ produces this strong phase-reversed reflection [7, 8]. On a reflection seismic profile the BSR is often a continuous reflector that mimics the seafloor and crosscuts stratigraphy at depths of a few hundred meters below seafloor (mbsf) [9]. In our work we focus on the amplitude variation with offset that is generally known as the AVO effect [10], as the magnitude of the BSR reflection amplitude increases with increasing angle of incidence if the reduction of Poisson’s ratio across the BSR is sufficient. That holds in the presence of free gas below or high amounts of gas hydrate above the BSR [11]. The AVO effect has been studied quantitatively with high effort on amplitude recovery and forward modeling of the reflection coefficient of the BSR in order to determine the physical properties of sediments containing hydrates, liquids, or free gas [12–14]. Offshore Costa Rica the amount of gas hydrate is considered to be too small to cause the AVO effect observed in common midpoint gathers. The observed AVO effect is rather reproduced assuming a small amount of free gas beneath the gas hydrate stability zone [13].The BSR has been mapped offshore western Costa Rica on 2-D reflection seismic profiles during cruise BGR99 in 1999 [13]. According to this interpretation (Figure 1), the BSR is not continuous over the entire margin. The aim of our study is to analyse the BSR by comparing near and far offsets stacks where it is rather weak at zero incident angles to constrain the depth of the base of the GHSZ. With this approach we consider the AVO behaviour of the BSR qualitatively.Figure 1
The reflection seismic profiles (light grey lines) offshore western Costa Rica were acquired during cruises BGR99 and SO81. The crosses mark the drill sites of ODP Leg 170 north-west of the Paleo Plate Boundary (PPB). The PPB refers to the boundary of the Cocos Plate’s two distinct origins, the East Pacific Rise, and the Cocos Nasca spreading centre. The Ridge Jump is characterized by a seamount chain [32]. The thick black line is profile BGR99-46, which is discussed in this work and is located offshore Nicoya Peninsula. The figure is modified from Müller et al. [13] who interpreted the BSR occurrence (dark grey shaded area) along the continental margin from onboard processed data.In order to predict a depth range for the BSR in the far offset stacks we estimate the base of the gas hydrate stability zone. The stability of gas hydrate at continental margins is controlled by pressure and temperature conditions as well as the composition of the gas and seawater salinity [15]. The heat flow trend [16–18] is directly related to the depth of the gas hydrate stability zone which is strongly temperature dependent. Therefore also the BSR seems to follow an isotherm beneath the seafloor as long as the water depth is relatively constant [8]. The BSR has been used to calibrate the base of the gas hydrate stability zone and to estimate the heat flow at continental margins [19–21]. In this study we calibrate the base of the gas hydrate stability zone at clear BSR reflections and calculate the geothermal gradient using the gas hydrate stability curve of Tishchenko et al. [22], which, compared to others, includes a variable salinity. We interpolate the geothermal gradient across areas with no clearly visible BSR and derive potential BSR depths in order to identify weak BSR reflections that are difficult to distinguish from sediment reflections otherwise. We compare the heat flow trend to in-situ drilling data and to trends predicted from thermal modeling to constrain the thermal regime of the continental slope.The thermal regime may be influenced by fluids from the subducting plate escaping at mud mounds at the middle slope [23, 24] and by recent slumping or erosion [25]. Thus, the estimated heat flow offers an insight into the active geological development of the margin.A notable difference between the heat flows of the area north-west of Nicoya Peninsula and the area to the south-east (towards Osa Peninsula) is related to the different origins of the subducting oceanic Cocos Plate [26, 27]. The oceanic plate subducting north-west of Nicoya Peninsula originates at the East Pacific Rise spreading centre and has a smooth surface. It is hydrothermally cooled via deep faults developed during plate convergence-related bending at the Middle American Trench [26, 28]. The south-eastern part of the subducting Cocos Plate is formed at the Cocos Nazca spreading centre, and its morphology is characterized by seamounts and the Cocos Ridge, which controls slope morphology [29], slope failure as, for example, the Nicoya Slide [30] (Figure 2), and basal erosion of the overriding continental plate [31]. The Paleo Plate Boundary between these two segments has been imaged via magnetic, bathymetric, and seismic studies [32] and trends towards the toe of Nicoya Peninsula. Across the plate boundary close to the Middle American Trench, Fisher et al. [27] detected an abrupt change in heat flow from 20–40 mW/m2 (East Pacific Rise) to 105–115 mW/m2 (Cocos Nazca spreading centre).Figure 2
Bathymetry close-up [64] of the area around line BGR99-46, the ODP Leg 170, the heat flow measurement locations S6 [18] and M54-43/113 [16], the Nicoya Slide, the Fisher Seamount, white encircled mud mounds [60], CMP 5690 (Figure 4), and the projected CMP 1410 (black dot) of profile SO81-09, where Pecher et al. [9] have applied a full waveform inversion (discussed in the text).Our work focuses on seismic line BGR99-46 that is located to the south-east of the plate boundary, north-west of the Nicoya Slide, and strikes perpendicular to the trench axis (Figure2). This line shows all typical characteristics of the continental margin of western Costa Rica (Figure 3). This includes its non-accretional frontal prism with mainly terrigenous sediments [33] that is influenced by compression [34], its extensional middle slope [35], the ophiolitic margin wedge with its rough surface [5], and sediment failures where no BSR is observed. The appearance of the BSR is thereby influenced by small-scale tectonics [9] and is ambiguous or even absent in some parts of the profile. In the following we show that seismic data processing with focus on comparing far to near offset CMP stacks and thus accounting for the AVO effect of the BSR qualitatively, a predictive 1-D forward modeling of the base of the gas hydrate stability zone, and the heat flow calculation over the entire profile improve the identification of the BSR against sediment reflections. This method leads to a stronger constraint on the heat flow profile across the continental slope and the thermal regime in the slope sediments. The resulting thickness of the GHSZ is an important input parameter for the quantification of gas hydrate inventories using non-seismic methods [36].Figure 3
Seismic profile BGR99-46. A BSR is patchy distributed (black bars) from the Middle American Trench (MAT) to its outcrop at~600 mbsl. The oceanic crust is characterized by a horst and graben structure. The margin wedge is strongly disrupted. The large slide mass, the surface erosion, and the mud mound are also marked on the bathymetry shown in Figure 2. The dashed box is the close-up shown in Figure 6.


## 2. Seismic Data and Processing

The seismic reflection data have been acquired offshore Costa Rica on cruise BGR99 in November 1999 by the Federal Institute for Geosciences and Natural Resources. The data acquisition parameters are shown in Table1. The long source-receiver offset of about 6 km and a high signal-to-noise ratio are appropriate preconditions to apply AVO analysis [13]. The mean slope angle on this profile lies between 1° and 4.5° (the maximum slope angle reaches 6°) resulting in a CMP smearing that is well within the Fresnel zone of ~370 m. The smearing effect has therefore not been considered within this study. The applied seismic data processing flow is shown in Table 2.Table 1
The data acquisition parameters of reflection seismic profile BGR99-46.
Data acquisition parametersCMP spacing6.25 mShot distance25 mFold120Number of airguns40Total airgun volume62.3 lAirgun array length15.56 mReceiver group length12.5 mMean frequency~50 HzNumber of channels480Maximum offset6099.3 mRecording length8.7 sSampling rate4 msTable 2
The processing flow applied to profile BGR99-46 is divided into two parts. The offset stacks are created after the application of a band-pass filter and amplitude recovery, while the full CMP stack is created after an additional deconvolution and automatic gain control (AGC). Finally, the profile is either migrated in depth (for retrieving BSR depths) or migrated in time (for BSR interpretation).
Processing stepsParametersOrmsby band-pass filter3-6-60-120 HzNormal Moveout (NMO) velocity analysisEvery 100 CMP (~625 m)True amplitude recovery (TAR): spherical divergence correctiong(t)=v2(t)t*/v02t0* [37]Offset Amplitude Recovery (first order correction of source and receiver array geometry)Mean frequency, RMS velocity function, gun array length, receiver group length (as listed above in Table1)Offset CMP stacksEvery 50 channels, offset range: 625 mPredictive deconvolutionOperator length: 50 ms, predictive distance: 16 msNMO correctionNMO stretch mute values: 30%, 100%Automatic gain control (AGC)Time window: 500 msCMP stackMeanKirchhoff time/depth migrationMax. dip:25°, max. frequency: 80 Hz, smoothed velocity fieldThe NMO velocity analysis can generally not resolve the gas layer beneath the BSR due to the vertical resolution limits at the depth of the BSR of around 13 m. However, the velocity increase with depth is reduced within the depth interval underneath the BSR. The typical velocity of this depth interval varies between 1750 and 1850 m/s. For comparison, a full-waveform inversion by Pecher et al. [9] on seismic profile SO81-9 (Figure 1 and CMP projected on Figure 2) resulted in a P-wave velocity reduction from 2000 m/s to 1720–1760 m/s for a 6–18 m thick gas layer beneath the gas hydrate stability zone. This low-velocity layer may create a thin bed tuning effect and artificial amplitude variations of the BSR reflection [38] but is not necessarily to be expected, as the base of the free gas layer is rather gradual [39]. We do not image the bottom of the free gas layer as a separate reflection, and our AVO response of the BSR does not show tuning-related amplitude distortions.The predictive deconvolution suppresses reverberations and compresses the wavelet [10]. In the case of an NMO stretch mute value of 100% the deconvolution is especially helpful to assure a clear image since the moveout correction has a similar effect to a low-pass filter by stretching the traces in a time varying manner. Frequency distortion especially increases at large offsets and shallow times [10]. Two different kinds of NMO stretch mute values are chosen for the CMP stack. For the BSR identification, signal periods with doubled length are preserved (NMO stretch mute 100%). In that way far offsets are included, thus providing a stronger image of the BSR (Figure 4). For interpretation of the seismic sections, a standard NMO mute of 30% is applied to ensure high vertical resolution. An AGC is applied before the full CMP stack in order to amplify the BSR and deeper structures in the migrated section. The section is then depth-migrated in order to obtain BSR depths from clear BSRs for calibration of the predictive model.Figure 4
CMP gather 5690 is an example for the AVO effect of the BSR. The black lines indicate the cutoffs for NMO stretch mute values of 30% and 100%. A late cutoff (NMO stretch mute 100%) includes far offsets and therefore stronger BSR amplitudes after CMP stacking. The upper box shows the absolute values of the amplitudes for the seafloor reflection (red dots) and for the BSR (blue dots).


## 3. The AVO Effect at the BSR on Far Offset CMP Stacks

The BSR at the continental margin of Costa Rica is characterized by an AVO effect that is caused by the presence of a few percent of free gas beneath the gas hydrate stability zone [13]. While the presence of gas hydrate increases P-wave velocity VP, free gas beneath the GHSZ reduces P-wave velocity and density [11]. The shear wave velocity VS, however, does not change significantly if only the water within the pores is substituted by free gas. Consequently, Poisson’s ratio that is related to the ratio of VP to VS drops from about 0.42 for marine hydrate-saturated sediments to nearly 0.1 for free-gas-saturated sediments [11]. This combination causes the BSR reflection to be phase reversed compared to the seafloor reflection and to increase in amplitude with increasing angle of incidence. The inversion of the reflection amplitude response for elastic parameters (e.g., the Poisson’s ratio) is non-unique [14], but for small incident angles a rule of thumb says that free gas induces a decrease of the reflection coefficient [40]. A few percent of free gas beneath the GHSZ causes a significant AVO effect [13], while higher concentrations of gas do not change this effect significantly [41]. CMP gather 5690 (Figure 4) shows the increase of the BSR amplitude with increasing offset while the seafloor reflection amplitude decreases towards zero while approaching the critical angle. This behaviour of increasing absolute values of reflection coefficients is also observed at gas-sand layering [42] and AVO analysis is thus a method commonly applied in hydrocarbon exploration [43]. For an AVO analysis it is generally necessary to restore amplitudes for source/receiver directivity, transmission, and absorption loss [13]. We applied a standard processing sequence correcting for source and receiver directivity provided by the seismic processing software. By stacking the far offsets the AVO effect of the BSR is examined here in an efficient way without applying full quantitative AVO analysis.Figure5 shows the magnitude of the theoretical P-wave reflection coefficient versus incident angle of the seafloor and the BSR. The absolute values of the reflection coefficients are calculated with the Aki and Richards [44] approximation of the non-linear Zoeppritz equations for P-waves.Figure 5
The magnitudes of the reflection coefficients for the seafloor reflection and the BSR are calculated for a model with~3% free gas below the gas hydrate stability zone [44]. Density and seismic velocities are taken from the full-waveform inversion of Pecher et al. [9]. The black boxes mark the incident angle ranges which are shown in Figure 6.Figure6 shows a section of the time-migrated seismic profile BGR99-46 (Figure 3) with an NMO mute value of 100% used for CMP stacking. Two close-ups reveal the benefits of comparing far to near offset CMP stacks. The insets on the left (1a and 2a) show CMP stacks of the first 50 channels (near offset stack), 141–766 m distance from the source. This offset range is related to incident angles at the BSR depth of about 1.5°–9.7°, estimated from a two-layer earth model. They provide a sharp image of the sediment layers, but the BSR amplitude is weaker than sediment reflections. In the far offset stacks on the right (1b and 2b) (2641–3266 m offset and incident angles of about 25.9°–37.4°) the seismic image loses its high resolution caused by the NMO stretch effect, the subsequent longer signal period, and interference with side reflections. These far offset CMP stacks though provide a stronger BSR. The BSR can be more easily distinguished from horizontally layered sediments.Figure 6
Centre: the time-migrated section between CMP 5300–9100 shows the middle slope of profile BGR99-46. Top and Bottom: close ups 1a and 2a are the near offset CMP stacks (141 m–766 m) with incident angles of~1.5°–9.7°. Close ups 1b and 2b are the far offset CMP stacks (2641 m–3266 m) with incident angles of ~25.9°–37.4°. The colorscale for the offset stacks is modified to emphasize negative reflections. The areas of interest are marked with black ellipses.The far offset CMP stacks therefore enhance weak (but continuous) BSR reflections, but weaken reflectors that result from the normal velocity increase with depth related to dewatering and compaction in water-bearing sediment layers. They also yield new pieces of BSR reflections that previously could not be identified. The far offset stacks alone do not give unambiguous evidence of the base of the GHSZ. An AVO effect can also be caused by other gas or fluid accumulations within the sedimentary column. To ensure that an identified anomaly is at a depth that corresponds to the base of the GHSZ, we additionally implement a 1-D forward modeling approach to estimate the depth of the base of the GHSZ from local temperature and pressure conditions.


## 4.1-D Forward Modeling of the Base of the GHSZ

The base of the GHSZ is presumed to be at BSR depth in reflection seismic profiles [6]. Gas hydrate in the sediment column of continental margins stabilizes due to high pressures and moderate temperatures. The temperatures required for gas hydrate destabilization are higher when additional components of CO2, H2S, and higher hydrocarbons are present in the gas phase, while a higher salinity has the opposite effect [4, 15, 21, 45]. Site 1041 from ODP Leg 170 offshore Nicoya Peninsula revealed a heterogeneous gas hydrate distribution [33]. Average concentrations of 1.64 vol.% of pore space were calculated by Hensen and Wallmann [46] from the sediment ability to produce methane from organic matter. The resultant gas hydrate is mostly composed of biogenic methane with locally varying portions of thermogenic gases (e.g., at mud mounds). In this work the gas hydrate stability curve of Tishchenko et al. [22] is used, which is based on an empirical algorithm for the stability of methane hydrate in seawater with a variable salinity. It deviates only slightly from approaches by Miles [47] or Brown et al. [4]. For Pacific bottom waters offshore Costa Rica a salinity of 3.46% is chosen [48]. The gas hydrate stability curve of Tishchenko et al. [22] agrees well with the laboratory data of Dickens and Quinby-Hunt [45] and with ODP drilling samples analysed by Grevemeyer and Villinger [20] as shown in Figure 7. The transformation of pressure to depth is done assuming hydrostatic conditions [47].Figure 7
Methane hydrate stability curve of Handa [65] for pure water and Tishchenko et al. [22] for bottom water with a salinity of 3.46%, laboratory data of Dickens and Quinby-Hunt [45] and ODP data [20].Our workflow is as follows: BSR depths at clear BSRs are taken from the depth-migrated section that are shown by dark blue lines in Figure8(a). A geothermal gradient for the upper sediment layers is calculated from the ratio of the temperature difference between the seafloor temperature, and the BSR temperature to the thickness of the upper sediment layers between the seafloor and the BSR.(a) The clear BSR (blue lines) is regularly distributed along the profile. After calculating its corresponding geothermal gradient and interpolating it, a modeled BSR (red lines) is derived. (b) The resulting smoothed heat flow trends with error bounds and an analytical heat flow curve from Molnar and England [51] are compared to in-situ data (Figure 3). The values north of the Paleo Plate Boundary near ODP Leg 170 and south of the plate boundary (S6) were derived from Langseth and Silver [18] and the values south of the plate boundary M54-43 and M54-113 were acquired by Grevemeyer et al. [16].
(a)(b)The temperature at the seafloor is taken from a water temperature curve that is fitted to seafloor temperatures at ODP Leg 170 drill sites [33] and to CTD measurements from cruise SO144 [49]. The temperature at the BSR is taken from the stability curve at the hydrostatic pressure that equates to the BSR depth.The geothermal gradient is then linearly interpolated along the seismic profile into areas of ambiguous or absent BSRs. Now potential BSR depths are calculated using the interpolated geothermal gradients, seafloor temperature, and the gas hydrate stability curve. They are marked red in Figure8(a). The observed and the modeled BSR depths are then plotted into the depth-migrated seismic profile as blue and red lines, respectively (Figure 9).(a) The clear BSR (blue thick line) and predicted BSR (red line) are plotted into the depth-migrated section of line BGR99-46. (b) The turquoise lines are former BSR interpretations for which heat flow values were calculated. The left-most turquoise reflection can be disqualified as a BSR because it exceeds heat flow boundaries and contradicts its trend. The heat flow in the middle belongs to the actual BSR identified in the far offset CMP stacks. The reflection marked turquoise on the right belongs to a sediment unconformity.
(a)(b)The heat flowQ is derived from the smoothed geothermal gradient and an extrapolated thermal conductivity of 1 W/mK [33] and is shown in Figure 8(b).Its trend is similar to heat flow data south of the Paleo Plate Boundary (Figure3), whereas the values cannot be directly compared as the heat flow at the plate boundary is slightly disturbed [16]. Fekete [50] has mapped an abrupt change in BSR depth on a crossline that results in a heat flow change from ~40 mW/m2 to ~75 mW/m2 16 km landward from the trench axis. It is also compared to an analytical curve from Molnar and England [51] with a heat flow of 110 mW/m2 for the subducting oceanic plate as measured by Fisher et al. [27] and calculated for a 22.7 Ma old plate by Stein and Stein [52]. The model has previously been applied by Grevemeyer et al. [53] to approximate shear stress and the seismogenic zone at the Nazca plate subduction zone. The discrepancy of the heat flow at distances of more than 30 km from the trench is explained by the increase of the subduction angle from 6° to 13° [54], which is not considered in the analytical model but would cause a reduction of the calculated heat flow.The 1-D modeling of the BSR results in a BSR depth distribution, its related geothermal gradient, and the heat flow trend for the entire seismic profile. The calculated BSR depth and its upper and lower boundaries offer the opportunity to fill BSR gaps by identifying pieces of high reflectivity in the far offset CMP stacks and confirming ambiguous reflections as BSR. Furthermore, BSR-like reflections may be reconsidered by comparing their associated heat flow to the calculated trend. In any case, the modeling predicts the base of the gas hydrate stability zone and may then, for example, be used for gas hydrate quantification using non-seismic methods that require the base of the GHSZ as input parameter. The derived heat flow trend can be further compared to analytical and numerical models for understanding the temperature field of the subducting slab by implying parameters like the shear stress [51, 53] and by discussing thermal anomalies that indicate fluid migration pathways [24].


## 5. Interpretation of Line BGR99-46

The tectonic situation offshore Costa Rica influences the BSR distribution in a couple of ways. Profile BGR99-46 is a representative example for the diverse BSR occurrence as shown in Figure6 at the slope basin between CMP 5300 and CMP 9100.Between the anticline and the basin (CMP 5800–7000) the BSR is interrupted. Basal erosion triggered by a subducting seamount may have caused the basin to subside and the GHSZ to deepen. If the subsidence happened recently so that the pore water is still unsaturated or the TOC value is too low to produce enough methane, a strong BSR may not be present. A phase-reversed reflector (marked with a turquoise line on the left of Figure9 bottom) can be disqualified as a BSR as it does not follow the calculated heat flow trend and even exceeds the upper limit of the uncertainties.Beneath an eroded seafloor the BSR is ambiguous as there are two parallel reflections (Figure6, 2a and 2b). They are more distinct in the far offset stack, but the true BSR cannot be determined without looking at the model that agrees in depth with the deeper reflector (Figure 9). The upper reflector may be a paleo-BSR that is preserved because the erosion has taken place recently and there is still free gas present. It takes a few ten thousands of years to reach new temperature equilibrium down to BSR depth as calculated by finite differences [55] or by the characteristic thermal diffusion distance [56] for 50 m of seafloor erosion. Foucher et al. [57] discuss an upward shift of the base of the gas hydrate stability zone due to seafloor warming or tectonic uplift and a leftover BSR that may keep its reflectivity up to 104 years due to free gas that slowly diffuses upwards. Considering a downward shift, we favour the interpretation that the former base of the gas hydrate has altered the sediments diagenetically and kept an impedance contrast and thus this reflection [3].Further landward along the slope, the BSR is interrupted again within a slide mass (Figure6: CMP 8500). A deep seated continuous reflector with a positive AVO effect can be identified as a discontinuity, marked as a turquoise line on the right in Figure 9, rather than as a BSR because the modeled BSR is ~150 m above that reflection. On the far offset stacks the reflector continuation underneath the slope basin can be observed which also contradicts a possible BSR.With respect to the heat flow, some other results are obtained (Figure8). The BSR-derived heat flow at the frontal prism is higher than the assumed heat flow of 110 mW/m2. It may indicate a heated prism, for example, by warm fluids that migrate upwards from greater depths or may be explained by the heat flow uncertainties. Other studies that support the primary assumption are summarized in the following. The frontal prism of the continental margin is non-accretionary as concluded by Kimura et al. [33]. A great amount of fluids is subducted. Resulting overpressure triggers a fracturing of the upper crust that provides fluid pathways and induces basal erosion [58]. Silver et al. [28] have analysed fluid seeps at the frontal prism and propose their origin in 10–15 km depth. Furthermore, Silver et al. [59] explain elevated heat flow at the frontal prism by diffuse fluid flow through at least the outer 3 km of the prism.The horizontally stressed frontal prism [34] forms an anticline (around CMP 5780) which is followed landward by the extensional slope basin. A mud mound (CMP 5420) that may provide active fluid seeps from the subducting slab [60] seems not to be active anymore, or upwelling fluids are of the same temperature as the surrounding sediments. The BSR and its affiliated heat flow are undisturbed.Within the huge slide mass (Figure3: CMP 8600–10800 and Figure 8:  km 30–45) the absence of a clear BSR can be explained by the disruption of the sediments with gas exposal and gas migration towards the upper slope. A significant part of the seismic energy is diffracted. Several BSR-like reflectors in this area belong to the base of slide masses where fluids accumulated. Anyway, there can be small amounts of gas hydrate present as the sediments are similar to those 6 km further upslope where a clear BSR indicates the existence of gas hydrate.The outcrop of the BSR (Figure3) in ~600 m water depth is characterized by a slump. Müller et al. [13] also detected slumps resulting from the outcrop of the BSR at the seafloor, probably connected with gas hydrate destabilisation and free gas depletion.


## 6. Discussion of the 1-D Modeling and the Heat Flow Calculation

A BSR depth uncertainty in the forward modeling of 35 m (~12%) is retrieved from the geothermal gradient’s undulation that is calculated from mapped BSR depths with the gas hydrate stability curve from Tishchenko et al. [22]. This uncertainty is introduced by the model’s seafloor dependency. It follows all small-scale undulations whereas the temperature field in greater depth and the subsequent BSR are considerably smoother. The uncertainty calculation does not take into account the error of the mapped BSR depth as it is the uncertainty of the BSR depth that is calculated from interpolated geothermal gradients and calibrated on mapped BSR depths. It is therefore the uncertainty that constrains the depth range in which we look for the BSR in the far offset stacks.The model assumes a conductive vertical heat transport. It does not include a horizontal component. This becomes relevant at steeper seafloor morphology and decreases the accuracy of the predicted BSR depth. An example is the modeled BSR in Figure9(a) between CMP 6500 and 7000 where the uncertainty to the BSR retrieved from far offset CMP stacks (turquoise line) reaches 50 m.Between km 30 to 45 in Figure8 the model is neither calibrated by any clear BSR nor proved by patchy reflections in the far offset CMP stacks and probably exceeds an uncertainty of 50 m.The model is based on 1-D heat conduction in a 3-D environment and a linear interpolation of the geothermal gradient. To interpolate the gradient linearly is a good approximation over short distances because the heat flow curve changes smoothly. It is however not appropriate in cases of thermal anomalies and greater distances. A regular distribution of clear BSRs along the profile or heat flow measurements for model calibration would keep the uncertainty small.The error for the calculated heat flow is estimated by the Gaussian error propagation for a water depth of 2200 m and becomes~22 mW/m2 (~44%). This estimation accounts for the following(i)
the errors of the pressure-depth transformation (~3%) which does not include lithostatic pressure and density variations as they are relatively small at BSR depths and near-hydrostatic conditions are assumed [61, 62]; Ganguly et al. [24] on the contrary assume lithostatic pressures within the upper sediments of the Cascadia margin that causes 8%–12% higher heat flow values that better fit to measured values;(ii)
the error of the geothermal gradient that includes the temperature variation at the BSR depth (~5%) that is the difference from the stability curve to measured in-situ values at the Pacific margin [20] and the error of the estimated seafloor depth of ~1% as the velocity within the water column varies only slightly between 1490 m/s and 1510 m/s; the error increases with increasing depth below the seafloor as it is a propagating error [63] but stays within a few percent in the upper sediment layer; we assume an error of ~3% for the determination of the BSR depth; the uncertainty in the estimation of the BSR depth is the main factor controlling the total uncertainty;(iii)
the error of the thermal conductivity (~5%) that is quite small because it is based on in-situ data from ODP Leg 170 [20, 33].


## 7. Conclusions

(i) The BSR is characterized by several properties in conventionally processed seismic sections like the phase reversal, the seafloor simulating characteristics, and the crosscutting of strata. Other characteristics require further processing efforts but deliver profound information, like the AVO effect. In this work we apply a fast qualitative analysis by comparing far to near offset CMP stacks that successfully reveals the BSR against sediment reflectors as can be seen in Figure6 (1a and 1b).(ii) Not all continuous reflectors with a phase reversal that may even show an AVO effect are reflecting the base of the GHSZ. Therefore the forward modeling of the BSR is a useful constraint for further interpretation on far offset CMP stacks. In this way both independent methods supplement each other: Figure6 (2a and 2b) reveals two phase-reversed reflectors beneath an eroded surface close to each other with strong amplitudes in the far offset stacks. We identify the BSR in the vicinity of these ambiguous reflections, calculate the geothermal gradient with the temperature at the BSR by the gas hydrate stability curve of Tishchenko et al. [22] and the temperature and depth at the seafloor, and interpolate it beneath the eroded surface. From the interpolated geothermal gradients we retrieve the potential BSR depth (Figure 9). The potential BSR depth goes along with the deeper of these two reflectors which is therefore identified as the BSR.(iii) The modeled BSR depth in this work has an uncertainty of about 35 m if it is calibrated by a regular spread of clear BSRs. For defining the depth of the gas hydrate stability zone the error of the traveltime-depth transformation of the seismic section has to be included.(iv) The calculated heat flow offers another possibility to verify or revise BSR interpretations. Reflectors that exceed its uncertainty boundaries are not confirmed as BSRs, though they may show some BSR characteristics. Instead, they may reflect fluid or gas accumulations along sediment discontinuities (e.g., turquoise reflector on the left of Figure9)(v) The heat flow trend at the Costa Rica continental margin that we derive from BSR depths has an inaccuracy of 44% and is mainly controlled by the uncertainty to estimate the BSR depth from seismic velocities. The heat flow trend fits well to the heat flow trend derived from in-situ data south of the Paleo Plate Boundary and an analytical solution [51]. Observed deviations are interpreted to be related to thermal anomalies, for example, in the frontal prism and the change of the subduction angle around 30 km landward of the trench.The broader implication of our work is that the depth of the base of the gas hydrate stability zone is determined with higher confidence, which is, for example, important for non-seismic methods to calculate gas hydrate concentrations. Additionally the derived heat flow trend and its variations from an expected trend support geological interpretation of the tectonic history of the margin.



---

*Source: 812713-2009-12-28.xml*