The Global Inventory of Methane Hydrate in Marine Sediments: A Theoretical Approach

The accumulation of methane hydrate in marine sediments is controlled by a number of physical and biogeochemical parameters including the thickness of the gas hydrate stability zone (GHSZ), the solubility of methane in pore fluids, the accumulation of particulate organic carbon at the seafloor, the kinetics of microbial organic matter degradation and methane generation in marine sediments, sediment compaction and the ascent of deep-seated pore fluids and methane gas into the GHSZ. Our present knowledge on these controlling factors is discussed and new estimates of global sediment and methane fluxes are provided applying a transport-reaction model at global scale. The modeling and the data evaluation yield improved and better constrained estimates of the global pore volume within the modern GHSZ ( ≥ 44 × 1015 m3), the Holocene POC accumulation rate at the seabed (~1.4 × 1014 g yr−1), the global rate of microbial methane production in the deep biosphere (4−25 × 1012 g C yr−1) and the inventory of methane hydrates in marine sediments ( ≥ 455 Gt of methane-bound carbon).


Introduction
The global abundance of methane hydrate in marine sediments is still poorly constrained.Estimates are based on extrapolation of field data [1,2] and geochemical transport-reaction modeling [3][4][5][6].They range over three orders of magnitude (500-55,700 Gt of C) and a clearly constrained consensus value has not emerged over the past decades.This is a major draw-back for gas hydrate research since the resource potential and the possible impact of methane hydrates on past and future climate change cannot be evaluated without a better constrained estimate of global methane hydrate abundance.Models and observations clearly show that methane hydrate formation is basically controlled by the following key parameters: (1) accumulation of particulate organic carbon at the seafloor; (2) kinetics of microbial organic matter degradation and methane generation in marine sediments; (3) thickness of the gas hydrate stability zone (GHSZ); (4) solubility of methane in pore fluids within the GHSZ; (5) sediment compaction; (6) ascent of deep-seated pore fluids and methane gas into the GHSZ.In the following, we present a comprehensive review on the key parameters listed above and identify major knowledge gaps.Subsequently, a geochemical transport-reaction model is applied to investigate the effects of the key parameters on gas hydrate accumulation in marine sediments on a quantitative basis.The model is finally used to constrain the global gas hydrate inventory in marine sediments.

Accumulation of Particulate Organic Carbon at the Seafloor
Methane production and methane hydrate formation are fueled by the accumulation of particulate organic carbon (POC) at the seabed.The mass accumulation rate of sediments (MAR in g cm −2 yr −1 ) and the POC accumulation rate (A POC in g cm −2 yr −1 ) are calculated as: where POC 0 is the concentration of particulate organic carbon in surface sediments (in wt.%), ρ S is the density of dry solids (in g cm −3 ), Ф is porosity and w is burial velocity (in cm yr −1 ).The density of dry solids in continental margin sediments varies from 2.65 g cm −3 for quartz to 2.1-2.76 g cm −3 for different types of clay minerals.Here, we adopt ρ S = 2.5 g cm −3 as mean density of dry solids at continental margins [6].The porosity of bioturbated surface sediments (0-10 cm sediment depth) usually falls into a range of 0.7-0.9 at most continental margins.Thus, a mean value of Ф = 0.8 is applied in the following to convert burial velocities of surface sediments into mass accumulation rates [6].

Holocene POC Accumulation Rates
POC concentrations in Holocene surface sediments have been measured at several thousand sites covering most depositional areas at the seabed [7,8].The available data have been compiled and evaluated to produce a global map of POC concentrations in marine surface sediments (Figure 1).For some parts of the global seafloor POC concentrations are not available (see white areas in Figure 1).In the following, we assume a POC value of 1 wt.% for these areas.Burial velocities of Holocene surface sediments have also been determined at several hundred locations [6].The published data clearly show that burial velocities are high on the continental shelf and decrease with increasing water depth.Most of the sediments accumulating at the seabed are derived from the continents.Moreover, the solid remains of marine plankton (carbonate, biogenic opal) are deposited at the deep-sea floor as pelagic sediments.The global accumulation rate of sediments at the modern seafloor amounts to about 19.6 × 10 15 g yr −1 .This number has been derived considering the transport of riverine, eolian and ice-rafted sediments to the ocean, the accumulation of neritic carbonates at the continental shelf and the accumulation of pelagic sediments at the deep-sea floor [6].The depth-distribution of sediment accumulation and the global accumulation rate can be used to constrain parameter values of a logistics equation defining the burial velocity as function of water depth [6]: with w 1 = 0.117 cm yr −1 , w 2 = 0.006 cm yr −1 , z 1 = 200 m, z 2 = 4000 m, c 1 = 3, c 2 = 10, and z in m.Equation ( 2) yields a global sediment accumulation of 19.6 × 10 15 g yr −1 for Ф = 0.8 and ρ S = 2.5 g cm −3  when applied to the global bathymetric data set given by Menard and Smith [9].The equation predicts that most of the global sedimentation takes place on the continental shelf (14.1 × 10 15 g yr −1 at 0-200 m water depth).These estimates are consistent with other data recently presented by Baturin [10].
Baturin derived a total sedimentation rate of 18.9 × 10 15 g yr −1 with a shelf and upper slope contribution of 13.6 × 10 15 g yr −1 .Sediments accumulating at shallow water depths (<200 m) cannot contribute to gas hydrate accumulation since marine methane hydrates are only formed at significantly higher pressures i.e., water depths.During the Holocene, the accumulation of methane hydrates is thus severely limited by the continental shelf acting as a trap for the continental sediment input.Holocene POC accumulation rates are calculated at global scale applying Equations (1,2) with Ф = 0.8 and ρ S = 2.5 g cm −3 and the POC concentrations shown in Figure 1.The resulting global map shows high A POC values at the shelf of productive continental margins where POC is rapidly buried in terrigenous sediments (Figure 2).Burial velocities vary over several orders of magnitude from less than 1 cm kyr −1 in Pacific deep-sea basins to more than 100 cm kyr −1 on the continental shelf while most of the POC concentrations fall into the range of 0.1 to 1 wt.%.The A POC variability is thus mainly induced by the strong contrast in burial velocity between pelagic and continental margin settings.The global POC accumulation rate in Holocene surface sediments calculated from the distribution shown in Figure 2 amounts to 1.37 × 10 14 g yr −1 .Since continental margins are the main locus for sediment accumulation, 88% of the global POC accumulation occurs within 500 km distance from the coastline (Figure 2).Our new estimate of Holocene POC accumulation (1.37 × 10 14 g yr −1 ) is very close to a previous estimate (1.40 × 10 14 g yr −1 ) given by Baturin [10] even though Baturin's global rate is based on an independent data set and modeling approach.This surprising conformity strongly supports the validity of the new approach presented by Burwicz et al. [6].

Quaternary POC Accumulation Rates
Most of the methane produced in marine sediments does not originate from Holocene surface sediments.The thickness of these Holocene deposits rarely exceeds 10 m while the gas hydrate stability zone extends several hundred meters into the sub-surface.Thus, mean POC accumulation rates averaged over a period of several million years are more appropriate for the prediction of hydrate accumulation than Holocene accumulation rates.During the Holocene, most riverine particles are deposited on the continental shelf because the shelf is not at isostatic equilibrium with the present sea level but is still affected by the much lower glacial sea level stand [11].The sea level was 120 m below its present value during the last glacial maximum [12] reducing the water-covered shelf and marginal sea areas by more than 50% [9,13,14].Under glacial conditions, the anomalous Holocene shelf accumulation rate may have been diminished by an order of magnitude [15] shifting the focus of sedimentation from the shelf to the continental rise and slope.Moreover, the transport of ice-rafted material and the deposition of eolian dust were strongly enhanced further increasing the accumulation rates at the margin seafloor [16].Hydrates only accumulate in margin sediments deposited at more than 200 m water depth at the continental slope and rise.The mean Quaternary sedimentation rate in this environment is certainly much higher than the Holocene value.The sediment accumulation at 200-3000 m water depth would rise by a factor of 5 if most of the Holocene shelf sedimentation was shifted to larger water depth under glacial conditions as proposed by Hay [15].Quaternary burial velocities were thus calculated increasing the sedimentation rate at the continental slope and rise (>200 m water depth) at the expanse of shelf deposits [6].The global sediment accumulation rate was maintained while the slope and rise burial velocities were enhanced by a factor of 5 over a 500 km wide zone around the continental margins [6].This 500 km zone includes the continental slope and a major portion of the continental rise.The resulting burial velocities were then applied to calculate Quaternary POC accumulation rates.These calculations were based on the Holocene POC distribution map (Figure 1) since mean Quaternary POC concentrations are not available at global scale.The error involved in the map of the Quaternary POC accumulation (Figure 3) thus exceeds the error of the corresponding Holocene distribution (Figure 2).Nevertheless, methane generation and gas hydrate accumulation in the modern ocean are controlled by the Quaternary mean values (Figure 3) rather than the anomalous Holocene accumulation rates (Figure 2).The global POC accumulation rate for Quaternary surface sediments calculated from the distribution shown in Figure 3 amounts to 1.55 × 10 14 g yr −1 .In the Quaternary reconstruction, 89.5% of the global POC accumulation occurs within a distant of 500 km from the coastline.The accumulation is, however, shifted from the shelf to the continental rise and slope.

Cenozoic-Cretaceous POC Accumulation Rates
Global accumulation rates of terrigenous sediments at the seafloor reconstructed over the last 150 Myr [17] show a strong increase towards the Quaternary induced by high rates of continental erosion (Figure 4).Quaternary erosion was accelerated by the waxing and waning of continental ice shields, rapid sea-level change and intense mountain building in the Himalayas and Tibetan High Plateau.Global rates of POC accumulation at the seafloor were reconstructed using the δ 13 C record of marine carbonates [17].They feature a maximum during the Quaternary induced by high rates of sediment accumulation (Figure 4).However, the contrast to the previous periods is not as strong as observed for terrigenous sedimentation.It thus seems that average concentrations of POC in marine sediments were higher in the geological past [Equation ( 1 )].Very high POC concentrations may have prevailed prior to the Paleocene-Eocene Thermal Maximum (PETM) at 55.5 Myr.These high concentrations may have favored the built-up of a large marine gas hydrate inventory which was possibly destabilized during the PETM to induce the large negative δ 13 C excursion observed in the geological record [18].

Microbial Degradation of Organic Matter and Methane Formation in Marine Sediments
Stable carbon isotope data show that methane hydrates are usually formed from biogenic methane produced by the anaerobic degradation of organic matter in the deep marine biosphere.It is thus important to understand the mechanisms and kinetics of organic matter degradation and microbial methane formation in the marine subsurface.Organic matter deposited at the seabed is degraded by microorganisms and other benthic biota using oxygen as terminal electron acceptor.The remaining organic matter is degraded by anaerobic bacteria and archaea using dissolved nitrate and manganese (+IV) and iron (+III) minerals as oxidizing agents [19].These additional electron acceptors are usually consumed within the bioturbated surface layer (0-10 cm sediment depth).Field data show that only a small fraction of the POC raining to the seafloor is buried below 10 cm depth [20].The data reveal a marked contrast between fine-grained continental margin and deep-sea sediments [21].While at continental margins about 10% of the POC raining to the seabed is conserved and buried below 10 cm sediment depth, this fraction is reduced to about 1% at the deep-sea floor [20].Due to the very low preservation of POC in open ocean environments, gas hydrates are usually not found in pelagic sediments but only at continental margins where a significant POC fraction is buried and therefore available for microbial methane formation in the deep subsurface.
POC buried below the bioturbated surface layer is further degraded by anaerobic bacteria and archaea using dissolved sulfate as electron acceptor: where C(H 2 O) represents sedimentary POC.Microbial methane production and accumulation starts below the depth of sulfate penetration.Several steps are needed before methane is finally produced as stable end product of anaerobic microbial POC degradation.In a first step, biogenic polymers are hydrolyzed and converted into monomers.These monomers (sugars, amino acids, lipids, etc.) are then fermented into CO 2 , H 2 and a number of organic acids.Methane is finally formed by methanogenic microorganisms converting CO 2 and H 2 into methane [22]: The stoichiometry of the overall process of organic matter degradation and methanogenesis is given by the following reaction: Equation ( 5) does not imply that methane is directly formed by POC degradation.It rather defines the substrate and the final products.The intermediate steps including fermentation and CO 2 reduction with H 2 are included in the overall stoichiometry [23].
Methane formed at depth is transported upwards into the sulfate-methane transition zone (SMTZ) via molecular diffusion and advection.Within this zone, methane is oxidized by consortia of bacteria and archaea using sulfate as terminal electron acceptor [24].The overall stoichiometry of anaerobic oxidation of methane (AOM) within the SMTZ is given by: Considering this complex network of microbial processes, the key questions that need to be answered to estimate the global rate of methane formation may be formulated as: How much POC is reaching the methanogenic zone and how much of that POC can be microbially converted into methane?To answer these questions not only the POC burial rates but also the down-core changes in POC reactivity need to be defined.Since more than three decades marine geochemists and microbiologists have observed that the POC reactivity decreases strongly with sediment age and depth.In his classical paper on POC degradation kinetics in marine sediments Middelburg showed that the down-core decrease in POC reactivity can be represented by the following equation [25]: age age k age (7) where k age is the age-dependent reactivity (in yr −1 ) and age 0 is the initial age of POC (in yr).The rate of POC degradation (R POC in g C g −1 yr −1 ) is calculated as: where C POC is the POC concentration (in g C per g of dry sediment).The equation was calibrated by a large number of rate measurements showing that the reactivity of POC decreases by 10 orders of magnitude going from fresh marine plankton to POC buried in deep-sea sediments.More recently, ODP drilling confirmed the presence of living microorganisms in pelagic and terrigenous sediments down to the underlying oceanic crust [26][27][28][29].The number of prokaryotic cells (bacteria and archaea) decreases from more than 10 8 cm −3 in the upper few meter of the sediment column to less than 10 6 cm −3 at depth [26,30].The logarithmic down-core decrease in cell numbers and the microbial POC degradation rates in the deep subsurface are broadly consistent with the rates predicted by the Middelburg model [27,31,32].
The main uncertainty of the Middleburg model is associated with the initial age parameter (age 0 ) defined as the average age of POC buried below the bioturbated zone.Initial ages of 300-30,000 yrs were applied to reproduce the pore water profiles in sediment cores retrieved from the Sakhalin slope [23].An even larger range of 800-180,000 yrs was needed to reproduce pore water profiles measured in ODP sediment cores [33].These large ranges probably reflect the ambient variability in burial velocity, the deposition of refractory organic matter and the possible loss of surface sediments during core retrieval [23].The largest initial ages were observed at Blake Ridge where drift sediments are deposited containing pre-aged organic matter of terrigenous and marine origin [23,33].
Equations (7,8) can be combined and solved to estimate the POC loss within the sulfate reduction zone and thereby the input of POC into the methanogenic zone.The fraction of POC reaching the methanogenic zone results in: where t exp is the sulfate exposure time and f MI is the concentration of POC in sediments reaching the methanogenic zone divided by the POC concentration in sediments buried below the bioturbated zone at 0.1 m depth.The sulfate exposure time depends on sulfate penetration depth (z S in m) and burial velocity (w in m yr −1 ): Burial velocities at the deep-sea floor typically range from 0.1 to 10 cm kyr −1 [Equation ( 2)] while sulfate penetrates usually more than 100 m into pelagic sediment.A compilation of sulfate penetration depths in open ocean sediments confirmed that sulfate is usually not fully depleted in these pelagic deposits [32].This observation is partly explained by the deep penetration of dissolved oxygen into pelagic sediments [29] and the delivery of sulfate through the base of the sediment column via sulfate-rich fluids circulating through the underlying oceanic crust.Deep-sea sediments are thus typically exposed to sulfate over their entire life span (up to ~150 Myr).Microbial rate measurements and pore water profiles obtained in pelagic sediments revealed that methanogenesis may occur already within the sulfate-bearing sediment horizons [27].Methane produced by this process is, however, consumed within the sulfate-bearing zone via AOM.Methane concentrations in pelagic sediments are therefore several orders of magnitude lower than the ambient saturation value with respect to methane hydrate [34].Due to this strong undersaturation of pore fluids, methane hydrates have never been documented in pelagic sediments.
In most continental margin settings, sulfate dissolved in ambient pore fluids is completely consumed at 10-100 m sediment depth [32,35].The sulfate exposure time results as t exp = 10 4 -10 6 yr for burial velocities of 10-100 cm kyr −1 .Thus, about 30-90% of the buried POC is reaching the methanogenic zone in continental margin sediments (Figure 5).The microbial degradation of POC delivered to the sulfate-free methanogenic zone may lead to the accumulation of methane in pore fluids and to the formation of methane hydrates pending on the rate of methane formation and the local pressure and temperature conditions.9) for initial ages of 1000 yrs, 10,000 yrs, and 100,000 yrs.Dissolved sulfate profiles often show an almost linear decrease with sediment depth suggesting that only a minor sulfate fraction is used for POC degradation while most of the sulfate is consumed via AOM in the underlying sulfate-methane transition zone [36,37].Equation (9) implies that only 10% of the POC buried in continental margin sediments is consumed via microbial sulfate reduction where pre-aged organic matter (age 0 = 10 5 yr) is deposited at the seabed (Figure 5).It is thus likely that linear sulfate profiles are induced by the deposition of old POC at the seabed.Terrestrial POC delivered to the margins by continental erosion is a likely source for this pre-aged material.However, sulfate profiles with significant curvature documenting strong sulfate depletion via anaerobic POC degradation are also commonly observed [23].They probably reflect the deposition of young marine organic matter since Equation (9) implies that up to 80% of this young POC (age 0 = 10 3 yr) is consumed via microbial sulfate reduction (Figure 5).
The kinetics of microbial methane formation in marine sediments was studied in detail by Wallmann et al. [23].Evaluating pore water and sediment data from the Sakhalin Slope and the Blake Ridge, the authors showed that the rate of methane formation within the methanogenic zone (R CH4 in g of methane carbon g −1 yr −1 ) can be calculated as:  (11) where K C is an inhibition constant (K C = 30-50 mM), C DIC is the concentration of dissolved inorganic carbon in ambient pore fluids and C CH4 the dissolved methane concentration (in mM).The factor 0.5 defines the fraction of POC being converted into methane.The remaining POC fraction is decomposed into carbon dioxide for sedimentary POC with an oxidation state of carbohydrates [Equation (5)].This extension of the Middelburg model was subsequently tested and applied at a number of ODP drill sites in prominent gas hydrate provinces [33].The evaluation of pore water profiles showed that kinetic Equation (11) gives indeed a good approximation of the down-core changes in microbial methane formation within the GHSZ.Equation (11) can be solved to calculate the fraction of POC being converted into methane carbon within the methanogenic zone: where t is the residence time of sedimentary POC within the methanogenic zone while age 0 is the initial age of POC entering the methanogenic zone from above.The organic matter entering the methanogenic zone of typical margin sediments has initial ages of 10 4 -10 6 yr considering ambient sedimentation rates and sulfate exposure times.Methane formation is further inhibited by the accumulation of dissolved metabolites (DIC and CH 4 ) within ambient pore fluids.Assuming typical concentrations of C DIC = C CH4 = 50 mM and an inhibition constant of K C = 40 mM, Equation ( 12) can be applied to calculate the POC fraction being converted into methane within the sulfate-free methanogenic zone (Figure 6).The residence time of sediments within the methanogenic zone is limited by the geothermal gradient since microbial methane formation cannot proceed at temperatures above ~100 °C.Assuming a geothermal gradient of 30-40 °C/km, the deep biosphere may extend to a sediment depth of about 3 km.At continental margins with typical burial velocities of 10-100 cm kyr −1 , the residence time within the methanogenic zone may thus fall into a range of t = 3-30 Myr. Figure 6 and Equation (12) show that about 10%-20% of the POC entering the methanogenic zone is microbially converted into methane under these assumptions (for age 0 = 10 5 yr).Since 30-90% of POC buried below the bioturbated surface zone of continental margin sediments is reaching the methanogenic zone (Figure 5), about 3%-18% of the buried POC is converted into methane by the deep biosphere residing in underlying sulfate-depleted sediment layers.
Considering the arguments and equations presented above it is now possible to estimate the rate of microbial methane formation in continental margin sediments at global scale (Table 1).This new estimate of microbial methane production in margin sediments (4-25 Tg C yr −1 , Table 1) is significantly lower than previous estimates of the global microbial methane production ranging from 60 Tg C yr −1 [38] to 240 Tg C yr −1 [39].These previous estimates approach or even exceed the global rate of POC burial below the bioturbated zone.Our best estimate of ~10 Tg C yr −1 would lead to an accumulation of up to 100,000 Gt CH 4 -C over an accumulation period of 10 Myr if none of the methane would be lost by upward diffusion, AOM, fluid flow and gas ascent.However, field data clearly show that a large fraction of methane is lost by upward diffusion and fluid flow into the overlying sulfate-bearing zone where methane is consumed by microbial consortia using sulfate as terminal electron acceptor [24].The biogenic methane inventory is thus certainly much smaller than the time-integrated rate of biogenic methane formation derived above.Nevertheless, the number of 100,000 Gt C is a useful upper limit for biogenic methane accumulation.Interestingly, the global methane hydrate inventory estimated by Klauda and Sandler [4] approaches this upper limit value.
The kinetic models of Middelburg and Wallmann can also be used to estimate the total POC degradation rate within the deep marine biosphere.The Middelburg model [Equation ( 8)] yields a rate of 1.2 × 10 14 g yr −1 applying a mean POC input flux of 1.4 × 10 14 g yr −1 , an initial age of 10 4 yr and a POC residence time in the deep biosphere of 10 Myr.It predicts that almost 90% of the buried POC would be consumed within the deep biosphere.Applying the same parameter values, an inhibition constant of 40 mM and mean methane and CO 2 concentrations of 50 mM, the Wallmann model [Equation (11)] predicts a POC consumption rate of 0.7 × 10 14 g yr −1 corresponding to ~50% of the burial flux.The remaining POC is further decomposed below the deep biosphere via chemical processes yielding thermogenic gas, oil, kerogen and other products of high-temperature (>100 °C) thermal maturation.It should, however, be noted that the numbers given above are very rough estimates since microbial organic matter degradation in the deep biosphere is a very complex and poorly understood process.

Thickness of the Gas Hydrate Stability Zone
Methane hydrates are only stable at low temperatures and high pressures.The stability field of structure type-I methane hydrate is well defined [40].Pitzer equations can be used to constrain the effects of seawater salinity and porewater composition on methane hydrate stability [34].The stability field for sulfate-free seawater with a salinity of 35 PSU is shown in Figure 7.The sharp phase boundary depicted in Figure 7 is however not strictly valid for sediments and other porous media.Capillary forces acting on gas hydrates and free gas residing in sediment pores of different sizes create a broad transition zone where free gas and gas hydrate coexist [41].The thickness of this zone has been estimated as 20-28 m for Hydrate Ridge and Blake Ridge [41].The seismic bottom simulating reflector (BSR) indicating the top of the free gas occurrence zone is therefore often located at shallower depths than predicted by the bulk phase boundary.Moreover, the composition of the hydrate-forming natural gas has a strong effect on the positioning of the phase boundary.Biogenic gas formed by the microbial degradation of organic matter is composed of methane and CO 2 and contains only trace amounts of higher hydrocarbons (C 2+ ).CO 2 hydrates are usually not formed in marine sediments since their solubility clearly exceeds the CO 2 concentration in ambient pore fluids.Methane hydrates are less soluble and therefore the major product of biogenic gas formation in marine sediments.The hydrate-bound gas thus typically contains more than 99.9% methane with only trace amounts of CO 2 and C 2+ [42].The low abundance of higher hydrocarbons in biogenic gas supports the formation of structure type I methane hydrate.However, thermogenic gas ascending from larger sediment depths contains significant amounts of C 2+ favoring the formation of structure type II and/or the inclusion of C 2+ -components in structure type I methane hydrate.These hydrates are stable over a significantly broader P-T range and thus induce a down-core displacement of the phase boundary [40].The thickness of the gas hydrate stability zone (GHSZ) in marine sediments is usually estimated considering the temperature (T) and pressure (P) conditions at the seabed and the increase in T and P with sediment depth.High-resolution data sets are available to define temperatures and pressures at the seafloor on global scale [6] while the conductive heat flux has been measured at more than 9000 sites [43].Most of these heat flow data have been generated by measuring the temperature gradient within the upper few meters of the sediment column.The conductive heat flux (q H ) is then calculated from the temperature gradient considering the ambient thermal conductivity of bulk sediments (λ B ): Heat flow data rather than temperature gradients (dT/dz) have been compiled at global scale [43].Geothermal gradients defining the thickness of the gas hydrate stability zone, thus, need to be derived from the heat flux data applying appropriate estimates for ambient thermal conductivity.Thermal conductivities of bulk sediments quantify the efficiency of heat transport which involves transport (1) from grain to grain; (2) from grain to liquid to grain; and (3) through pore-filling liquid [44].Rather than calculate the contribution of each heat transport path explicitly, thermal conductivity is often estimated using a two-phase mixing model to combine the thermal conductivities of the sediment grains with the pore fluid.The thermal conductivities of methane hydrate and water differ by <10% at the temperatures found in hydrate-bearing sediments [44].For this reason, first-order thermal conductivity estimates can neglect the presence of methane hydrate and assume that the sediment pore space contains only water [45].The following two-phase mixing model provides a reasonable approximation for thermal conductivity in the uppermost 500-1000 m of the sediment column [45]: where λ B , λ W , and λ S are thermal conductivities of bulk sediment, seawater and solids, respectively.The thermal conductivity of seawater increases with pressure and temperature [46].It varies from 0.56 to 0.60 W m −1 K −1 over the P-T range encountered within the GHSZ [46].A mean value of λ W = 0.58 W m −1 K −1 is thus representative for the marine hydrate stability zone.Quartz has a conductivity of 7.7-8.4W m −1 K −1 [44] while calcite has an average conductivity of 3.6 W m −1 K −1 [47].Thermal conductivities of clay minerals vary from 5.2 W m −1 K −1 for chlorite to 1.8-2.7 W m −1 K −1 for smectite, illite and kaolinite [47].Data measured in ODP cores retrieved at various continental margin sites indicate a mean value of λ S = 2.5 W m −1 K −1 for the dry solids in fine-grained margin sediments [45].The average porosity within the upper few meter of the sediment column is typically close to 0.75 [45].Applying Equation ( 14), the thermal conductivity of clay-rich surface sediments results as λ B = 0.836 W m −1 K −1 for λ W = 0.58 W m −1 K −1 and λ S = 2.5 W m −1 K −1 .The geothermal gradient in surface sediments is thus derived by dividing the ambient heat flux data by λ B = 0.836 W m −1 K −1 [Equation (13)].
Thermal conductivity changes, however, with sediment depth due to the compaction of sediments inducing a down-core decrease in porosity.Porosity change is often described by the following equation [19]: where Φ 0 is the porosity at the sediment surface, Φ f the porosity at the base of the sediment column and p is an attenuation coefficient.The mean porosity profile observed in ODP cores can be approximated applying Φ 0 = 0.74, Φ 0 = 0.3 and p = 1/600 m −1 [48].
The steady-state temperature profile in marine sediments is usually controlled by heat conduction.Neglecting heat production and convection, it can be described by the following differential equation [49]: This equation can be solved applying the following boundary conditions at the top of the sediment column at z = 0: where T Grad is the temperature gradient recorded in surface sediments, q H is the heat flow calculated from the gradient and the thermal conductivity of surface sediments while T BW is the temperature of ambient bottom water.The solution of this equation for constant λ B results as: The equation can also be solved for variable λ B considering the down-core changes in porosity and λ B defined in Equations (14,15).The resulting temperature profile is depicted in Figure 8 for typical continental margin sediments (solid black line in the right panel).It shows a clearly non-linear shape since the down-core decrease in porosity induces a significant λ B increase with depth (left panel in Figure 8).The linear profile calculated from the λ B value in surface sediments [λ B = 0.836 W m −1 K −1 , red dotted line in the left panel, Equation ( 16)] overestimates sediment temperatures at depth and thus underestimates the thickness of the GHSZ by almost 150 m (Figure 8).It is thus important to consider the non-linear nature of sedimentary temperature profiles.The temperature profile and GHSZ thickness can, however, be approximated applying a constant λ B of 1 W m −1 K −1 (blue dotted-broken line in right panel of Figure 8).This mean value is thus often applied to estimate the vertical extent of the GHSZ in marine sediments [45].
Global maps of the GHSZ thickness were produced by Burwicz et al. [6] and Pinero et al. [50] assuming constant temperature gradients and applying the same global data sets for water depth, bottom water temperature, heat flow and sediment thickness.The base of the GHSZ is defined by the intersection of the sedimentary temperature profile and the phase boundary for sulfate-free pore water (Figure 8).Pinero et al. [50] assumed a constant thermal conductivity of λ B = 1.0 W m −1 K −1 over the entire ocean.The map produced under this assumption gives a good approximation for the GHSZ thickness in fine-grained continental margin sediments (Figure 9).The vertical extend of the GHSZ is, however, underestimated where sediments are dominated by sand, chlorite and calcite since sedimentary bulk conductivity is strongly enhanced by these components [47].Moreover, the map is based on a thermodynamic model valid for pure methane hydrate of structure type I [34].It may thus underestimate the vertical extent of the GHSZ where methane hydrate is stabilized by C 2+ components entrained into the GHSZ by the ascent of thermogenic gas (Figure 7).Burwicz et al. [6] applied a constant λ B value of 1.5 W m −1 K −1 which over-predicts the GHSZ thickness in margin sediments composed of smectite, illite and kaolinite but may give a better approximation for pelagic sediments.Both maps show an extended and deep-reaching GHSZ at high latitudes and around passive continental margins whereas a thin GHSZ is derived for the open ocean where the thickness of the GHSZ is limited by the thin sedimentary cover.15)], thermal conductivity of bulk sediment [λ B , Equation ( 14)] and temperature [T, Equation ( 16)] in typical continental margin sediments (Φ 0 = 0.74, Φ f = 0.3 and p = 1/600 m −1 , q H = 0.03 W m −2 , λ W = 0.58 W m −1 K −1 , λ S = 2.5 W m −1 K −1 ).The phase boundary for methane hydrate (structure type I) shown in the right panel (black broken line) is calculated for sediments deposited at 2000 m water depth, hydrostatic pressure conditions and sulfate-free sediment porewater with a salinity of 35 [34].The global volume of the GHSZ is estimated to be 67 × 10 15 m 3 integrating the data shown in Figure 9.This value should be considered as a minimum estimate since the thickness of the GHSZ may be underestimated in those areas where solid phases with high thermal conductivity and C 2+ -bearing gas hydrates are deposited.The pore space available for gas hydrate accumulation amounts to 44 × 10 15 m 3 considering the decay of porosity over the sediment column [Figure 8, Equation (15)].Up to 4 × 10 21 g CH 4 -C could reside in the GHSZ if the entire pore space (44 × 10 15 m 3 ) would be filled by methane hydrate.This maximum estimate clearly exceeds the total methane pool generated via microbial POC degradation over a period of 10 Myr (1.0 × 10 20 g CH 4 -C, see section 3).Thus, only 2.5% of the global GHSZ pore space could be occupied by gas hydrates if the entire biogenic methane pool would be fixed as gas hydrate.

Solubility of Methane in Sedimentary Pore Fluids
Methane produced by the microbial degradation of organic matter is dissolved in ambient pore fluids.The solubility of methane in water is however limited by the low polarity and hydrophobic nature of methane molecules.Excess methane formed by pervasive methane production thus precipitates as gas hydrate or forms a free gas phase pending on the local P-T conditions.Methane hydrate precipitates when the methane concentration in the pore fluids exceeds the solubility of methane hydrate.The equilibrium concentration with respect to methane hydrate is restored by hydrate precipitation extracting excess methane from solution.The solubility of methane hydrate has been measured experimentally and calculated using various thermodynamic approaches [44].It increases with temperature and decreases with pressure and salinity whereas the solubility of methane gas in water is enhanced under high pressure and reduced by an increase in temperature and salinity.The down-core rise in sediment temperature thus promotes an increase in dissolved methane concentrations within the GHSZ and a decrease in dissolved methane in the underlying free gas zone (Figure 10).Capillary forces affect the stability and solubility of gas hydrate and gas in fine-grained marine sediments.A transition zone is formed where the solubility ranges of gas bubbles and hydrate crystals residing in sediment pores of different sizes overlap [41].The discontinuity at the base of the GHSZ (Figure 10) is thus replaced by a smooth and continuous transition in methane solubility [41,44].

Sediment Compaction, Fluid Flow and Gas Ascent
Gas hydrates are formed from methane being produced within the GHSZ (section 3) and by methane ascending into the GHSZ from below.Upward flow of methane is induced by numerous processes including sediment compaction, tectonic over-pressuring and the buoyancy of methane gas and methane charged fluids.Sediment compaction is driven by the overburden load of younger sediments being continuously deposited at the seabed [48].In most marine settings, it induces an exponential decrease in porosity (Φ) and water content with sediment depth (Figure 8).Therefore, exponential equations such as Equation ( 15) are used to describe the compaction-driven down-core decrease in porosity.Applying Equations (1,15), the burial velocity of pore water (v B ) and solids (w B ) are calculated as [19]: The equations above are valid for steady state conditions, i.e., for homogenous deposits formed by continuous sedimentation.The porosity of marine sediments deposited on oceanic crust is usually not completely eliminated by compaction.The residual porosity at the base of the sediment column (i.e., the crust-sediment interface) is given by Φ f while w f is the burial velocity of water and solids which is attained after the porosity has been reduced to Φ f via steady-state compaction.Compaction promotes the burial of solids while pore water burial is diminished since water is expelled and transported upwards due to the down-core reduction in pore space (Figure 11).However, steady-state compaction does not induce a flow of water into the overlying water column: The pore fluids do not escape into the bottom water since the pore water burial velocity induced by the continuous deposition of young sediments at the seabed exceeds the velocity of the compaction-driven upward fluid flow.Seepage of pore fluids is frequently observed at continental margins where sediments are compressed by tectonic forces [52] since marine sediments subducted below continental crust and accreted at active continental margins are over-pressured and dewatered to such an extent that pore fluids escape into the bottom water [53][54][55][56].Salt diapirism, sediment erosion and slope failure induce additional fluid seepage at passive and active continental margins.Upward fluid flow is ultimately driven by excess pressures as described by the Darcy equation [57]: where the Darcy flow u (m 3 m −2 yr −1 ) is the volume of water (m 3 ) being transported through a cross-sectional sediment area (m 2 ) over a period of time (yr).Upward fluid flow occurs when the sedimentary pressure gradient ( z P   in Pa m −1 ) exceeds the hydrostatic pressure gradient defined as product of water density (ρ W in kg m −3 ) and gravitational acceleration (g = 9.81 Pa m 2 kg −1 ).The velocity of fluid flow is modulated by the intrinsic permeability of bulk sediments (k in m 2 ) and the viscosity of water (µ w ≈ 10 −3 Pa s).The intrinsic permeability of marine sediment varies over several orders of magnitude.It increases with porosity and grain size: Terrigenous marine sediments containing fine-grained clay have typical permeability values of 10 −18 -10 −14 m 2 while clean sands have values of 10 −12 -10 −10 m 2 [57].
The sedimentary pressure gradient is enhanced in compressive tectonic settings and may approach the lithostatic pressure gradient defined as the product of g and the bulk density of sediments.The bulk density of sediments  B is calculated as: The Darcy flow at lithostatic pressure u L is thus calculated as: The inter-granular fluid flow velocity v is related to the Darcy flow velocity u as: Both velocities, u and v, attain a positive sign for a downward movement with respect to the sediment-water interface while a negative sign indicates the ascent of pore fluids towards the sediment surface.Fluids ascending through a constantly growing homogenous sediment deposit subject to steady-state compaction have an inter-granular fluid flow velocity of: Fluids are thus expelled across the sediment-water interface into the overlying water column when the velocity of fluids infiltrating the sediment column from below (v) exceeds the burial velocity of pore fluids (v B ). Upward fluid flow velocities (v) have been measured at various cold seep sites and mud volcanoes located at active and passive continental margins [53,54,[58][59][60][61][62][63].They typically range from 0.1 to 100 cm yr −1 .These flow velocities are consistent with u L values predicted for over-pressured clay-bearing sediments deposited at continental margins (Figure 12).Gaseous methane occurs below the GHSZ and in the transition zone at the base of the GHSZ.The bottom simulating reflector (BSR) observed in many seismic surveys marks the boundary between hydrate-bearing sediments and the underlying gas-bearing sediment zone.The density of methane gas is lower than the density of pore fluids.Therefore buoyancy forces may induce an upward flow of gas into the GHSZ.For sediments containing both gas and water, Darcy's law is formulated as [57]: where the subscripts W and G indicate water and gas, respectively.Relative permeability appears as additional non-dimensional parameter in the equations above (k RW and k RG ).Various equations have been proposed to quantify the relative permeability of water and gas in marine sediments at different saturation values [44,57].One of the more commonly used forms is given by [64]: where m is the pore-size-distribution index.The normalized water saturation S WN is defined as: where S W is the water saturation, i.e., fraction of pore space occupied by water (S W = 1 − S G ).The irreducible water saturation has a typical value of S WC = 0.3 while the corresponding gas saturation falls into the range of S GC = 0.01-0.1 [65].The irreducible gas saturation is the most critical parameter for the prediction of upward gas migration in marine sediments since gas migration occurs only when the gas saturation S G exceeds S GC .A high S GC value of 0.1 indicates that gas migration is inhibited if less than 10% of the pore space is occupied by gas [S G < 0.1, S WN = 1, k RG = 0, q G = 0, s.Equations (25)(26)(27)].This would imply that gas migration across the BSR is effectively suppressed since the methane gas produced by microbial POC degradation below the GHSZ rarely builds up to saturations of S G > 0.1.The BSR is produced by the strong seismic contrast between gas below the GHSZ and gas hydrate within the GHSZ.Seismic surveys show a wide-spread occurrence of the BSR indicating that a significant fraction of methane gas is not mobile but retained and buried in sediments underlying the GHSZ.However, seismic chimneys indicating the ascent of free gas are also frequently observed at continental margins and are often associated with high gas hydrate saturations in the overlying GHSZ.It thus seems that the free gas is usually immobile but locally mobilized and transported into the GHSZ.The controls on gas migration into the GHSZ are poorly understood and this lack of physical understanding causes a large degree of uncertainty in the prediction of gas hydrate accumulation in marine sediments.
The gas pressure P G is related to the water pressure P W as [57]: where P C is the capillary pressure.P C may be calculated as [66]: where P CC has a value of 13,300 Pa in low permeability sediments while npc takes a value of 1.7 [66].
Darcy gas flow at constant gas saturation ( 0 ) and hydrostatic pressure ( ) can be calculated as function of gas saturation (S G ) and intrinsic permeability to quantify the buoyancy driven gas ascent into the GHSZ: Figure 13 shows this buoyancy driven flux for typical continental margin sediments deposited at 2000 m water depth.Applying a geothermal gradient of 30 °C/km, the temperature at the base of the GHSZ is 19.1 °C while the hydrostatic pressure at the phase boundary amounts to 257 bar (Figure 8).The densities of methane gas and water are ρ G = 198 kg m −3 [67] and ρ W = 1036 kg m −3 [68] under these P-T conditions while the dynamic viscosity of methane gas amounts to µ G = 11 × 10 −6 Pa s [69].The gas flow is completely inhibited at S G < S GC and increases rapidly with increasing S G at S G > S GC .Buoyancy forces would thus induce a rapid ascent of gas into the GHSZ if the gas saturation exceeds the critical value.Under these conditions the methane hydrates buried and dissociated below the GHSZ would not be lost from the system but would be constantly recycled into the GHSZ via upward gas migration.The methane hydrate inventory would expand over time until the pore space available for gas migration is clocked by methane hydrate.Very high methane hydrate saturations could be attained by this recycling mechanism.It is important to note that gas ascent is a much more efficient mechanism for the transport of methane into the GHSZ than fluid flow.For the conditions shown in Figure 13, an upward Darcy gas flow of 1 cm yr −1 introduces a methane flux of 1980 g m −2 yr −1 (for ρ G = 198 kg m −3 ) while a corresponding water flow of 1 cm yr −1 induces a methane flux of only 24 g m −2 yr −1 at a dissolved methane concentration of 150 mM (Figure 10).
The POC accumulation rates that are fueling methane production at depth rarely exceed 10 g m −2 yr −1 (Figure 3).The entire methane pool potentially formed by the microbial and thermal decomposition of buried POC would thus be transported into the GHSZ at gas and water flow velocities as low as u G ≈ 0.002 cm yr −1 and u W ≈ 0.2 cm yr −1 .
A better understand of the physical controls on gas and water ascent into the GHSZ is thus a high priory requirement for the prediction of gas hydrate accumulation at local and global scale.

Controls on Gas Hydrate Accumulation in Marine Sediments
In the following section we will present a series of model runs obtained with the transport-reaction model introduced by Wallmann et al. [23] to investigate the controls on gas hydrate accumulation.We will start this analysis by presenting the results of a reference simulation.The sensitivity tests presented in the following sections are based on this reference case.

Reference Case
The reference case presents a continental margin site at about 2000 m water depth with moderately high concentrations of POC and very high sedimentation rates.Porosity and temperature profiles applied in the standard case are shown in Figure 8, the corresponding methane solubility values are depicted in Figure 10 while the burial velocities of solids and pore fluids are shown in Figure 11.Parameter values used in this reference model run are listed in Table 2.The model simulates the reactive transport of three dissolved species (DIC, sulfate, methane) and two solid species (POC and methane hydrate).Molecular diffusion, burial and compaction are considered as transport processes.The simulated reactions include POC degradation via sulfate reduction, AOM, methanogenesis, gas hydrate formation and gas hydrate dissolution.Fixed concentration values are applied at the upper boundary (Table 2) while zero gradients are used as lower boundary conditions for all dissolved species.The model domain excludes the upper bioturbated zone of the sediment column and includes the transition zone at the base of GHSZ where the saturation concentration of dissolved methane is assumed to be constant over a depth interval of 23 m (see section 4) [41].All model results shown in the following figures and tables are valid for steady state conditions which were reached after a simulation time of ~3 Ma.The reference model does not consider gas and fluid ascent into the GHSZ from below.
The dissolved inorganic carbon concentrations calculated in the reference run show a continuous increase with sediment depth reflecting the microbial degradation of POC (Figure 14).Dissolved methane production starts at ~8 m sediment depth after sulfate is almost completely consumed by organic matter degradation and AOM (Figure 15).Dissolved methane increases down-core and reaches the saturation value with respect to methane hydrate at a sediment depth of ~86 m.At this point gas hydrate starts to form in the sediment column (Figure 14).Methane hydrate formation continues since the microbial production of methane within the GHSZ over-compensates the down-core increase in methane hydrate solubility.Close to the base of the GHSZ at ~545 m sediment depth, methane hydrate dissolves since the saturation level is no longer maintained by the sluggish methane production from aging POC.The saturation level is restored by methane hydrate dissolution in the deeper portion of the GHSZ (545-567 m) and by a combination of dissolution and dissociation in the underlying transition zone (567-590 m).Methane hydrate buried below the model domain to sediment depths >590 m dissociates to form dissolved and gaseous methane.The modeling shows that only ~22% of the POC being deposited below the bioturbated zone (729 mmol m −2 yr −1 , Table 3) is degraded within the GHSZ.The remaining portion is buried to larger sediment depths.Sulfate reduction contributes ~26% to the total POC degradation while most of the POC degradation (~74%) occurs within the methanogenic zone of the GHSZ.Sulfate diffusing into the sediment from above (52.8mmol m −2 yr −1 ) is consumed via AOM (~60%) and POC degradation (~40%).The slight curvature of the sulfate pore water profile (Figure 15) is induced by sulfate reduction via POC degradation occurring within the top ~8 m of the sediment column.Methane produced via POC degradation within the GHSZ (60.4 mmol m −2 yr −1 ) is consumed by AOM (~52%) and buried below the GHSZ (39%).Only a small portion (~9%) is permanently fixed in gas hydrates (5.5 mmol m −2 yr −1 , Table 3).The resulting gas hydrate saturations are very low.The maximum saturation of 0.35 vol-% is reached at ~470 m sediment depth.At this depth 99.65% of the pore space is filled by water while only 0.35% is occupied by methane hydrate.The depth-integrated inventory of carbon bound in methane hydrate amounts to 62 kg m −2 .Most of the carbon resides in the sedimentary organic carbon pool and in the pore water as dissolved methane and DIC (Table 3).In the following Sections (7.2-7.6)selected parameter values of the reference simulation are varied to explore the controls on gas hydrate accumulation in marine sediments.

Effects of POC Concentration, Burial Velocity and GHSZ Thickness on Gas Hydrate Accumulation
Model runs were executed to study the sensitivity of gas hydrate accumulation to the POC concentration and the burial velocity at the upper boundary of the model (Figure 16).Methane hydrate accumulation within the GHSZ is enhanced under high POC since more substrate is available for microbial methane production (left panel in Figure 16).The relationship between gas hydrate accumulation and burial velocity (w 0 ) is more complicated.At w 0 values of 10-150 cm kyr −1 , gas hydrate accumulation increases with w 0 since the supply of POC into the methanogenic zone is promoted by increasing w 0 values.Moreover, POC reactivity is enhanced due to the decrease in burial age.However, a further increase in burial velocity does not induce an expansion of the gas hydrate inventory (Figure 16).A decline in the gas hydrate inventory is observed at extremely high burial velocities (>150 cm kyr −1 for 2 wt.%POC; >250 cm kyr −1 for 1 wt.%POC).These trends are induced by the accelerated burial of gas hydrate and dissolved methane below the GHSZ.They are further supported by the decrease in the residence time of POC within the GHSZ which reduces the time-integrated rate of methane production within the GHSZ.It is important to note that model runs with identical values for the POC accumulation rate may yield very different gas hydrate inventories.Thus, the hydrate inventory increases to 427 kg m −2 at POC = 5 wt.% and w 0 = 20 cm kyr −1 while a value of only 154 kg m −2 is attained at POC = 2 wt.% and w 0 = 50 cm ky −1 even though the POC accumulation rate amounts to 6.5 g m −2 yr −1 in both model runs [Equation ( 1 )].The POC accumulation rate is thus not a useful master parameter for the calculation of gas hydrate inventories especially at high POC concentrations and burial velocities.An expansion of the GHSZ induces an exponential increase in gas hydrate inventories since the prolonged residence time of sediments within the GHSZ allows for more microbial methane formation and gas hydrate accumulation (Figure 17).Results of the individual model runs (Figures 16 and 17) were analyzed to derive a simple transfer function:  (31) where GHI is the inventory of methane hydrate within the GHSZ (in kg methane-C m −2 ), L GHSZ is the thickness of the gas hydrate stability zone (in m) while POC 0 (in wt.%) and w 0 (in cm kyr −1 ) are the POC concentration and the burial velocity of solids in surface sediments.Parameter values of the equation were determined by a minimization procedure using the FindMinimum object in MATHEMATICA 8.0.The best fit to the individual model runs depicted in Figures 16 and 17 was obtained with (parameter values ± standard deviations): a = 0.003171 ± 0.00042; b = 1.862 ± 0.020; c = 54.83 ± 6.9; d = 1.092 ± 0.039; e = 1.183 ± 0.022; f = 0.01537 ± 0.0043.The results of the more than 100 model runs are well reproduced by the transfer function (Figure 18).The error in the gas hydrate inventory introduced by the transfer function may be quite significant at low GHI values but is always <10% at GHI ≥ 200 kg m −2 (see residuals in the right panel of Figure 18).The function is valid for the parameter values applied in the reference simulation (Table 2).It can be used to estimate GHI from ambient L GHSZ , w 0 , and POC 0 values for GHI ≥ 0 over a range of L GHSZ = 0-1000 m; w 0 = 10-600 cm kyr −1 , and POC 0 = 0.5-5 wt.%.It is more reliable than our previously published transfer function using POC accumulation rates as master variable for the prediction of gas hydrate inventories [33].
Controls on the gas hydrate inventory (GHI) are shown in Figure 19.GHI increases with POC 0 and L GHSZ .Maximum values of >1000 kg m −2 are obtained at POC 0 > 3.5 and L GHSZ > 750 m.GHI also increases with burial velocity at w 0 = 10-100 cm kyr −1 while a further w 0 increase beyond 100 cm kyr −1 has no significant effect on gas hydrate accumulation.GHI values are zero for POC 0 < 0.5 wt.%.

Effects of Sediment Compaction on Gas Hydrate Accumulation
The accumulation of gas hydrates is also strongly affected by sediment compaction since the burial velocity of solids and pore water is not only controlled by the mass accumulation rate but also by the down-core decrease in sediment porosity.Compaction accelerates the velocity of sediment and gas hydrate burial while the burial of pore fluids is diminished (Figure 11).Gas hydrate accumulation is favored by a decrease in pore water burial since less dissolved methane is lost to the sediments underlying the GHSZ while a more rapid burial of solids decreases the residence time of POC and gas hydrate within the GHSZ and thereby the accumulation of gas hydrates.However, compaction also increases the solid to water ratio within the GHSZ.Thus, for a given unit volume of bulk sediment more methane is released from the decaying POC while less methane is needed to saturate the pore fluids.The overall effect of compaction is thus a strong increase in the gas hydrate inventory within the GHSZ (Figure 20).
The porosity data measured within the GHSZ give some constraints on the degree of compaction.However, different interpretations of the porosity profile are possible since the base of the sediment column is usually located below the GHSZ.Sedimentary deposits at continental margins are often several kilometers deep such that the base of the sediment column is not reached by conventional scientific drilling.The porosity at the base of margin deposits is thus poorly constrained.It may be significantly lower than the value assumed in the reference simulation (Φ f = 0.3). Figure 21 shows that a complete compaction at large sediment depth is also consistent with the porosity observations within the GHSZ at continental margins.Complete compaction (Φ f = 0) is however inhibiting the burial of pore fluids [Equation (19) and Figure 21] and thus greatly enhancing the accumulation of gas hydrates within the GHSZ.Complete compaction may occur at the base of thick sedimentary deposits and in compressive tectonic settings.This situation may be more representative for many continental margin sites than the conditions applied in the reference simulation.Therefore, an additional set of model runs was performed to explore gas hydrate accumulation for a complete compaction scenario applying the parameter settings defined in Table 2 with Φ f = 0.0, p = 1/1300 (Figure 21).The simulations show that the gas hydrate inventories are significantly enhanced by complete compaction (Figure 22).In these simulations dissolved methane is only lost by upward diffusion into the AOM zone while the loss via burial of dissolved methane at the base of GHSZ is inhibited (Figure 21).The largest increase in GHI is observed at moderately low burial velocities and mass accumulation rates.Thus, the GHI at w 0 = 50 cm kyr −1 and POC 0 = 1 wt.% amounts to 55.8 g C m −2 at complete compaction (Figure 22) while a value of only 2.8 g C m −2 is calculated for the reference scenario (Figure 16).The burial velocity in surface sediments has only a small effect on gas hydrate accumulation for fully compacting sediments.At POC 0 = 2 wt.% almost the same GHI value is calculated over a w 0 range of 50-150 cm kyr −1 (Figure 22).
where GHI is the inventory of methane hydrate within the GHSZ (in kg methane-C m −2 ), L GHSZ is the thickness of the gas hydrate stability zone (in m) while POC 0 (in wt.%) and w 0 (in cm kyr −1 ) are the POC concentration and the burial velocity of solids in surface sediments.Parameter values of the equation were determined by a minimization procedure using the FindMinimum object in MATHEMATICA 8.0.The best fit to the model runs was obtained with: a = 0.002848 ± 0.00049, b = 1.681 ± 0.027, c = 24.42 ± 7.2, d = 0.9944 ± 0.10, e = −1.441± 0.19, f = 0.3925 ± 0.032.
The results of the 80 model runs are well reproduced by the transfer function (Figure 23).The function is valid for the parameter values applied in the reference simulation (Table 2 with Φ f = 0.0, p = 1/1300) for GHI ≥ 0 over a range of L GHSZ = 0-1000 m; w 0 = 10-150 cm kyr −1 , and POC 0 = 0.2-5 wt.%.Controls on GHI are shown in Figure 24.The contour plots look similar to those produced for the reference simulation (Figure 19).However, significant differences arise at w 0 = 20-100 cm kyr −1 and POC 0 = 1-2 wt.%where the function for full compaction predicts significantly higher gas hydrate inventories in margin sediments.These differences are relevant since the POC concentrations and burial velocities observed in margin sediments often fall into these w 0 and POC 0 ranges.

Effects of Fluid and Gas Flow on Gas Hydrate Accumulation
The accumulation of gas hydrates is strongly supported by the ascent of pore fluids into the GHSZ.Model runs under reference conditions (Table 2) show that the gas hydrate inventory is greatly enhanced already at low upward fluid flow velocities (Figure 25).The inventory increases by two orders of magnitude at an upward Darcy flow of u W = 1.5 cm yr −1 .The gas hydrate saturation at the base of the GHSZ increases to ~90% at u W = 2.2 cm yr −1 .At this point the infiltration of fluids into the GHSZ is effectively suppressed since the pore space is almost completely clocked by methane hydrate.The effect of gas inflow was estimated assuming that the ascending gas is completely dissolved at the base of the GHSZ (Figure 25).The modeling showed that under this assumption an almost complete saturation of the pore space is achieved at an upward gas flow of u G = 0.025 cm −1 yr −1 .Fluid and gas infiltration are thus effectively supporting the accumulation of gas hydrates in marine sediments.The model runs based on the reference scenario produced a maximum gas hydrate saturation of only ~3% without gas and fluid flow.It is thus very likely that the much higher saturations that have been observed at many continental margin sites have been created almost exclusively by the ascent of pore fluids and/or methane gas into the GHSZ.

Effects of POC Degradation Kinetics on Gas Hydrate Accumulation
The initial age (age 0 ) and the inhibition constant (K C ) are used in the model to define the reactivity of POC [Equations (7,11)].Both parameters are varied to study their effect on gas hydrate accumulation applying the parameter settings of the reference simulation (Table 2).The model runs show that GHI increases with age 0 until a maximum is reached at age 0 = 5 kyr (Figure 26).A further increase causes a significant contraction of the GHI until gas hydrate formation is almost completely suppressed at age 0 = 100 kyr.The initial increase is related to a better preservation of POC within the sulfate reduction zone supporting methanogenesis in the underlying sediments.A further increase in age 0 reduces the rate of methanogenesis within the GHSZ and thereby the accumulation of gas hydrates.The dissolved sulfate profile attains an almost linear shape at age 0 = 100 kyr since most of the sulfate is consumed via AOM rather than POC degradation under this parameter setting (Figure 27).Similar sulfate profiles have been observed at many continental margin sites.The modeling strongly suggests that they are related to the deposition of old POC at the sediment surface.The high initial age of this material (100 kyr) may reflect the deposition of POC previously eroded at other locations and subsequently re-deposited at the study site (i.e., drift deposits such as Black Ridge).It may also reflect the deposition of terrigenous POC eroded from adjacent continental areas.The inhibition constant (K C ) has a very strong effect on gas hydrate accumulation.Methanogenesis is strongly suppressed at low K C values while high values allow for significantly more methane production and gas hydrate accumulation (Figure 26).The down-core increase in dissolved ammonium, bromide, iodide, total alkalinity and DIC observed in ODP cores provides important constraints on the POC degradation rates and can thus be used to constrain the ambient K C value.The evaluation of these pore water profiles at major gas hydrate provinces showed that K C usually falls into a range of 30-50 mM [23,33].
The Middelburg model [Equations (7,8)] was integrated into the transport-reaction model to further test the effects of reaction kinetics on gas hydrate accumulation.It allows for more gas hydrate formation than the model used in this paper [23] since methanogenesis is not inhibited by DIC and methane accumulating in ambient pore fluids (Figure 28).The GHI value of 910 kg C m −2 calculated with Middelburg's kinetic rate law is almost identical to the GHI values calculated with Wallmann's kinetic equation at high K C values (GHI = 907 kg C m −2 at K C = 100,000 mM, s. Figure 26).A new kinetic rate law for POC degradation was recently introduced by Burdige [70].Burdige proposed that microbial POC degradation rates increase with sediment temperature and that this increase can be described by a simple Arrhenius equation: where k(T) is a temperature-dependent kinetic constant, E a is the activation energy, T is sediment temperature and R is the gas constant (R = 8.314472 J mol −1 K −1 ).The corresponding non-dimensional constant describing the down-core increase in reactivity is defined as: where T BW is bottom water temperature.We introduced this constant into our rate law to calculate gas hydrate accumulation for the reference case considering the effects of both metabolite inhibition and temperature on POC degradation: (35) Applying an activation energy of E a = 200 kJ mol −1 [70] the value of the kinetic constant k T increases exponentially from 1 at the sediment surface to 197 at the base of the sediment column.The rates of POC degradation are thus strongly enhanced in the lower section of the model column (Figure 28).The gas hydrate inventory is expanded to 520 kg C m −2 due to the temperature effect.Hydrate saturations calculated with the Middelburg model [Equations (7,8)] and the Burdige model [Equation (35)] vary within the same order of magnitude (Figure 28).It should, however, be noted that the strong down-core decrease in POC and increase in DIC produced by the Middelburg and Burdige models are not consistent with observations.The POC profiles measured in ODP cores from major gas hydrate provinces show only a minor decrease with sediment depth while the DIC concentrations do not exceed 300 mM [23,33,71].The rates of POC degradation, methane production and gas hydrate formation are thus strongly over-estimated by the Middelburg and Burdige models.The Middelburg equation is purely empirical while the Arrhenius equation applied by Burdige is valid for elementary reactions but not for complex multi-step process such as the microbial degradation of POC.Both approaches are thus not based on solid theoretically concepts.The chosen apparent activation energy of 200 kJ mol −1 has been determined for abiotic hydrocarbon generation from kerogen while values reported for microbial POC degradation in surface sediments span a large range between 23 and 132 kJ mol −1 with values clustering between 50 and 60 kJ mol −1 and 70 and 90 kJ mol −1 [72,73].
In the following we use OPD Site 799A (Leg 128) featuring a strong geothermal gradient and high-quality dissolved ammonium data to test our model at elevated temperatures of up to 50 °C.The shape of the ammonium depth profile is giving us important constraints on the down-core change in organic matter reactivity since the biogeochemistry of ammonium is very simple in anoxic sediments: Ammonium is released into the porewater by the microbial degradation of organic matter while a minor portion is adsorbed on the solid phase until an adsorption/desorption equilibrium is attained.The stoichiometry of ammonium release is constrained by field data (POC/PON ratios) while the adsorption constant has been defined experimentally [74].OPD Site 799A is located in the Japan Sea, a marginal sea in a continental arc setting [75].Sediments were taken at about 2000 m water depth from the Kita-Yamato Trough, a narrow graben structure filled with sediments featuring a very high geothermal gradient of about 100 °C km −1 .Sediments cored in Hole 799A include biosiliceous sediments and fine-grained detrital sediments intercalated with carbonate-rich intervals and numerous ash layers.They accumulate at an average rate of about 70 m Myr −1 and contain abundant organic matter of mostly marine origin with particulate organic carbon values (POC) ranging from 0.24% to 5.66%.Porosity values scatter around a mean value of 0.75.Both, POC and porosity show no significant down-core trends.Pore water data indicate that sediments retrieved from Hole 799A are subject to in-situ diagenetic processes including microbial degradation of sedimentary organic matter.Pore-water gradients and concentrations are not significantly affected by upward fluid flow [75].Our transport-reaction model [23] was applied to site 799A to simulate the microbial degradation of POC and particulate organic nitrogen (PON) and the turnover of ammonium.The reaction network of the model includes degradation of particulate organic matter via sulfate reduction and methanogenesis, anaerobic oxidation of methane, the release of dissolved ammonium from degrading PON, and ammonium adsorption.The effects of porosity, tortuosity and temperature on diffusion coefficients of dissolved sulfate, methane, ammonium and dissolved inorganic carbon (DIC) were considered in the model applying appropriate empirical equations [76].Constant concentrations were prescribed at the upper boundary of the model column while zero gradients were used as lower boundary condition.The following parameter values were applied in the modeling: porosity: 0.75 (constant), molar PON/POC ratio: 0.1, burial velocity: 70 m Myr −1 , geothermal gradient: 100 °C km −1 , K C : 40 mM, age 0 : 10 kyr, adsorption constant for ammonium: 1.7 cm 3 g −1 [74].
The model results clearly show that the kinetic model for POC degradation introduced by Wallmann et al. [23] is able to reproduce the data measured at OPD site 799A (s. Figure 29).The ammonium profiles calculated by the model are very close to observations and the model results are consistent with solid phase data (POC and total N).ODP site 799A was also simulated applying the original Middelburg model without inhibition term.The model results shown in Figure 29 confirms that the Middelburg model overestimates rates of organic matter degradation at large sediment depths since the predicted ammonium concentrations clearly exceed the measured values.This observation is consistent with a recent study on microbial sulfate reduction rates in fjord sediments showing a more pronounced down-core decrease in POC reactivity than predicted by the Middelburg model [77].The data were used to further test the Burdige model for a moderate activation energy of only E a = 50 kJ mol −1 .Even with this modest temperature sensitivity the model strongly over-predicts ammonium release and POC degradation rates at large sediment depths (s. Figure 29).
Microbial organic matter degradation in the deep biosphere is a very complex and poorly understood process.Nevertheless, our model is clearly giving a good approximation to the general trends observed in geochemical ODP data and can thus be used to predict organic matter degradation and microbial methane formation in anoxic marine sediments to sediment temperatures of up to 50 °C.The good correspondence between observations and model results is very encouraging but also surprising since our model does not explicitly consider effects of temperature on microbial degradation rates even though the microbial degradation of labile organic matter in marine surface sediments is promoted under high temperatures [72,73].A weak temperature effect is implicitly considered in our model: Rates are suppressed by the built-up of high metabolite concentrations in ambient pore fluids which is partly controlled by the diffusive flux of these species towards the sediment surface.Since diffusion coefficients increase with temperature, elevated temperatures promote the diffusive loss of metabolites and thereby the rate of organic matter degradation.Thus, the rate of POC degradation amounts to 5.9 µmol m −3 yr −1 at the base of the model column (450 mbsf) in the model run depicted in Figure 29 and decreases to 5.5 µmol m −3 yr −1 when the sediment temperature is kept at a constant value of 2 °C throughout the model column.Considering this very modest temperature sensitivity, our model should significantly underestimate degradation rates at about 50 °C if microbial turnover rates in the deep biosphere would strongly respond to temperature.It thus seems that temperature is not having a strong effect on the kinetics of refractory organic matter degradation in the deep biosphere.Rates are rather controlled by the strong down-core decrease in reactivity induced by the preferential degradation of labile organic matter at shallow depth levels and the inhibition of degradation processes by the accumulation of dissolved metabolites.Considering the long time periods available for adaptation in slowly accumulating sediments, it might be possible that each sediment level is harboring microorganisms being adapted to ambient temperature.This would suggest a down-core succession of microbial communities along the geothermal gradient with each community operating close to its temperature optimum.In this scenario, ambient temperature would have a negligible effect on the microbial degradation of refractory organic matter in the deep biosphere.However, much more basic research on microbial processes, diversity, and adaptation in the deep biosphere is needed to test this hypothesis.

Constraining the Global Methane Hydrate Inventory in the Modern Ocean
In the following section, the transfer functions derived in Sections 7.2 and 7.3 are applied on a global grid to estimate the abundance and distribution of methane hydrate in marine sediments.The global grid was defined by the POC concentrations in marine surface sediments (POC 0 ) and burial velocities (w 0 ) reported in Burwicz et al. [6] and the GHSZ map (L GHSZ ) shown in Figure 9.
Transfer function Equation (31) which is valid for sediments experiencing a low degree of compaction was applied with Holocene burial velocities to provide a minimum estimate for the accumulation of gas hydrates in marine sediments.The global inventory of methane hydrates in marine sediments amounts to only 3 Gt of C under these boundary conditions.Low-grade gas hydrate occurrences are predicted for some high productivity areas with depth-integrated hydrate inventories of less than 40 kg m −2 .These low values are consistent with the previous model results for Holocene boundary conditions presented by Burwicz et al. [6] yielding a global gas hydrate inventory of only 4 Gt of C.
The transfer function valid for low compaction [Equation ( 31)] was also applied on a grid with Quaternary burial velocities to consider the down-slope transport of sediments during glacial sea-level low-stands (Figure 30).The global inventory of methane hydrates in marine sediments expands to 122 Gt of C due to the increase in burial velocity.The map shows extended hydrate-bearing areas at productive continental margins.Elevated inventories of up to 230 kg m −2 are found at the highly productive margin off Peru.Our best estimate of the global distribution of marine sediments is shown in Figure 31.It is based on the assumptions that margin sediments are completely compacted at large sediment depths and that shelf sediments are eroded and deposited on the continental slope and rise under glacial conditions.The global inventory of gas hydrates derived from our best estimate map (Figure 31) is 455 Gt of C. The map predicts gas hydrate accumulations of ~10 kg m −2 over extended margin areas with elevated values of >100 kg m −2 in high productivity zones at the western continental margins and in the Arabian Sea.High values are also calculated at the Argentinean margin and north of Madagascar.The background concentration of methane hydrates in margin sediments of ~10 kg m −2 corresponds to very low gas hydrate saturations of less than 0.1 Vol.% which are not resolved by geophysical methods and drilling.The true distribution of gas hydrates is probably patchier than the distribution shown in Figure 31 since sediments are not always experiencing full compaction such that dissolved methane is buried and lost below the GHSZ.Fluid flow is implicitly considered in the map since pore fluids are not buried below the GHSZ.This approach is probably valid for passive continental margins.However, it does not consider the additional fluid flow at active margins where pore fluids of the incoming oceanic plate are mobilized and released into the GHSZ.The gas hydrate abundance at active margins is thus systematically underestimated in Figure 31.It should also be noted that steady-state conditions were applied in all model runs.The models were run over a period of several million years until the depth-integrated inventory of gas hydrates reached a constant value.It was assumed that geothermal gradients and POC accumulation rates at the seafloor have been constant over the entire model period.Thus, the model runs do not consider subsidence, thermal evolution and the changes in POC accumulation during basin formation.
Gas ascent and focused fluid flow at active continental margins are not considered in the map.These processes are responsible for the formation of high-grade gas hydrate accumulations which have been documented at an ever increasing number of sites [78,79].The global picture that emerges from our modeling exercise may be summarized as: Gas hydrates occur at many continental margin sites at very low saturations of less than 0.1%.Higher saturations of 1%-3% appear at productive continental margins where abundant organic matter is delivered to the seabed by the sinking of biogenic detritus.High-grade deposits with saturations of >3% are formed only where fluid and/or gas have migrated into the GHSZ.These deposits are not shown on the map since it is currently not possible to constrain the global distribution of these high-grade gas hydrate deposits.
Our estimate of the global abundance of marine gas hydrates is certainly too low since gas ascent and focused fluid flow are not considered.Fluid and gas flow have an enormously strong effect on gas hydrate accumulation.The gas hydrate inventory is enhanced by about two orders magnitude already at very low velocities of fluid and gas ascent (u w < 2 cm kyr −1 ; u G < 0.02 cm kyr −1 , s. Figure 25).It is currently not possible to constrain the global number and abundance of these hot spots for gas hydrate accumulation.However, the modeling presented in this contribution suggests that the global inventory of methane in these hot spot areas may exceed the global hydrate inventory of low-grade deposits if only 1% of the global seafloor would be affected by upward gas or fluid flow.
The gas hydrate distribution shown in Figure 31 may be used for prospection purposes while Equation ( 32) can be applied to estimate the gas hydrate inventory where ambient data for w 0 , L GHSZ and POC 0 are available.It should, however, be noted that detailed local surveys will almost certainly reveal inventories differing considerably from the predicted values since the regional basin evolution and the local pathways for fluid and gas flow are not resolved in our global estimates.Seismic surveys and drilling data are needed to constrain local gas hydrate inventories.Moreover, new 3-D basin modeling tools have recently been developed that can be used to better constrain gas hydrate inventories in well characterized sedimentary basins [80].These models can be used to simulate in detail the thermal and sedimentary history and the ascent of gas and fluids through faults and fractures at selected sites.However, it will not be possible to apply these basin models at global scale in the foreseeable future.
The first model-based estimate of the global methane hydrate inventory under Holocene boundary conditions was presented by Buffett and Archer [3].They used the POC rain rate to the seabed as a major external driving force for the simulation of hydrate formation in marine sediments.The rain rate was calculated as a function of water depth and applied as an upper boundary condition for an early diagenesis model simulating the degradation of POC in the top meter of the sediment column [81].The POC burial rate at 1 m sediment depth calculated with this "muds" model was applied as upper boundary condition for the simulation of methane turnover in the underlying sediment sequence.The later simulations consider the microbial degradation of POC via sulfate reduction and methanogenesis and the anaerobic oxidation of methane within the sulfate-methane transition zone [82,83].POC was separated into an inert fraction and a labile fraction that was degraded over the top km of the sediment column following first order reaction kinetics (rate constant 3 × 10 −13 s −1 ).The model also considers an externally imposed rate of upward fluid flow.It was calibrated using field data obtained at Blake Ridge (a passive margin site) and the Cascadia margin (active continental margin).The best fit to observations was obtained assuming that the 25% of the total POC is labile and using upward fluid flow velocities of 0.23 mm yr −1 (passive margin) to 0.4 mm yr −1 (active margin).Applying these parameter values on a global grid and assuming a compensating downward fluid flow over 50% of the global seafloor area resulted in a total methane hydrate inventory of 3000 Gt C. The hydrate inventory was to a large degree controlled by the velocity of upward fluid flow.Without imposed fluid flow, the global hydrate inventory was reduced to 600 Gt C. Subsequently, the authors discovered an extrapolation error in the calculation of POC rain rates as function of water depth [5].The model-based estimate of the global hydrate inventory was reduced from 3000 Gt C to approximately 700-900 Gt C after correction for this error [5].Klauda and Sandler [4] used a modified version of the Davie and Buffett model [82,83] to estimate the global marine hydrate inventory under Holocene boundary conditions.The entire POC pool was assumed to be completely degradable with a reduced decay constant of 1.5 × 10 −14 s −1 .The upper boundary of the model domain was constrained by field data.POC concentrations measured in surface sediments and Holocene sedimentation rates averaged over the major ocean basins were applied to define the POC burial flux at the sediment surface.The global hydrate inventory estimated with this approach was surprisingly high (55-700 Gt C) even though the model was run without imposing upward fluid flow.Unfortunately, Klauda and Sandler [4] did not consider POC degradation via sulfate reduction.The large difference to our estimate may be explained by the neglect of microbial sulfate reduction and AOM since the accumulation of methane in pore fluids and gas hydrates is strongly diminished by these processes.
Burwicz et al. [6] performed the first simulation of gas hydrate formation under Quaternary boundary conditions at global scale.Their estimate of the global gas hydrate inventory in marine sediments (995 Gt C) is larger than our Quaternary estimates for low and complete compaction (122 Gt C and 455 Gt C, respectively).Burwicz et al. [6] applied a different compaction law and assumed very low porosities for GHSZ sediments.Moreover, they used a significantly higher value for the mean thermal conductivity of marine sediments (1.5 W m −1 K −1 ) resulting in a very strong expansion of the GHSZ with respect to our modeling (s.discussion in section 4).To test the effect of the GHSZ expansion, we applied their GHSZ grid in an additional model run.We obtained a global GHI value of 937 Gt C for Quaternary conditions and full compaction [Equation (32)] which is close to the value derived by Burwicz et al. [6] (995 Gt C).We conclude that the difference between our best estimate (455 Gt C) and the previous estimate (995 Gt C) by Burwicz et al. [6] is mainly caused by the different choice of thermal conductivity values.The maps presented by Burwicz et al. [6] were produced by performing a full-scale 1-D transport-reaction simulation at each grid point.They are thus more precise than our estimates which are affected by additional errors introduced by the transfer functions.However, the available data strongly suggest that the thermal conductivity of most continental margin sediments is smaller than 1.5 W m −1 K −1 and close to the value applied in our modeling (1.0 W m −1 K −1 ).We thus conclude that the gas hydrate distribution shown in Figure 31 and our new best estimate of the global gas hydrate abundance (455 Gt C) are more realistic than previous results.

Accumulation of Gas Hydrates under PETM Boundary Conditions
The strong negative isotopic carbon excursion observed during the PETM (55.5 Ma ago) has been ascribed to the dissociation of methane hydrates in marine sediments [18,84,85].About 1000 Gt of methane carbon release are required to explain the isotope values observed in marine carbonates.The conditions at the PETM were very different from today.The bottom water temperatures were close to 10 °C, the burial velocities of terrigenous sediments were reduced by a factor of ~4-5 while the POC concentrations in surface sediments were probably much higher than today (Figure 4).
To test the potential for gas hydrate formation prior to onset of the PETM, the standard simulation with full compaction was run at a bottom water temperature of 10 °C, a POC concentration of 4 wt.% and a burial velocity of w 0 = 33.7 cm kyr −1 (Figure 32) The thickness of the GHSZ was reduced to ~265 m due to the elevated bottom water temperature while methane production rates were enhanced due to the increase in POC concentration.The resulting gas hydrate saturation values were similar to the saturations obtained under modern conditions (Figure 32).However, the overall effect was a decrease in the GHI from 99 kg C m −2 in the modern ocean to 19 kg C m −2 in the PETM.These model results suggest that the potential for gas hydrate accumulation before the onset of the PETM was lower than today.Wide-spread ascent of pore fluids and/or methane gas into the GHSZ has to be postulated for the period prior to the PETM to explain the isotopic excursion by gas hydrate dissociation.

Conclusions and Open Questions
A number of improved global estimates can be derived from the data and model results presented above:  The global POC accumulation rate in the modern ocean amounts to ~1.4 × 10 14 g yr −1 . The global rate of microbial methane formation in marine sediments is lower than previously assumed and falls into the range of 4 − 25 × 10 12 g C yr −1 . The pore volume within the GHSZ is extensive and can be constrained as ≥44 × 10 15 m 3 .
 The global inventory of methane hydrates in marine sediments is ≥455 Gt C.
The global distribution maps of methane hydrate presented in this contribution are representative for methane hydrate being formed by the microbial degradation of POC within the GHSZ.The maximum pore space saturation reached under these conditions is 3 vol.%.However, sediment cores and logs obtained at an ever growing number of sites have clearly shown that the saturation values are much higher than 3 vol.%under favorable geologic conditions.Our model results clearly show that these highly concentrated gas hydrate deposits have been formed by the migration of dissolved or gaseous methane into the GHSZ.Potentially exploitable gas hydrate reservoirs are thus only formed where gas or water has been transported into the GHSZ via upward advection.Fluid and gas flow have an enormously strong effect on gas hydrate accumulation.The gas hydrate inventory is enhanced by about two orders magnitude already at very low velocities of fluid and gas ascent (u w < 2 cm kyr −1 ; u G < 0.02 cm kyr −1 ).
The global distribution of methane hydrates is probably characterized by the wide-spread occurrence of low-grade hydrate saturations (<0.1 vol.%) along continental margins superimposed by elevated saturations in highly productive regions (1-3 vol.%).Observational data reveal a large number of "sweet spots" where saturations are enhanced by methane migration (>3 vol.%).However, it is currently not possible to constrain the global number and abundance of these high-grade deposits by numerical modeling.
The Quaternary is characterized by very low bottom water temperatures and extremely high rates of POC accumulation at the seabed.It offers unusually favorable conditions for gas hydrate accumulation since the inventory of gas hydrates in marine sediments increases under low temperature and high POC input.The conditions were much less favorable during most of the Cenozoic and Cretaceous.This is also true for the PETM where large-scale gas hydrate dissociation has been evoked to explain a global negative excursion observed in the δ 13 C record.Thus, wide-spread ascent of pore fluids and/or methane gas into the GHSZ has to be postulated for the period prior to the PETM to explain the isotopic excursion by large-scale gas hydrate dissociation.
More work is needed to better constrain the global gas hydrate inventory in the seabed.This work should address at least some of the following open questions:  What are the major geochemical and biological control mechanisms regulating the rate of microbial methane production in the deep biosphere? What precisely are the physical and geological controls on gas transport into the GHSZ?Is gas ascending into the GHSZ already at low gas saturation values?What defines the critical gas saturation value that needs to be overcome to allow for gas migration? What fraction of the global seabed is affected by upward gas and fluid flow? How does the hydrate inventory evolve under non-steady state conditions?and his team for their excellent editorial handling of our manuscript.The thoughtful comments by the two anonymous reviewers greatly helped to improve our manuscript.

Figure 1 .
Figure 1.Concentration of particulate organic carbon (POC) in Holocene surface sediments [7,8].The color coding indicates POC values in wt.%.Regions of the seafloor where data are not available are shown in white.

Figure 4 .
Figure 4. Accumulation of POC and terrigenous sediments over the last 150 Myr.The solid line shows the global marine POC accumulation rate normalized to the Quaternary value while the dotted line gives the global rate of terrigenous sedimentation at the seafloor normalized to the Quaternary [17].The Paleocene-Eocene Thermal Maximum (PETM, at 55.5 Myr) is indicated by the vertical arrow.

Figure 5 .
Figure 5. Fraction of buried POC entering the methanogenic zone (f MI ) as calculated from Equation (9) for initial ages of 1000 yrs, 10,000 yrs, and 100,000 yrs.

Figure 6 .
Figure 6.Fraction of POC converted into methane carbon within the methanogenic zone of continental margin sediments (f M ).The marked area indicates the range of POC residence times (t = 3-30 Myr) within the methanogenic zone at continental margins.

Figure 7 .
Figure 7. Phase diagram for methane hydrate (structure type I) for sulfate-free seawater with a salinity of 35 [34].Arrows indicate the formation of a transition zone by capillary forces and the broadening of the gas hydrate stability zone by higher hydrocarbon compounds (C 2+ ).

Figure 9 .
Figure 9. Thickness of the gas hydrate stability zone calculated from global heat flux data [43] considering bathymetry, bottom water temperature, sediment thickness data and a bulk thermal conductivity of 1.0 W m −1 K −1 [50].

Figure 10 .
Figure 10.Solubility of methane in marine sediments deposited at 2000 m water depth.Solubility is calculated for methane hydrate (blue circles) and methane gas (red crosses) in sulfate-free seawater with S = 35 assuming hydrostatic conditions, a bottom water temperature of 2 °C and a linear geothermal gradient of 0.03 °C m −1[34,51].The black line indicates the dissolved methane concentration attained in methane-saturated pore fluids (in mmol per kg solution).The base of the gas hydrate stability zone (GHSZ) is located at the intersection of the gas hydrate and free gas solubility curves.

Figure 12 .
Figure 12.Velocity of Darcy fluid flow at lithostatic pressure (u L ) in marine sediments as function of intrinsic permeability.Darcy flow is calculated applying density values of ρ W = 1036 kg m −3 and ρ S = 2500 kg m −3 and a porosity of Φ = 0.5.

Figure 13 .
Figure 13.Relative permeability of gas [k RG , upper panel, Equation (26)] and buoyancydriven gas flow (u G-H , lower panel, Equation 30) in fine-grained marine sediments as function of gas saturation (S G = 1 − S W ). k RG is calculated applying Equations 26 and 27 with S WC = 0.3, S GC = 0.01 and m = 0.75.The gas flux is calculated for hydrostatic pressure conditions applying Equation 30 with ρ W = 1036 kg m −3 , ρ G = 198 kg m −3 , k = 10 −15 m 2 , µ G = 11 × 10 −6 Pa s and the relative gas permeability depicted in the upper panel.

Figure 14 .
Figure 14.Concentrations of dissolved inorganic carbon (DIC), dissolved methane (CH 4 ), particulate organic carbon (POC) and methane hydrate saturation (GH-Sat.) in the reference model run.The dotted line in the upper right panel indicated the concentration of methane in equilibrium with methane hydrate.

Figure 15 .
Figure 15.Concentrations of dissolved sulfate (solid line) and dissolved methane (dotted line) in surface sediments calculated in the reference model run.The shaded bar indicates the sediment layer where methane is oxidized by microbial consortia using sulfate as terminal electron acceptor (AOM zone).The AOM zone is equivalent to the sulfate-methane transition zone (SMTZ).

Figure 16 .
Figure16.Effect of POC concentration and burial velocity on methane hydrate accumulation in marine sediments.POC 0 and w 0 are the concentration of particulate organic carbon and the burial velocity of solids in surface sediments while GHI is the depth-integrated inventory of carbon bound in methane hydrate.R indicates the reference case.

Figure 17 .
Figure 17.Effect of GHSZ thickness on methane hydrate accumulation in marine sediments.L GHSZ denotes the thickness of the gas hydrate stability zone.It was varied by changing the heat flow applied in the reference model run over a range of q H = 0.02-0.1 W m −2 corresponding to temperature gradients of 20-100 °C km −1 .In the left panel, the burial velocity at the interface was maintained at the reference value (w 0 = 134.6 cm kyr −1 ).The right panel shows model results with diminished w 0 values.R indicates the reference case.

Figure 18 .
Figure 18.Gas hydrate inventories (GHI in kg m −2 ) within the GHSZ calculated by individual runs of the transport-reaction model (GHI M , s. Figures 16,17) and with the transfer function [GHI T , Equation (31)] for the same set of parameter values.The solid line shown in the left panel has a slope of unity and serves as a reference for the transfer function values.The residuals shown in the right panel are calculated as GHI T -GHI M .
Figure 18.Gas hydrate inventories (GHI in kg m −2 ) within the GHSZ calculated by individual runs of the transport-reaction model (GHI M , s. Figures 16,17) and with the transfer function [GHI T , Equation (31)] for the same set of parameter values.The solid line shown in the left panel has a slope of unity and serves as a reference for the transfer function values.The residuals shown in the right panel are calculated as GHI T -GHI M .

Figure 19 .
Figure 19.Controls on gas hydrate accumulation as predicted by the transfer function [Equation (31)].The contour lines indicate gas hydrate inventories (GHI) in kg C m −2 .The contour plot on the left hand side is calculated for w 0 = 100 cm kyr −1 while the contour plot on the right shows the results for L GHSZ = 400 m.

Figure 20 .
Figure 20.Effect of sediment compaction on gas hydrate accumulation.The porosity after compaction (Ω f ) and the attenuation coefficient for the down-core decrease in porosity (p) are varied while the porosity at the sediment surface is maintained at the reference value of Ω 0 = 0.74 [s.Equation (15)].Compaction decreases towards the right in both plots.R indicates the reference case.

Figure 22 .
Figure 22.Effect of POC concentration and burial velocity on methane hydrate accumulation in continental margin sediments experiencing complete compaction (Φ f = 0.0, p = 1/1300 m −1 ).See legend of Figure 16 for further information.

Figure 23 .
Figure 23.Gas hydrate inventories (GHI in kg m −2 ) within the GHSZ calculated by individual runs of the transport-reaction model (GHI M ) and with the transfer function [GHI T , Equation (32)] for fully compacting margin sediments (see legend of Figure 18 for further information).

Figure 24 .
Figure 24.Controls on gas hydrate accumulation as predicted by the transfer function for margin sediments subject to complete compaction.The contour lines indicate gas hydrate inventories (GHI) in kg C m −2 .The contour plot on the left hand side is calculated for w 0 = 100 cm kyr −1 while the contour plot on the right shows the results for L GHSZ = 400 m.

Figure 25 .
Figure 25.Effect of fluid flow (u W ) and gas flow (u G ) on gas hydrate accumulation.R indicates the reference case.

Figure 26 .
Figure 26.Kinetic controls on gas hydrate accumulation.age 0 is the initial age of POC while K C is a constant defining the inhibition of POC degradation induced by the accumulation of dissolved metabolites (DIC, CH 4 ) in ambient pore fluids.R indicates the reference case.

Figure 27 .
Figure 27.Dissolved sulfate and methane profiles at age 0 = 100 kyr.See legend of Figure 15 for further information.

Figure 29 .
Figure 29.Model results for ODP site 799A.Concentrations of dissolved ammonium (NH 4 ), particulate nitrogen (N), particulate organic carbon (POC) and POC degradation rates obtained applying the kinetic rate laws introduced by Wallmann et al. [Equation (11)], Middelburg [Equations (7,8)] and Burdige [Equation (35)] with an activation energy of E a = 50 kJ mol −1 .N concentrations are the sum of particulate organic nitrogen (PON) and adsorbed ammonium.Symbols indicate the data measured at ODP site 799A while lines indicate model results.

Figure 30 .
Figure 30.Global distribution of methane hydrate in marine sediments.The map was calculated applying low compaction [Equation (31)] and Quaternary burial velocities.The color coding indicates depth-integrated hydrate inventories in kg C m −2 .

Figure 31 .
Figure 31.Global distribution of methane hydrate in marine sediments.The map was calculated applying complete compaction [Equation (32)] and Quaternary burial velocities.The color coding indicates depth-integrated hydrate inventories in kg C m −2 .

Figure 32 .
Figure32.Concentrations of dissolved inorganic carbon (DIC), dissolved methane (CH 4 ), particulate organic carbon (POC) and methane hydrate saturation (GH-Sat.)calculated for the modern ocean and the PETM.The modern ocean simulation was done for the reference case with complete compaction (Φ f = 0.0, p = 1/1300 m −1 ).The PETM simulation was based on the reference case with complete compaction applying a bottom water temperature of 10 °C, a POC concentration of POC 0 = 4 wt.% and a reduced burial velocity of w 0 = 33.7 cm kyr −1 .

Table 1 .
Quaternary POC and methane balance for global continental margins.

Table 2 .
Parameter values and upper boundary conditions applied in the reference case.

Table 3 .
Fluxes at the upper boundary of the model column, depth-integrated rates and carbon inventories calculated in the reference simulation.