Distribution of gas hydrates on continental margins by means of a mathematical envelope: A method applied to the interpretation of 3D Seismic Data

A 3D seismic volume from the Nankai Trough accretionary wedge (SE Japan) is used to evaluate the subsurface distribution of gas hydrates as a function of structural and stratigraphic complexity, variable heat flow patterns and the presence of subsurface fluid conduits. Eleven equations were modified for depth, pressure, and temperature, modeled in 3D, and compared with the distribution of Bottom‐Simulating Reflections (BSRs) offshore Nankai. The results show that the equations produce overlapping—and thus potentially consistent—predictions for the distribution of BSRs, leading us to propose the concept of a “BSR Stability Envelope” as a method to quantify the subsurface distribution of gas hydrates on continental margins. In addition, we show that the ratio (R) between shallow and deep BSRs of seven subenvelopes, which are defined by BSR stability equations, indicates local gas hydrate equilibrium conditions. Values of R < 1 relate to cooler regions, whereas when R > 1 the majority of BSRs are located in warmer structural traps. The method in this paper can be used to recognize any divergence between observed and theoretical depths of occurrence of BSRs on 3D or 4D (time lapse) seismic volumes. In the Nankai Trough, our results point out for equilibrium conditions in BSRs located away from the Megasplay Fault Zone and major thrust faults. This latter observation demonstrates the applicability of the method to: (a) the recognition of subsurface fluid conduits and (b) the prediction of maximum and minimum depths of occurrence of gas hydrates on continental margins, under distinct thermal and hydrologic conditions.

Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Introduction
[2] The subsurface distribution of gas hydrates in offshore basins is important not only in economic terms, but also for geohazard and climate change studies [Harvey and Huang, 1995;Clennell et al., 1999;Bouriak et al., 2000;Davie and Buffet, 2001;Sultan et al., 2004a]. Information on the thermal conditions of continental margins can also be derived from understanding the conditions in which gas hydrates are formed in the subsurface [Trehu et al., 2004]. As they represent a transient surface between solid-liquid and free gas, gas hydrates are known to form within specific subsurface pressure and temperature conditions [e.g., Kinoshita et al., 2011], as shown offshore Nankai (SE Japan) at IODP Sites C0001 and C0008 [Kimura et al., 2008;Harris et al., 2011] and after direct seafloor heat flow measurements [Yamano et al., 2003]. On seismic data, most gas hydrate accumulations are recognized as negative-polarity reflections that are subparallel to the seafloor. These reflections are named Bottom-Simulating Reflections (BSRs), and correspond to the base of a gas hydrate stability surface (BGHS) where solid-liquid gas clathrates change into free gas [Shipley et al., 1979;White, 1979;Singh et al., 1993]. The presence of reflections of negative polarity indicates that an interval of lower velocity directly underlies higher velocity sediment above. For instance, vertical reflection profiles acquired at ODP Site 997 indicate the velocity of sediments overlying a BSR to approach 1850 m/s, whereas a velocity of 1400 m/s is recorded beneath the BSR due to the presence of free gas [Holbrook et al., 1996].
[3] Published data reveal that complex parameters control BSR generation, some of which are responsible for systematic errors when modeling the sub-surface distribution of gas hydrates [e.g., Rempel and Buffet, 1997;Xu and Ruppel, 1999;Bouriak et al., 2000;Wood and Ruppel, 2000;Grevemeyer and Villanger, 2001;Chand and Minshull, 2003]. In the Nankai Trough the depth of BSRs has been considered as responding primarily to regional heat-flow trends, but an uncertainty of $20% is observed due to recent horizontal shortening . This latter parameter is seldom considered in mathematical models.
In parallel, discrepancies in the estimated volume of gas hydrates in the Gulf of Mexico and China Sea were found to result from variations in water depth, temperature and subsurface heat flow [Milkov and Sassen, 2001;Trung, 2012]. These discrepancies were addressed by producing Gas Hydrate Stability Zone (GHSZ) thickness maps, based on previous work recognizing the existence of mathematical limits to gas hydrate subsurface distribution, but not modeling it in 3D [e.g., Sloan, 1998;Helgerud et al., 1999].
[4] Notwithstanding the latter efforts, no systematic methods have yet been proposed to quantify subsurface distribution of BSRs on 3D seismic data. Particularly important is to mathematically transform GHSZ thickness maps in BSR distribution envelopes capable of representing the upper and lower depth limits in which BSRs are expected to form. When precisely defined, these mathematical envelopes can be tied to elasticwave data used to quantify the relative percentages of gas hydrates in the subsurface, thus helping to understand gas hydrate distribution in structurally complex areas [Berge et al., 1999;Helgerud et al., 1999;Pecher et al., 2005;Hornbach et al., 2008;Kinoshita et al., 2011].
[5] Acknowledging the importance of defining such mathematical envelopes, this paper models the theoretical depths of BSRs in the Megasplay Fault Zone (MFZ) and upper Imbricate Thrust, Zone (ITZ) of the Nankai Trough, both comprising regions of active thrust faulting [Martin et al., 2004;Kinoshita et al., 2011] (Figures 1a-1c). The MFZ and ITZ have been recognized as regions of distributed deformation, generated during the Quaternary, in which BSR distribution clearly depends on local variations in heat flow, tectonic uplift, and sedimentation rates [Ashi et al., 2009;Kinoshita et al., 2011]. We demonstrate that BSRs offshore Nankai are distributed within a BSR  stability envelope that is quantifiable on 3D seismic data [e.g., Dubrule et al., 1998]. The ratio (R) between shallow and deep BSRs within the stability envelope is proposed in this work as a method to identify local variations in temperature, pressure, lithology, pore pressure, salinity, tectonic uplift, and subsidence.
[6] The paper starts with a description of the parameters governing the subsurface distribution of gas hydrates. This is followed by: (a) a description of data and methods utilized, (b) an overview of the geology of the Nankai Trough, and (c) a summary of the lithostratigraphic units drilled by the Integrated Ocean Drilling Program (IODP) in the area of interest to this study ( Figure 1a). We show gas hydrate distribution to follow a BSR stability envelope influenced by geological and physical processes. The methodology behind the definition of this stability envelope is then described in detail. In essence, our method comprises the compilation of eleven equations for Base of Gas Hydrate Stability (BGHS) surfaces, which are used to analyze the distribution of BSRs offshore Nankai. At the end of the paper, we assess the effect of local variations in P, T, and salinity on the BGHS surfaces and discuss the concept of a BSR stability envelope on continental margins. In summary, this paper aims to: [7] 1. Test the validity of known empirical solutions for the stability of gas hydrates in the subsurface.
[8] 2. Investigate to what degree stability equations vary with respect to specific geological conditions, such as lithological types, structure, and thermal regimes.
[9] 3. Understand the potential role of subsurface heat flow, tectonic structure, and surface processes when characterizing BSR distributions.
[10] 4. Revise the use of singular boundaries for gas hydrate stability on continental margins and to propose a region/zone concept for the estimation of BSR stability envelopes.
[11] A major advantage of the study area is that most BSRs are discordant in relation to subsurface strata ( Figure 2). BSRs and gas-related anomalies occur on the crest or flanks of thrust anticlines, forming clear flat spots and anomaly clusters of different geometry to involving strata ( Figure 2). In addition, folds, reactivated rollover anticlines and thrusts are truncated by unconformities, and correlative stratigraphic surfaces, cored by IODP expeditions 314, 316, and 333 (Figures 1d and 2).
Hence, the method in this paper can be used to investigate the control of parameters such as salinity, temperature, and hydrostatic pressure on the subsurface distribution of BSRs and associated gas hydrates. It will assist the interpreters in finding the subsurface envelope in which gas hydrates are expected to form if subsurface T, P, and petrophysical conditions are adequate.

Parameters Governing the Subsurface Distribution of Gas Hydrates
[12] The subsurface distribution of gas hydrates is governed by a set of parameters matching the conditions required to achieve a three-phase equilibrium state where solid clathrate, free gas, and water coexist [e.g., Grevemeyer and Villinger, 2001;Kinoshita et al., 2011]. Water depth, temperature at the seafloor, and geothermal gradients are the primary controls of this lower limit of stability [Wang et al., 2006]. Mixtures of methane and other gas components, pore-water salinity, in situ processes of methane activity, pore size, and pore distribution are just a few other parameters to consider when formulating an approach to stability (see Appendix I in supporting information 1 ).
[13] Dickens and Quinby-Hunt [1994] showed that variations in BSRs subsurface distribution for similar P-T conditions can also derive from the diffusion of dissolved ions, the incorporation of trace gases in BSRs, and the impact of surrounding sediments. Dickens and Quinby-Hunt [1994] solution is helpful to determine the depth of a gas hydrate stability zone governed by hydrostatic pressure at any subsurface point. However, this same equation is only applicable to a small range in pressures (2.5-10 MPa) [Sloan, 1998] and to pure waterpure methane systems, which are seldom recorded in the subsurface. Subsequently, Peltzer and Brewer [2000] concluded that a second-order polynomial equation fitted to the Dickens and Quinby-Hunt [1994] stability data does a better job analyzing gas hydrate formation at relatively high-temperature ranges.
[14] More recent work conducted in the Gulf of Mexico by Milkov and Sassen [2001] presented gas-hydrate thickness maps derived from logarithmic functions in Sloan [1998]. In contrast to Dickens and Quinby-Hunt [1994] solution, Sloan [1998] work predicts stability pressures for known temperature ranges. Best-fit stability curves for pure methane and two methane and other gas mixtures (structures I, II, and H), together with a solution for sediment temperature at a depth below sea floor, were also compiled by Milkov and Sassen [2001]. This latter approach followed that of Sloan [1998] in which thickness maps were calculated from the intersection of the stability curves with sediment temperature data using the Newton-Raphson method of solving the equations simultaneously for depth. Trung [2012] followed the same method as Milkov and Sassen [2001] to calculate a BGHS in the South China Sea at a pore water salinity value of 3.5%. An alternative geothermal calculation was also used based on well data from the South China Sea.
[15] These studies highlight the existence of two groups of mathematical and empirical models used to characterize the distribution of gas hydrates in marine sediments: (a) models considering equations solved for the conservation of momentum, fluid mass, and energy for transient and steady-state systems [Xu and Ruppel, 1999;Rempel and Buffet, 1998;Egeberg and Dickens, 1999;Buffet, 2001] and (b) simpler models taking into account the energy conservation of an essentially transient regime [Chaouch and Briaud, 1997;Delisle et al., 1998;Sultan et al., 2004aSultan et al., , 2004bLu and Sultan, 2008;Sun et al., 2011].
[16] The first group of models uses data on the conservation of methane and supply of gas into hydrate systems. For instance, Xu and Ruppel [1999] and Rempel and Buffet [1998] studied the rate of hydrate formation through advectivediffusive flow coupled with heat transfer whenever methane is present in the zone of stability. In addition, Egeberg and Dickens [1999] took into account the effects of chloride in pore fluid when determining hydrate formation as a function of depth. Davie and Buffet [2001] combined both chlorinity and methane content within pore fluid to show methane solubility profiles can influence over formation and distribution of gas hydrates. They have also shown that chlorinity values are influenced by increasing fluid flow rates.
[17] The second group of equations, assuming the conservation of energy in methane systems, emphasizes the control of latent heat on the formation and dissociation of gas hydrates [Chaouch and Briaud, 1997;Delisle et al., 1998]. Early models, however, were regarded as inadequate by neglecting the effects of gas composition and concentration. Sultan et al. [2004a] first addressed these inadequacies by developing a model based on the enthalpy law that accounted for pore water chemistry and sediment pore-size distribution. In Sultan et al. [2004a], the formation and dissociation of gas hydrate responds chiefly to changes in pressure, temperature, and gas concentrations. Based on this latter work, Lu and Sultan [2008] established empirical expressions for the formation, fraction, and density of gas hydrate for a range of systems that considers gas composition, salinity, and pore size as main controlling parameters. More recently, Sun et al.[2011]producedanempirical formula for the stability of gas hydrates in the South China Sea from seafloor sediment and water samples analyzed in multistep experiments. Sun et al. [2011] results showed that stability is mainly affected by pore water salinity due to a comparable shift in dissociation temperatures in seafloor water and sediment tests relative to a pure water system.
[18] The interpretation that different parameters have different degrees of control over the dissociation of gas hydrates in the subsurface is, however, the subject of much debate. Mixtures of methane and other gas components, pore-water salinity, in situ processes of methane activity and the pore size and distribution are just a few of the considerations when formulating an approach to stability. The local geothermal gradients and heat flows are also important factors. However, in many studies local variations in these values due to tectonic structure and advecting fluids are either considered linear with depth or not known at all [e.g., Pecher et al., 2005;Hornbach et al., 2008].

Seismic Data Interpretation and BSR Envelope Modeling
[19] The interpreted 3D seismic volume was acquired by Petroleum Geo-Services (PGS) in 2006 using two arrays of 28 Sodera air guns totalling 51 l in volume [Moore et al., 2007[Moore et al., , 2009. Four receiver streamers were used, each with a length 4500 m and a spacing of 150 m. Individual hydrophone cables carried 360 receiver groups at 12.5 m spacing. With this array setup of 4-streamer/2-source were produced eight source-receiver common midpoint (CMP) profiles per sail line, each with a spacing of 37.5 m and nominal 30-fold data. In total, the interpreted 3D volume covers a total area of 585 km 2 , with the longer axis of the survey orientated at 330.1 (Figure 1b). Data processing included 3D prestack depth migration (PSDM), which provides more detailed imaging of faults and small-scale structures [Moore et al., 2009].
[20] In order to map the observed BSRs, the 3D seismic volume was interpreted every 5 inlines, or every inline in the regions of higher structural complexity. Mapped BSRs represent the transition from concentrated solid hydrate to free gas . This change in bulk pore fill is often displayed as high amplitude reflections of reverse polarity that crosscut recognized bedding planes . Highamplitude BSRs or other acoustic anomalies are thus shown as subgeometric features of high reflection strength (Figures 2 and 3).
[21] BSR envelope modeling, and resulting surface area analyses, were undertaken in Petrel© by converting the anomalies into polygons (Figures 2  and 3). This was followed by a surface conversion with a zero expansion boundary of the polygon, allowing the interpreted BSRs to be quantified without overestimating or altering BSR surface area. The surfaces were then categorized into three structural regions on the Nankai continental slope ( Figure 1c). Graphical solutions showing the distribution of anomalies across the three zones were prepared using Grapher© 8 software.

Quantification of Gas Hydrate Distribution on Seismic Data
[22] In order to investigate the existence of a stability envelope in Nankai, eleven (11) equations were used and recalculated for T, P, and depth in this work ( Figure 4). The hydrostatic pressure components of the eleven chosen equations were calculated using the density-gravity-depth relationship in equation (1), at distinct points in the grid: [23] P h 5 Hydrostatic pressure (MPa) [24] q 5 Density of seawater 5 1024 kg m 23 [25] h 5 Height of water column above seafloor (m) [26] g 5 gravity constant 5 9.81 m s 22 [27] Details on the mathematical transformations undertaken in this work are given as supporting information (Appendix I in supporting information).
[28] The last step of our method included 3D seismic modeling on Petrel @ , in which the different BGHSs equations were plotted in 3D against the interpreted seismic data, and individually in a graphing system to show the variation in their depth and linear shape. In practice, this signifies that equations were computed in Petrel © as 880 grid-point surfaces to define eleven base gas hydrate stability surfaces (BGHS). Relative depths and BSR surface areas were then calculated for each BGHS (Figures 4 and 5). The temperature value utilized in the equations were averaged from Harris et al. [2011] for Expeditions 315 and 316, and estimated to be 4 C at the depth of occurrence of the BSRs.

Geological Evolution of the Nankai Trough
[29] The Nankai Trough, located offshore southeast Japan, is part of a subduction complex developed between the Philippine Sea Plate and the Japan island arc system (Figure 1a). In the study area, the Philippine Sea Plate underthrusts the Eurasian Plate at a rate of $4.1-6.5 cm/y [Seno et al., 1993]. Subduction of the Philippine Sea Plate generated the Nankai Trough and its frontal accretionary prism [Taira et al., 1992] (Figure 1c).
[30] To the North of the subducting plate, the Shikoku Basin was formed by back-arc spreading processes behind the Izu-Bonin arc complex (25-15 Ma) [Okino et al., 1994]. This seafloor spreading event occurred in three phases, separating the older Kyushu-Palau Ridge from the presently-active Izu-Bonin arc (Figure 1a). Intraplate volcanism in the Shikoku Basin, associated with late stage rifting until $7 Ma, formed major seamount chains [Chamot-Rooke et al., 1987;Ishii et al., 2000], which have a significant impact on the tectonic dynamics of the subduction zone [Ike et al., 2008a[Ike et al., . 2008b].

Main Stratigraphic Units
[31] The lithostratigraphy of the study area is known from several IODP Sites. Site C0001 is situated on a small plateau seaward of the Kumano Basin ( Figure 1b). Sediment cores revealed two lithological units at this site. Unit I was divided into three mud-dominated subunits (Ia-Ic) based on their petrological character [Ashi et al., 2009). Unit I has been interpreted as belonging to a Pleistocene-Holoceneslopeapron.UnitIIislatePliocenetolate Miocene in age and comprises green-gray bioturbated muds with zones of darker green sediments indicative of an higher clay content. Local wavy laminae are also present. Unit II represents a change into upper accretionary prism deposits.
[32] Site C0008 and C00018 drilled hemipelagic mud in Unit Ia, intercalated with thinly bedded sand inclusions and volcanic ash layers containing significant glass and pumice (Figure 1d). Unit Ib crosses a series of sandy turbidite layers intercalated with silts and clays. Clayey gravel with mudstones and volcanic pumice separate Units Ia and Ib (Expedition 333 Scientists, 2012) (Figure 1d). Unit II is a sand unit with minor gravel constituent.
[33] At Site C0004, Unit I comprises a calcareous nannofossil-rich mud with lesser amounts of  siliceous material and ash [Screaton et al., 2009]. Unit IIa was described as a mass transport complex (MTC) of brecciated gravels in a hemipelagic muddy background with minor sand input. Unit IIb is dominated by dark greenish gray mud with rare turbidite sands and volcaniclastics.

Structural Setting
[34] Faults are essentially found in clusters below 500 m below the seafloor (mbsf), above which smaller, rarer occurrences are observed on seismic data ( Figure 2). Faults are also clustered below 900 mbsf, in Unit IV. From this, Screaton et al.
[2009] interpreted three separate deformation phases from kinematic studies of planar and linear features: (a) thrust faulting with a potential strike-slip element, reflecting northwest-southeast shortening early in the evolution of the Nankai region, followed by (b) two Late Pliocene-Quaternary phases of normal faulting. Normal faulting is the dominant structural characteristic of the slope basin, whereas reverse and strike-slip faults govern the structure of the deeper formations (Figures 2b  and 2c).

Distribution of Acoustic Anomalies and BSRs
[35] In the study area, we defined three distinct zones based on their characteristic subsurface structure (Figure 1c). Zone 1 shows complex bathymetric features in the form of prominent ridges running along the upper continental slope. These ridges are related to thrusting of the Megasplay Fault (MSF) over accreted strata (Figures 1c and  2a). In contrast, Zone 2 extends across a perched slope basin on the upper part of the accretionary prism, onto a relatively flat seafloor section of the upper ITZ (Figures 1c and 2b). Zone 3 is defined by the shallow expression of thrust anticlines that form the seaward margin of the ITZ (Figures 1c  and 2c).

Zone 1
[36] The most distinguishable BSRs appear in Zone 1. For example, Figure 2a shows continuous BSRs exhibiting high amplitude, reverse polarity reflections cutting across dipping strata. This BSR extends for $1.9 km in a NW-SE direction, with reflections becoming discontinuous toward the splay fault zone. Although relatively flat, the shape of the BSR changes relative to the seafloor bathymetry (Figure 2a). At places, BSRs also appear vertically stacked as in Figure 2c, in which a BSR with lower amplitude is observed $ 120 m below a high-amplitude BSR.
[37] Shallower sets of high amplitude seismic anomalies, some of which likely to comprise gas hydrates, appear seaward of the seafloor ridge (Figure 2a). They are positioned at the boundary of forearc basin sediments and appear to follow this same boundary, $200 m below the seafloor. BSRs appear more scattered to the Northwest, with several reflections vertically stacked (Figure 2a). Where the seafloor shows a prominent ridge, a $100 m downward shift of BSRs is observed. Southeast of this ridge BSRs show a laterally continuous pairing of high-amplitude reflections (Figure 2a).

Zone 2
[38] Anomalies in Zone 2 occur atop and at the flanks of thrust anticlines (Figure 2b). They are also observed in the upper ITZ basin where basinfill sediments are much thicker and display tight folds (Figures 1c and 2b). BSRs on the flanks of thrust anticlines have varied distributions and sizes, following changes in the geometry of thrust folds. There are rare occurrences of flat, BSRs crosscutting background strata in Zone 2, a character likely related to thrusting and anticline uplift. In most of Zone 2, BSRs tend to favor the flanks of thrust anticlines (Figure 2b).

Zone 3
[39] Zone 3 is dominated by the presence of two anticline folds. Figure 2c shows the landward limb of an anticline as it enters the study area. Vertically stacked, high-amplitude BSRs cut across the dipping reflections and follow a relatively flat seafloor bathymetry at $350 mbsf. Some of these BSRs can be laterally traced into Zone 2 ( Figure  1d). In Figure 2c, the high-amplitude anomaly to the Northwest of the section does not appear to cross the inner part of the imaged thrust anticline (see arrow to the left of Figure 2c).

Depth Distribution
[40] Mathematically, a central cluster of BGHS is observed in our models (Figure 4). Within this central cluster are equations (2)-(4) and the different transformations of equations (5), (6), and (9). Significantly, a combination of equations (2) and (10) comprises the deepest and stability envelope boundary, while equation (11) forms its upper boundary (Figure 3). In order to analyze BSR distribution within the calculated stability envelope, we divided the subsurface strata into seven subenvelopes, A to G limited by eight principal empirical equations (Figures 5 and 6). These eight empirical equations divide the BSR stability Equation 4 9 (with Equation 8)

Equation 2
Equation 10 (n = 140) Figure 6. Graph showing the surface area distribution of BSRs within the gas-hydrate stability envelope. Note the larger surface area of BSRs in subenvelopes F and G in Zone 2. Zone 3 shows the larger surface area in the shallower subenvelopes. Zone 1 shows a relative homogeneous BSR distribution. See text for a discussion on surface area distribution of BSRs in the study area. envelope in relatively evenly spaced subenvelopes ( Figure 6).
[41] After mathematically transforming the eleven Base Gas Hydrate Stability equations, or BGHS, we represented them in 3D and analyzed BSR distribution within the interpreted seismic units and calculated BGHS in Nankai (Figures 3 and 5). The overall shapes of the BGHSs show a thickening of the stability envelope with depth, a characteristic anticipated from published data. At greater water depths, BGHSs spacing is also increased, a character denoting wider regions where BSRs can occur ( Figure 5). Minor seismic anomalies, likely not related to gas hydrate accumulations as they follow stratigraphic boundaries along the structure of thrust anticlines, occur outside the stability envelope defined ( Figure 5). The only exception is the large BSR mapped under the prominent ridge to the northwest of MFZ (Zone 1), where small portions of the BSR's landward and seaward edges intersect the computed BGHS ( Figure 5). In the remainder of Zone 1, BSRs are more evenly distributed within the modeled stability envelope.

3D Surface Area Distribution
[42] The surface area of BSRs is plotted graphically in Figure 6 for each subenvelope. BSRs in Zone 1 yield relatively small surface areas compared to the other two zones. Importantly, the relatively cold Zone 2 peaks in cumulative surface area in subenvelope F ( Figure 6). Zone 3, in contrast, shows a large value in subenvelope B and then follows a decreasing trend upward toward the deeper subenvelopes ( Figure 6).
[43] An important observation in this work is that zones with relatively high thermal conductivity, as Zones 1 and 3 [see Kinoshita et al., 2011], show the largest BSRs in the shallow subenvelopes A to C ( Figure 6). This is particularly the case for Zone 3, where BSRs are three times larger in subenvelope B than in the bottom half of the stability envelope. The colder Zone 2, in contrast, shows increased volumes of BSRs in the lower subenvelopes E and F ( Figure 6).
[44] This study thus proposes the concept of a stability envelope, defined by specific mathematical boundaries, to characterize the distribution of gas hydrates in the subsurface. By definition, this envelope can be divided into subenvelopes, each representing different conditions of pressure and temperature derived from empirical equations. The application of this concept to the outermost forearc and upper accretionary wedge region of Nankai results in measurable distributions of BSRs across seven subenvelopes of a BSR stability envelope (Figures 5 and 6).
[45] Based on this latter method, we propose that the ratio (R) between the surface area of upper and lower BSRs within a predefined stability envelope to indicate colder or warmer BSR traps. Values of R < 1 indicate lower hydrothermal gradients, whilst R > 1 are observed within warmer regions of the MFZ and over thrust anticlines eroded during the Quaternary, i.e., with warmer temperatures ad lower pressures due to active erosion and uplift [e.g., Hornbach et al., 2008] (Figure 7). Zone 1 shows a relatively homogeneous surface area distribution, but R > 1 as predicted by our model. In contrast, surface area measurements in Zone 2 indicate a stable state for most of the mapped BSRs, conditions that are promoted by the higher sedimentation rates recorded in Zone 2, which result in the localized cooling of the ITZ ( Figures  6 and 7). Our data also show that the majority of BSRs in Zone 3, located in structural traps of thrust anticlines that tend to a state of equilibrium with the surrounding topography , are stable under the present topographic conditions (see section 7 and Figures 2c and 7). [2011] interpreted high-amplitude reflections on seismic profiles as representing potential gas migration around structurally deformed zones. In addition, Crutchley et al. [2011] interpreted a dipping zone of high amplitude in a ridge as a phase boundary between solid clathrate and underlying free gas because of the stronger low velocity reflection. A BSR was not observed in the immediate region but was recognizable either side of the ridge at a lower depth.

Local Factors Affecting BSR Distribution
[47] In the interpreted seismic volume, the largest number of BSRs and associated seismic anomalies is observed within or around thrust anticlines. In Figures 2b and 2c, for instance, BSRs observed on the edge of thrust anticlines. However, after tracking the BSRs across the strike of the margin, we observe that the BSRs amplitude diminishes and that they shift position into dipping strata in the fold. Thus, as in Liu and Flemings [2007] data from New Zealand, we can identify dipping reflections in Nankai that are likely related to the presence of permeable layers, stratigraphic or faults and fractures focusing (free) gas into the zone of stability (Figures 2b and 2c). Blanked zones on seismic data are also observed as focused zones of fluid migrating into the BSR stability envelope [see Crutchley et al., 2011] (Figures 2b and 2c).
[48] The observation that BSRs in Zone 1 are mainly focused in a region landward from the MFZ is proposed here to relate to the advection of fluids along the MFZ from deeper strata in the subduction zone. The advection of hot fluids into Zone 1 should raise the local geothermal gradient, and focus the bulk of anomaly in the shallower subenvelopes A and B. Deeper BSRs occurring in Zone 1 are plausible evidence for paleo-BSRs or for gas hydrates currently in their early stages of migration following a change in the thermal structure ( Figure 2a). As BSR surface area decreases in the uppermost subenvelopes of Zone 2, a likely explanation to this decrease is that an equilibrated return to the background geotherm occurs on the Nankai slope with (a) increasing distance from the main heat flow anomaly induced by the MFZ, and (b) with increasing sedimentation rates in the ITZ (Figures 1c and 2b).
[49] By analyzing BSR distribution within the stability envelope, inferences can be thus made about the mechanisms controlling the formation of hydrate in Zones 1-3. Zone 1 is likely dominated by upward hydrothermal migration through the stability zone, precipitating hydrate at shallower depths. Zones 2 and 3 are inferred to be more dependent on the richness of gas supply through permeable conduits, both stratigraphically and structural, forming free gas vents. Zone 3, in particular, is recognized as a region of marked seafloor erosion associated with thrust fault uplift-a factor that locally increases subsurface temperature by allowing seawater to percolate through deeper strata and by reducing confining pressures [Pecher et al., 2005;Hornbach et al., 2008].
[50] The importance of undertaking an envelope analysis as the one in this paper is that variations in the distribution of BSRs inside the stability envelope can be used to identify different processes influencing BSR formation and subsurface fluid migration. In contrast to Kinoshita et al. [2011] work, we postulate that BSRs are essentially inside their stability envelope in Zones 2 and 3, and likely represent gas hydrates in relative equilibrium within structural and sedimentary traps. This seems to be also the case for most of the BSRs located in thrust anticlines ( Figure 2c) and in pinch-out structures on the flanks of these same anticlines (Figure 2b). Such an interpretation suggests Kinoshita et al. [2011] variations in the depth of BSRs across the Nankai slope as solely related to localized temperature and pressure effects, namely along faults focusing fluids accumulated deeper in the accretionary wedge into the surface. Careful analyses should be undertaken when assuming these local effects relate to active thrusting and erosion alone, as the anomalies are still within the stability envelope calculated for such depth ranges.

Physical Parameters Affecting the Depth of BGHSs
[51] Geological and geophysical data reveal complex heat flow dissipation processes in the Nankai accretionary wedge. For instance, Marcaillou et al. [2012] stressed the importance of variable trench sedimentation rates to thermal structure of a subducting slab. They suggested hydrothermal warming through fluid circulation near the trench to be insignificant, but of noticeable importance at basement highs and active fault regions of the arc. In contrast, variable heat-flow values were observed by Ashi et al. [1999] approximately 15-25 km landward of the deformation front. The authors proposed advective heat transfer within Nankai's accretionary wedge focused by upward fluid flow through the MFZ. Sharp variations in heat flow were found to coincide with particular chemosynthetic communities and cold seeps around splay-fault scarps [Toki et al., 2004].
[52] Recent surface deformation and depositional/ erosional processes were also seen as possible mechanisms for changes in the thermal structure of the slope around the MFZ [Ashi et al., 2002;Hamamoto et al., 2011]. Ashi et al. [2002] suggested weak or absent BSRs in slope basins, steep slopes and deep-sea canyons to be associated with rapid depositional processes, i.e., a scenario in which there was not enough time for methane gas production in infilling slope basins. On steep slopes, the erosion of sediments would result in a shift of local BGHSs downwards, preventing the accumulation of free gas and hydrates. Alternating high-permeability sands and impermeable muds can also focus or hinder the vertical migration of methane in certain regions. Variable heat flow patterns offshore Nankai were more recently documented in Harris et al. [2011Harris et al. [ , 2013. Heat flow studies along the NanTroSEIZE drill sites showed significantly lower thermal conductivity values than studies using near seafloor measurements [see Hamamoto et al., 2011]. Harris et al. [2011] put the difference in results down to variations in bottom water temperatures at water depths of around 3 km. The authors used bathymetric corrections for better confidence in heat flow determination to conclude that heat flow slightly decreases at Sites C0008A and C0008C, i.e., on the distal slope. At Site C0001, located on the prominent seafloor high adjacent to the Kumano Basin, there was a 16% increase in heat flow that contrasts with C0008A and C0008C. After correcting the model results for sedimentation processes, Harris et al. [2011] revealed that all the sites on the accretionary wedge, where known recent erosion and accumulation has occurred, yielded a 5-10% increase in heat flow values.
[53] Comparisons between the transformed empirical equations in Appendix I and the observed BSR distribution along the Nankai slope provide important insights into how certain parameters affect gas hydrate stability depths on the larger scale of the accretionary wedge. Figure I shows Dickens and Quinby-Hunt's [1994] original stability equation (equation (2)) for seawater solved for depth versus the second polynomial fit that was applied by Peltzer and Brewer [2000] (equation (3)). The polynomial fit, believed by Peltzer and Brewer [2000] to comprise a better representation for higher pressure and temperatures outside of Dickens and Quinby-Hunt's [1994] limited range, as predicted by the former authors.
[54] In Figure 4, the range of equations put forward by Lu and Sultan [2008] is shown to plot as a clustered group (equations (5) and (6)). This result confirms that variable water salinity values and pore radius have small impacts on the depth of BGHSs. In order to illustrate this latter point, Figure IIa shows two variations of equation (5) for different salinity values. A change in salinity from 2% to 3.35% had little effect on the depth of the two calculated BGHSs. Figure IIb shows three different pore radius sizes and their impact on computed BGHSs. Once again, the difference in the three computed equations, based on a known range of values for gas bearing sands [Nelson, 2009], is minimal compared to the range of equations used.
[55] Figure III displays the 3D modeling of equations (10) and (11), which are derived from studies in the Gulf of Mexico and South China Sea, respectively. Each equation used two equilibrium temperatures for the dissociation of gas hydrate (equations (7) and (8)). The use of two distinct equilibrium temperatures with a difference of %2.5 C (equations (7) and (8)) illustrates how the background temperature of the area affects the depth of stability. Higher temperatures result in thicker subenvelopes, a surprizing result as anomalous BSR observations are usually associated with local increases in temperature. This result is due to the inverse relationship between geothermal gradient and depth in the equations that govern stability thickness below the seafloor. Thus, variations in the geothermal gradient and bottom water temperatures result in significant stability depth differences between the stability surfaces ( Figure  III).

The BSR Stability Envelope Concept
[56] How does the distribution of mapped BSRs in this study compare with previous data from Nankai and other continental margins? The analysis in section 8.1 shows that differences of $2.5 Ci n background temperature can significantly affect the stability depth of gas hydrates. Temperature differences, usually occurring in association with BALE ET AL.: GAS HYDRATES' MATHEMATICAL ENVELOPE variable geothermal gradients or bottom water temperatures, also result in significant stability depth differences between modeled BGHSs.
[57] Dickens and Quinby-Hunt's [1994] original stability equation (equation (2)) for seawater shows a tendency not to increase satisfactorily with depth as much with the increase of hydrostatic pressure. In fact, the polynomial fit from equation (3) was considered as a better representation of the BGHS for higher pressure and temperatures ( Figure 3). In addition, equations (5) and (6) put forward by Lu and Sultan [2008] showed weaker variations in BSR depth than expected with changes in salinity and size of the pore radius (see Figures I, II, and III in Appendix I). Also, the effect of three different pore radius sizes, based on known ranges of values for gas bearing sands [Nelson, 2009], is minimal compared to the range of equations used.
[58] These comparisons stress that structural complexities and localized fluid flow on continental margins cannot fit in individual stability equations. This is mainly due to the fact that geothermal component of the computed stability equations only considers variations with increasing water depth, i.e., do not account for local changes in heat flow. In addition to this limitation, thrust faulting and local slope depocentres result in large variations in sediment unit thickness within the study area. Associated lithological variations will adjust the depth of stability from the assumed calculated depth or focus the formation of hydrate in the most favorable conditions, not where the calculations solve for [Harrison et al., 1982;Ginsburg et al., 2000].
[59] A significant result in this work is that each structural domain (Zones 1-3) shows distinct BSRs distributions in relation to BGHSs. BSRs thus occur within well defined subenvelopes, with BSR relative distribution and surface areas indicating if gas hydrates are in equilibrium with local topography, thermal, salinity, and depositional conditions. These subenvelopes can be mathematically modeled on 3D seismic data and the ratio R calculated by considering the relative distribution of BSRs. In the study area, the shallow BSRs within thrust anticlines of the Imbricate Thrust Zone (ITZ) relate to relatively high heat flows, as changing seawater salinity and temperatures were ruled out as controlling parameters . Deeper BSRs in these structures were interpreted to represent paleo-BSRs or regions where stable gas hydrates have not yet been established, the primary mechanism for such anomalies being the erosion of hanging-wall blocks due to active uplift over the last 10 kyr [Kinoshita et al., 2011, Figure 10]. In contrast, regions showing the larger sedimentation rates are potentially colder. Thus, when calculated across the study area, ratios  (10) and (11), and corresponding BSRs and fluid-related anomalies on the flanks of thrust anticlines. White patches correspond to interpreted BSRs. BSRs are located within the stability zone with the exception of Zone 1, where BSRs crosscut the lower boundary close to the seafloor ridge shown in Figure 2a.
[60] Our approach has the advantage of delimiting the distribution of BSRs within specific surfaces predefined by empirical and experimental equations (Figure 8). By juxtaposing and comparing attribute data to the geometry of BSRs, it will be possible to apply geological templates to model the subsurface gas hydrate distribution [e.g., Bryant et al., 2000], allowing the identification of any anomalous regions where gas hydrates are in disequilibrium. This method has also the advantage of allowing the recognition of time-dependent variations in topography, thermal conductivity, salinity, and sedimentation rates on the distribution of gas hydrates, which are reflected as changes in the position and surface area of BSRs within their stability envelope. Apart from a ''static'' analysis based on 3D seismic data, our method can be applied to gas hydrates (or other fluid-related acoustic anomalies) in 4D, using time-lapse data (Figure 8).

Conclusions
[61] Bottom-simulating reflections (BSRs) representing the transition from stable solid hydrate to free gas were mapped adjacently to a Megasplay fault branch in SE Japan (Nankai Trough). In the study area, the positions of BSRs do not conform to a single calculated Base of Gas Hydrate Stability surface (BGHS). Instead, the concept of a BGHS is extended into a BSR stability envelope, defined by multiple BGHS surfaces. This BSR stability envelope can be used to quantify subsurface BSR distribution more accurately, recognizing at the same time nonuniform stability regimes. The surfaces that define the envelope of stability represent conditions from laboratory studies or specific continental margins in which defining parameters such as geothermal gradient, pore size, and bottom water temperatures differ.
[62] The quantification in this work shows that the distribution of anomalous reflectors varies in depth with changing geological and oceanographic conditions, including [63] 1. The added influences of local heat flow variations; [64] 2. Stratigraphic discontinuities, and [65] 3. Differences in free gas sourcing and plumbing systems in an area governed by active tectonics.
[66] This paper shows that the ratio R between the volume of uppermost and lowermost BSRs in the stability envelope is of use to understand their equilibrium conditions. Ratios R > 1 indicate higher (local) heat flow, whilst R < 1 are associated with colder regions of the continental slope, seaward from the MFZ. Our data indicate the MFZ as a major heat conduit to the surface, with gas hydrates approaching equilibrium conditions the further away they are from the MFZ and uplifting structures. Higher temperatures also result in a thicker stability zone. Increases in chloride content shifts the depth of stability upward, whilst a relative increase in pore radius also decreases the thickness of the BSR stability envelope. Seismic and future borehole data will be paramount to better estimate the volumes of methane accumulated in the ITZ of Nankai, and thus correlate the results of our model with the causes of methane formation and dissociation.
[67] On the regional scale of the Nankai Trough, the results in this paper show that BSRs are in relative equilibrium within the calculated stability envelope on the Nankai continental slope. Our results indicate that the variations in the depth of BSRs across the study area  as, essentially, related to local thermal effects. They do not reflect major tectonic events that might have affected the modern Nankai slope, as the mapped BSRs are under equilibrium conditions within the stability envelope. Our analysis thus replaces ''static'' analyses of BSR distributions based on 3D seismic data, and can be applied to gas hydrates (or other fluid-related acoustic anomalies) in 4D, using time-lapse data.