Heat Flow Measurements over Bottom Simulating Reflectors, Offshore Southwestern Taiwan

1 Institute of Oceanography, National Taiwan University, Taipei, Taiwan, ROC 2 Macronix International Co., Ltd., Science Park, Hsinchu, Taiwan, ROC * Corresponding author address: Prof. Chuen-Tien Shyu, Institute of Oceanography, National Taiwan University, Taipei, Taiwan, ROC; E-mail: ctshyu@oc.ntu.edu.tw The submarine accretionary wedge offshore of southwestern Taiwan consists of an intensively deformed upper slope and a lower slope characterized by mostly west-vergent anticlines and thrusts. The fast accumulation of sediment covers and seals organic material before it is oxidized within the accretionary wedge. The bottom simulating reflectors (BSRs) associated with the accretionary prism and imbricate wedge are often observed in marine seismic reflection data from continental slopes and rises. A total of 23 heat flow measurements were carried out along seismic profiles with obvious BSRs offshore of southwestern Taiwan. On average, through sedimentation-correction heat flows have increased by 10%, and the base of gas hydrate stability zone (BGHS) predicted by geotherm has up-shifted by 11%. The average depth of BGHSs is 304 mbsf versus 358 mbsf for BSRs. The average heat flow from measurements is 64 mW m, while that estimated from BSRs is 54 mW m (16 percent less) which is reasonable compared with past results of other researchers. They are both increasing from coastal waters southward to the front of the accretionary wedge. One interesting fact found is that one of our heat flow measurements showed no indication of impact friction raised temperature on data recordings. The friction may have been absorbed by gas hydrate for dissociation. (


INTRODUCTION
Gas hydrate is an ice-like compound of water with inclusions of methane, ethane, propane, hydrogen sulfide, carbon dioxide, and salt (NaCl) (McIver 1974;Yefremova and Gritchina 1981;Brooks et al. 1984;Ginsburg et al. 1992;Kvenvolden 1995;Kastner et al. 1998). High pressures and low temperatures are required for the formation of a gas hydrate. These requirements restrict its natural occurrence to deep ocean bottom sediments and areas of thick permafrost. The stability of gas hydrate is dependent on pressure, temperature, and the solubility of gas as a function of pressure and temperature in the system (e.g., Handa 1990;Zatsepina and Buffett 1997). The gas hydrate stability is significantly more susceptible to changes in temperature than pressure. There are two methods to estimate the boundary of methane hydrate stability for all combinations of temperature, pressure and salinity. The first method is based on experimental recordings of the critical temperature and pressure at the onset of gas hydrate (e.g., Handa 1990;Quinby-Hunt 1994, 1997;Brown et al. 1996). The second method is based on minimizing the Gibbs Free Energy of the system. Computer programs to determine the gas hydrate stability boundary at varied temperatures and pressures may be available (Sloan 1990(Sloan , 1998Buffett 1997, 1998;Buffett and Zatsepina 1999). Gas hydrate forms wherever appropriate physical conditions, moderately low temperature, and moderately high pressure exist. These conditions are usually found on sea floor sediments at depths greater than about 500 m (depending on the water temperature). Minimum temperature normally occurs at the sea floor; however, downward through the sediments, the temperature rises along the geothermal gradient (geotherm). At the base of the gas hydrate stability zone (BGHS), or at a depth where the gas hydrate begins to dissociate due to temperature and pressure conditions no longer meet the requirement for methane hydrate to be stable (e.g., too higher temperature). Hydrate is confined to the uppermost layer (above BGHS) of sediments, while free gas is stable below the BGHS. When the BGHS between gas hydrate above, and nonhydrated (and possibly gas-charged) sediments below, exists within the sedimentary column, the decrease in velocity and density makes for an acoustic impedance (density multiplied by sound velocity) contrast (reflectivity) at the BGHS to be prominent and consequently a negative-polarity reflector. This contrast may produce rather striking bottomsimulating reflectors (BSRs) on seismic reflection profiles. In general, BSRs are more or less parallel to the sea floor and exhibit increasing depth with water depth. This is predicted from the pressure-temperature conditions for stable gas hydrate, where increasing hydrostatic pressure and decreasing water-bottom temperature will make hydrates stable to greater depths in the sediments. The depth and in situ temperature at BGHS can be estimated when seismic velocity, pressure, and gas composition are known. Shipley et al. (1979) estimated geotherms from the depth of the BGHS on continental slopes and rises. Townend (1997) concluded that BSR data is very useful in determining offshore heat flow around New Zealand, and suggested that heat flow estimates need to be corrected for thermal effects of ongoing sediment deposition. Ganguly et al. (2000) analyzed the heat flow variations from data of BSR on the continental slope of the northern Cascadia margin, and accentuated that carrying out topographic corrections in regions of significant relief such as the continental slope are important to heat flow values. Though the heat probe measured heat flows and BSR heat flow estimates are often combined to interpret the regional thermal patterns and trends (Hyndman et al. 1993), heat flow data measured by heat probe is not always consistent with that estimated from BSRs (Davis et al. 1990). In the present study, we examined the detailed variation of heat flows measured with heat probes and compared them with those estimated from BSRs on the seismic profiles for offshore southern Taiwan. The importance of this study is: (1) that the heat flow values, variations, and trends of the continental slopes and rises of offshore southwestern Taiwan may greatly affect the gas hydrate distribution patterns in this region; and (2) that the derived BGHSs either from geotherms of heat probe data or from BSR estimates are very useful for the assessment of gas hydrate deposits, which could be commercially developed to supply huge volumes of energy to a needy world far into the future.

TECTONIC SETTING
The studied areas off southwestern Taiwan in the northern end of the South China Sea are at the junction between the Chinese continental margin and the accretionary wedge of southern Taiwan. The wedge was formed by the eastward subduction of the South China Sea lithosphere under the Luzon arc during the Cenozoic time (Big 1972;Bowin et al. 1978;Teng 1990;Huang et al. 1997). Normal faults which form horsts and grabens are observed in the shelf and slope regions of the passive Chinese continental margin located in the northwestern corner of the junction area. The frontal portion of the imbricated thrusts and related folds, observed in the lower slope region of the accretionary wedge which extended north-north-westerly toward the Chinese continental margin, were gradually buried and terminated as they reached the Chinese continental slope (Liu et al. 1997). The Chinese continental margin in this area is characterized by a host-and-graben system related to the extensional tectonic of a passive margin (Sun 1982(Sun , 1985Hu 1988). Seismic profiles indicate that large amounts of incoming orogenic sediments from the Taiwan mountain belt have caused outward growth of the accretionary wedge and westward migration of the deformation front (Liu et al. 1997;Yu and Huang 2006). The submarine accretionary wedge of offshore southwestern Taiwan consists of an intensively deformed upper slope, and a lower slope characterized by mostly west-vergent anticlines and thrusts. Many mud volcanoes and diapirs were reported both on-land and submarine southwestern Taiwan (Chow et al. 2000;Yang et al. 2004;Chiu et al. 2006). The trend evolution of the frontal folds and thrusts in the lower slope was described in detail by Liu et al. (2005).
The fast accumulation of sediment covers and seals the organic material before it is oxidized in the accretionary wedge. BSRs are often observed in marine seismic reflection data from continental slopes and rises, many of which are associated with the accretionary prism and imbricate wedge (Katz 1982;Lewis and Pettinga 1993;Ashi et al. 2002). Intensive BSRs have been observed in off southwestern Taiwan (Chi et al. 1998;Liu et al. 2006). Some abnormally geological and geochemical features (Chao and You 2006;Chuang et al. 2006;Horng and Chen 2006;Jiang et al. 2006;Lin et al. 2006;Yang et al. 2006) have also been reported. It infers that gas hydrates are widely distributed in the region.

Determination of the Temperature and Thermal Conductivity of Sediments
The heat flow data presented in this paper was collected aboard the R/V Ocean Researcher I from 2000 to 2003. Most of the heat flow sites were chosen along pre-existing seismic profiles of obvious BSRs (Figs. 1, 5). The Lister-type high-resolution heat probe used for heat flow measurements was developed at the Institute of Oceanography, National Taiwan University. It allows temperature resolution to be better than 0.1 mK in the range of -1 to 25°C. Although the heat probe has the capability of multi-penetration, due to the time constraint, only one penetration was attempted for all measuring sites. During each penetration, the instrument remains undisturbed in the seafloor for about 30 minutes to measure temperature and thermal conductivity in situ (12 minutes for the geotherm; 8 minutes for the thermal conductivity). The equilibrium ambient temperature of sediment was obtained by fitting the temperature decay curve with the thermal decay model of a cylinder (Bullard 1954;Carslaw and Jaeger 1959), and regressing with respect to parameters of ambient temperature, initial temperature rise, thermal conductivity, and heat capacity of the sediment (Shyu and Chang 2005). The same regressing process was applied to the heat pulse decay and extrapolated to infinite time (asymptotic solution) for solving the thermal conductivity of sediment (Hyndman et al. 1979). The data at shallow penetration sites generally reveals abnormal geotherms compared to other sites. This may be related to bottom water temperature variations in this area; consequently, these data were excluded from analysis. The resulting temperatures and thermal conductivities at 23 sites are displayed in Fig. 2 and listed in Table 1.

Estimation of the BGHS from Geotherm
The pressure at the depth of the BGHS is the sum of the pressure of the seawater column overlying the seafloor and the pressure of the sediment column overlying the BGHS. Although hydrostatic pressure (from the pore-water of the sediment) of the sediment column could be used for the calculation, lithostatic pressure (from the gross density of the sediment), on the other hand, may give heat flow values which are closer to those determined by seafloor heat probe measurement (Davis et al. 1990;Ganguly et al. 2000). Nevertheless, a suitable criterion for the model selection (hydrostatic or lithostatic model) has not been proposed. In the study area, the pressure is close to lithostatic due to fluid expulsion along faults or the mud diapir pierced the sea floor Chiu et al. (2006). In the following, we convert the BGHS to pressure by taking both hydrostatic pressure for the water column and lithostatic pressure for the sediment column.
Sediments on the lower slope of the accretionary wedge offshore southwestern Taiwan consist of silts and clays. In our pressure calculation for the methane hydrate stability curve, the density for the water and the sediments in the study area were assumed to be 1.025 and 1.7 g cm -3 (Chen 1997), respectively. The depth (in mbsf) at BGHS for each heat flow site is determined at the intersection between the geotherm extrapolated from measured temperatures and the gas hydrate boundary curve. The pressure-temperature empirical equation for methane   hydrate stability (gas hydrate boundary curve) in seawater of 35 0 00 salinity is given in the following (Dholabhai et al. 1991;Dickens and Quinby-Hunt 1997): where T is the temperature in K and P is the pressure in MPa. A typical example in Fig. 3 shows a simple geometrical procedure of obtaining the depth of BGHS. The resulted depths of the BGHS for all heat flow measurement sites are summarized in Table 2, most of which are around 250 ~ 400 mbsf.

RESULTS AND DISCUSSION
The discrepancy between the theoretical BGHS (predicted by the extrapolation of the geotherm) and the observed base of gas hydrate bearing sediments may be related to miscellaneousneous mechanisms such as the effect of the interstitial water chemistry (Lu and Matsumoto 2001), the capillary effect (Handa and Stupin 1992), the gas composition (Waseda and Uchida 2004), the uplift of the BGHS (Foucher et al. 2002), and the change of thermal gradient with depth. We will discuss the heat flow distribution, trend, and the comparison of the BGHSs predicted by geotherms with the observed BSRs on seismic profiles below.

Geotherms and Thermal Conductivities
In general, the geotherms show less departure from a straight line than those of thermal conductivities. Average geotherms are determined from the best linear fit of equilibrium tem- peratures for each thermistor versus depth. To assure a better fit of the regression line, we did not take those data points that deviated far away from the regression line into account. The extremely low or slightly negative geotherm that occurred at site No. 15 ( Fig. 2; Table 1) is of the most interest to us because there is no indication of friction raised temperature on temperature recording when the heat probe impacted the ocean floor (Fig. 4). This may happen if the Fig. 4. Temperature recordings at site No. 15. Note: (1) There was no recorded abnormal raised temperature from friction when the heat probe was inserted in the seafloor; and (2) irregular oscillation temperature was recorded by the sensor #7 on the seafloor surface.
ocean floor material is very soft and very close to a fluid state (e.g., gas in seawater). According to the present data collected on location (Fig. 5) we conjecture that for the reduction of friction, the material is likely to be an outcrop of powdery gas hydrate disseminated over the sediments. Thus a small amount of friction generated heat by the probe penetration is taken into the gas hydrate dissociation process (endothermic process). That is, the powdery ice-like solid gas hydrate absorbed the friction heat and changed to methane at the same temperature. Additional auxiliary evidence is the temperature fluctuation of bottom water which could be affected by the upward dissociation methane flux (Fig. 4 ). Most heat flow measurement sites are located on the muddy floor. Figure 6 (upper panel) shows how the heat flow increases to the southwest with water depth; this trend could be explained by the lower geotherm close to the coast at the top of the sea floor having been reduced by a high rate of sedimentary deposition while the higher geotherm, in the front of the Parenthetical values at sites are the sedimentation-corrected BGHS (BGHSc) in mbsf. BGHSc at site 714ht2 is not predicted due to extremely low geotherm ( 0). accretionary wedge, came about as a result of the northeastward subduction of the China Sea lithosphere. Sediments are derived mainly from Shou-Shan Gully and Kaoping Canyon passages. Abundant terrigenous clastic materials continuously pour into the shelf and slope by river propagation and longshore erosion. In general, sedimental skewness is greater than 0.1 which implies that coarse sediments are dominant in this area (Chen 1997). This may explain partly why the deviation of thermal conductivity with respect to a line is obviously greater than that of temperature (Fig. 2). The conductivities are generally uniform with a mean value of 1.00 ± 0.10 Wm -1 K -1 which is higher than that of typical marine sediments (~0.9 Wm -1 K -1 ). Experiments show that pure gas hydrate has a low thermal conductivity of about 0.5 Wm -1 K -1 . In contrast, water ice, a substance to which hydrate is often compared, has thermal conductivity of 2.23 Wm -1 K -1 (Sloan 1998). Therefore, our result of slightly higher thermal conductivity may imply that the sediment contains hydrate and ice mixtures. However, it is impossible to determine whether the sediments contain gas hydrate, or what the percentage of the gas hydrate is contained by the level of thermal conductivity. Thermal conductivity generally increases to the southeast. Although there are some differences in mean conductivity between sites, variation with respect to depth is not obvious at most of the sites. Because of the lower conductivity values of 0.87 and 0.84 Wm -1 K -1 , and temperature gradients of 0.039 and 0.042 K m -1 at sites No. 3 and No. 4, respectively, heat flows at these sites are lower (Table 1).

Effects of Sedimentation
Since heat flow measurements are made by a Lister-type heat probe, which can only penetrate no deeper than 4.5 mbsf, the temperature gradient across the thin upper layer of sediments may be depressed due to sedimentation. Each grain of sediment arrives at the sea floor at the temperature of the sea bottom water through which it has fallen, which is normally lower than the temperature of those sediments already on the seabed. As the water-sediment boundary is slowly rising, the temperature at every point below it must adjust upward towards a new equilibrium value. This means that part of the conducted heat must be used to warm the material through which it is passing, and the heat flow value at the sea floor is reduced.
The fine particle turbidity sediments of offshore southwestern Taiwan are mainly from erosion of the southern end of the western foothills and western flank of the Central Range. Sedimentation rates estimated from seismic profile analysis range from 0.1 to 0.5 cm yr -1 (Liu et al. 1997). With this extremely high sedimentation rate, the temperature and geotherm disturbance becomes significant in comparison with measurement errors. A better way to exclude the effects of sedimentation is to estimate the thickness of recently deposited sedimentary layers from seismic reflection profiles and the sedimentation rates at each site. However, sedimentation rate may vary greatly depending on site locations. In the following calculation, according to Liu et al. (1997), we have adopted an average sedimentation rate of 0.3 cm yr -1 and estimated the thickness of the recently deposited sedimentary layers at each site. The correction equation is given by Jaeger (1965;adapted by Kappelmeyer and Haenel 1974): where GRAD is the corrected thermal gradient, GRADS the thermal gradient with sedimentation, and Ω is given by: where v s is the mean uncompacted sedimentation rate, τ is the duration of sedimentation, k is the thermal diffusivity , the ratio of thermal conductivity (K) to heat capacity ( ρc ) of the sediment, and "erf" denotes the error function. In our calculation, we have assum ρc = (5.79 -3.67K + 1.016K 2 )10 6 Jm -3 K -1 as that proposed by Hyndman et al. (1979). The uncorrected and sedimentation-corrected geotherms, BGHSs, and heat flows are given in Table 2, Figs. 5 and 7. The geotherms and heat flows are both raised to about 10% on average. Consequently, the intersection of geotherms and gas hydrate stability boundary curves become relatively close to the seafloor surface. That is, the sedimentation-correction leads to an up-shift of the BGHS by about 11% of the uncorrected BGHS. As we expected, the sedimentation-corrected BGHS (BGHSc) is decreasing with heat flow (Fig. 8).
The depths of BSRs in Table 2 were converted by the following empirical depth-sound travel time relation given by Hamilton's geoacoustic model (1980), which was based on data from 20 areas of terrigenous sediments in the North Pacific and adjacent areas (applicable to silt, turbidites, mudstone shale): where t is the seismic P-wave one-way time in sec and Z is the depth below seafloor in mbsf. The theoretical BGHS on the Blake Ridge in the western Atlantic is deeper than the actual observed gas hydrate based on ODP Leg 164 drilling well data (Paull and Matsumoto 2000). Henry et al. (1999) attributed the downward shift of the BGHS (about 9 ~ 25% deeper) to capillary effects. The coring data taken from survey wells in the eastern Nankai Trough, on the other hand, indicate that the predicted depth of the BGHS is shallower than the base of the gas bearing sandstones (less than about 15%) and the depth of the BSRs (less than about 14%) (Matsumoto et al. 2004). They proposed more complicated possible mechanisms to explain the discrepancy. The present data (Table 2; Fig. 9) has shown that most of the BGHScs (17 of 23) are shallower than the depth of BSRs. The average depth of BGHSc is 304 mbsf which is about 16% shallower than that of BSR (358 mbsf). Beside the above mechanisms which may lead to error, the differences may be relevant to the uncertainties of the sedimentation rate used in Eq. 2 and the BSR depth conversion Eq. 4 (due to uncertainty of seismic velocity). We would rather have an average depth close to 304 mbsf than 358 mbsf for BSR. This will be discussed in the following section.
Contrary to the common belief that the BGHS is theoretically increasing with water depth, we found that this relation is not obvious in our examples (Fig. 10). This may be related to the fact that the seabed temperature is in a non-equilibrium state due to continuous crustal movement.      Yamano et al. (1982) first applied the base depth of BSRs to estimate heat flow in the Nankai accretionary prism, and found that heat flows estimated by the methane hydrate stability boundary diagram were almost consistent with the values measured by heat probe. The same method has been applied by others (e.g., Hyndman et al. 1993;Townend 1997;Ashi et al. 2002). The difficulties are the estimation of parameters such as the density of hydrate bearing sediments for pressure calculations and thermal conductivity for heat flow calculations. For a more detailed estimation of these parameters, refer to the article by Ashi et al. (2002). The method used here is basically similar to that initially proposed by Yamano et al. (1982) and discussed in detail by Minshull and White (1989). In the estimation of heat flows from BSRs, all previous used parameters like the temperature at the seafloor surface, thermal conductivity, density of gas hydrate bearing sediments, and the methane hydrate stability equation were adopted except the depths of BSRs being derived from Eq. 4. Accordingly, estimated heat flows increase southward and reach a maximum in a zone which is slightly eastward of that from direct measurements (Figs. 6, 11). The average heat flow estimated from the BSRs is 54 mW m -2 , which is slightly below that from measurements (58 or 64 mW m -2 for sedimentationcorrected heat flow, Table 2; Fig.12). They both vary positively with water depth (Figs. 6,11,12,and 13). Since the dominant structural patterns in the study are fold-and-thrust structures or products of continuous collision, higher heat flows are supposed to occur reasonably. For example, if we expect to have a consistent heat flow value from heat probe measurement and from BSRs, then according to Table 2 we need to have a greater geotherm derived from BSRs or a thinner thickness of BSRs for geotherms. BGHSc (bottom) on plots of water depth versus BSR depth below seafloor. Four curves with different heat flow values represent theoretical relationships between BSR depths and water depths (assuming average temperature on the seafloor = 2.850°C and average thermal conductivity of the sediments = 1.00 W mK -1 as shown in Table 1). Fig. 13. Trends of measured heat flows (top) and estimated heat flows from BSRs (bottom) with water depth.

CONCLUTIONS
Twenty-three heat flow sites along pre-existing seismic profiles with obvious BSRs of offshore southwestern Taiwan were surveyed by a Lister-type high-resolution probe. Heat flows at these sites were also estimated from multi-channel seismic reflection data. The depths of BGHSs have been predicted by the intersections between the geotherms and the gas hydrate boundary curves, which have been compared with that of BSRs. The conclusions are summarized as follows: 1. Through sedimentation-corrected geotherms and heat flows, both are raised to about 10% on average. The correction leads to an up-shift of the BGHS by about 11% on average. 2. The average heat flow estimated from the base depth of BSRs is 54 mW m -2 , which is below the measured heat flow of 64 mW m -2 with sedimentation-correction, on average. We would rather have a greater geotherm derived from a BSR or a thinner thickness of BSR for consistency considerations because 54 mW m -2 is too low on a continuous collision zone with fold-and-thrust structures. Nevertheless, they are both increasing approximately southward. 3. When comparing the depths of sedimentation-corrected BGHSs and that of BSRs, we may roughly estimate the range of BGHSs to be between about 300 to 360 mbsf. The average depth of BGHSs is 304 mbsf, which is about 16% shallower than that of BSRs in the study area. This amount of difference seems reasonable compared to previous studies. 4. Contrary to the common belief that the BGHS is theoretically increasing with water depth in our examples, we found that their relation is not obvious. This may be related to the fact that the seabed temperature is in a non-equilibrium state due to continuous crustal interaction and movement. 5. One heat flow measurement showed that there was no indication of friction heat raised temperature when the heat probe impacted on the seafloor. We surmise the ice-like solid gas hydrate has absorbed the small amount of friction heat for dissociation and is changed to a methane state. However, the statement is based on data from one measurement only, further experimentation is highly desirable.