References

The accumulation of gas hydrates in marine sediments is essentially controlled by the accumulation of particulate organic carbon (POC) which is microbially converted into methane, the thickness of the gas hydrate stability zone (GHSZ) where methane can be trapped, the sedimentation rate (SR) that controls the time that POC and the generated methane stays within the GHSZ, and the delivery of methane from deep-seated sediments by ascending pore fluids and gas into the GHSZ. Recently, Wallmann et al. (2012) presented transfer functions to predict the gas hydrate inventory in diffusion-controlled geological systems based on SR, POC and GHSZ thickness for two different scenarios: normal and full compacting sediments. We apply these functions to global data sets of bathymetry, heat flow, seafloor temperature, POC input and SR, estimating a global mass of carbon stored in marine methane hydrates from 3 to 455 Gt of carbon (GtC) depending on the sedimentation and compaction conditions. The global sediment volume of the GHSZ in continental margins is estimated to be 60–67 × 1015 m3, with a total of 7 × 1015 m3 of pore volume (available for GH accumulation). However, seepage of methane-rich fluids is known to have a pronounced effect on gas hydrate accumulation. Therefore, we carried out a set of systematic model runs with the transport-reaction code in order to derive an extended transfer function explicitly considering upward fluid advection. Using averaged fluid velocities for active margins, which were derived from mass balance considerations, this extended transfer function predicts the enhanced gas hydrate accumulation along the continental margins worldwide. Different scenarios were investigated resulting in a global mass of sub-seafloor gas hydrates of ~ 550 GtC. Overall, our systematic approach allows to clearly and quantitatively distinguish between the effect of biogenic methane generation from POC and fluid advection on the accumulation of gas hydrate, and hence, provides a simple prognostic tool for the estimation of large-scale and global gas hydrate inventories in marine sediments.


Introduction
Submarine gas hydrates (GH) have been recovered in more than 40 regions worldwide and their presence has been deduced from geophysical, geochemical and geological evidences at more than 100 continental margin sites (e.g. Mazurenko and Soloviev, 2003;Milkov and Sassen, 2002;Lorenson and Kvenvolden, 2007). They have been recognized as a key factor with respect to past and future climate change (e.g. Hester and Brewer, 2009;Dickens, 2001a;Adam, 2002;Archer, 2007;Lunt et al., 2011;Biastoch et al., 2011) and are considered as a new unconventional resource of natural gas (e.g. Sloan, 2003;Boswell, 2009). Models and observations clearly show that GH formation is essentially controlled by (i) the accumulation of particulate organic carbon at the seafloor (POC), (ii) the sedimentation rate, (iii) the kinetics of microbial degradation of organic matter and its related generation of methane, (iv) the thickness of the gas hydrate stability zone (L GHSZ ), (v) the solubility of methane in pore fluids within the gas hydrate stability zone (GHSZ), and (vi) the ascent of methane-rich pore fluids and gas from deeper sediment strata into the GHSZ.
In contrast to the general progress in gas hydrate research made so far, the global abundance of GH in marine sediments still remains poorly constrained. Estimates are generally based on extrapolation of field data (e.g. Kvenvolden and Claypool, 1988;Dickens, 2001b;Milkov, 2004) and geochemical transport-reaction modelling (e.g. Buffett and Archer, 2004;Klauda and Sandler, 2005;Archer et al., 2009;. Although some of the earliest estimates were based on unrealistic conditions (e.g. Dobrynin et al., 1981) and can likely be excluded, they still range over three orders of magnitude (1 × 10 2 -5.6 × 10 4 GtC; cf. Fig. 1) and a clearly constrained consensus value has not emerged over the past decades. This uncertainty implies a major drawback since the resource potential and the possible impact of methane hydrates on past and future climate change cannot be evaluated without a well constrained estimate of global GH abundance and its distribution. The availability of global data sets and the performance of global models greatly improved during the past decade. Grid-based transport-reaction modelling is thus a promising tool for the evaluation of global gas hydrate occurrences. Total estimates of GH based on globally gridded data were first published by Gornitz and Fung (1994) who used available information on bathymetry, ocean bottom temperature, organic carbon concentrations as well as averaged values of heat flow and thermal conductivity in marine sediments. The simplified assumption of a GH saturation of 5-10 vol. % of pore space, was somehow balanced by a relatively small potential volume for GH, however, it lead to a prediction of a huge global inventory of marine gas hydrates ranging from 14 200 to 74 700 GtC. Klauda and Sandler (2005) deduced the same range of GH saturations from systematic model runs with a modified version of the transport-reaction code of Davie and Buffett (2001) and predicted a mass of 55 800 GtC. However, Klauda and Sandler's estimate is unrealistic, because they assumed the entire POC pool to be degradable at a constant reactivity (1.5 × 10 −14 s −1 ) -this is in contradiction to the general knowledge of organic matter degradation slowing down with time/age (e.g. Middelburg, 1989;Westrich and Berner, 1984;Wang and Van Cappellen, 1996;Wallmann et al., 2006). Moreover, important microbial processes such as the degradation of organic matter via sulfate reduction and anaerobic methane oxidation (AOM) are not considered in their model. The studies of Buffet and Archer (2004) and, more recently, Archer et al. (2009) consider sulfate reduction and AOM as defined in the original model of Davie and Buffett (2001). For their best fit they assumed that 25 % of the total organic carbon is available for degradation, as well as distinct fluid flow velocities at passive and active margins, which were defined as a function of the sedimentation rate and compensated by downward fluid flow over 50 % of the global seafloor area. The total GH inventory of 3000 GtC (Buffet and Archer, 2004) was however affected by an extrapolation error and subsequently corrected downward to about 1000 GtC (Archer et al., 2009). However, the model by Buffett and Archer (2004), Archer (2007) and Archer et al. (2009) Wallmann et al. (2006) for the formation of biogenic methane. Their budget ranges from 4 to 995 GtC and covers sedimentation rate values from the Holocene to enhanced accumulations during glacial times, but does not explicitly consider the effect of fluid or gas ascent.
The goal of the present study is to provide a systematic analysis that clearly distinguishes between the effects of POC, SR and fluid advection and to offer a simple, but realistic tool for the prediction of GH inventories. In a first step we apply a transfer function presented by Wallmann et al. (2012), which considers GH formation from biogenic methane in diffusion-controlled systems, on gridded global data sets of the control parameters. Subsequently, an extended transfer function is developed, considering upward fluid advection and consequently providing a more realistic global budget of gas hydrates.

Data sources and methods
A transfer function for the quantification of methane hydrates in marine sediments was recently developed by Marquardt et al. (2010) based on the diagenetic transport-reaction model formulated by Wallmann et al. (2006). Marquardt et al. (2010) proposed that the accumulation rate of particulate organic carbon (POCar), a parameter combining the POC content and the sedimentation rate (SR), and the L GHSZ , are the two most important and independent parameters controlling the formation of GH from biogenic methane formed within the GHSZ. However, the POC content and the sedimentation rate can have opposite effects for the methane production. As example, two geological settings with the same POCar but a different combination of POC and SR will produce a significantly different GH content. A low GH content would result for low POC and high SR, while higher GH contents would result from a higher POC and lower SR, due to the increasing organic matter input and residence time within the GHSZ (e.g. Bhatnagar et al., 2007;Wallmann et al., 2012). Subsequently, Wallmann et al. (2012) presented a set of new functions considering the three independent parameters (POC, SR and L GHSZ ), for two different scenarios: normal compaction, and full compaction. The compaction controls the GH formation in two different ways. Higher compaction will produce more hydrates due to the decrease in porosity and the subsequent increase in CH 4 concentration within the pore water. A higher compaction will also decelerate the loss of methane-rich pore fluids below the GHSZ, due to the increasing difference in burial velocity of pore water and bulk sediment. In general, the new fits by Wallmann et al. (2012) produce more realistic results for low POC and SR, and thin L GHSZ than the previous transfer function by Marquardt et al. (2010). Moreover, Wallmann et al. (2012) consider a full compaction scenario, which may be especially convenient for thick sedimentary deposits in compressive tectonic settings.
For normal compaction, the following equation was defined (Wallmann et al., 2012) where GHI is the inventory of methane hydrate within the GHSZ in kg m −2 of C, L GHSZ is the thickness of the GHSZ in m, POC is the particulate organic carbon in wt. % and SR is the burial velocity of solids in surface sediments in cm kyr −1 , and where a 1 = 0.00317(42); a 2 = 1.862(20); a 3 = 54.8(6.9); a 4 = 1.092 (40); a 5 = 1.183(22); a 6 = 0.0154(43). The errors of the fit parameters are given as 1-sigma standard deviations, and noted according to S.I. convention. The linear correlation matrix of fit parameters can be found in Table A1. For full compaction, the following equation (Wallmann et al., 2012) was applied where b 1 = 0.00285(49); b 2 = 1.681(27); b 3 = 24.4(7.2); b 4 = 0.99(10); b 5 = −1.44(19); b 6 = 0.393(32). The linear correlation matrix of fit parameters can be found in Table A2. The approaches and global data sets used to describe the three main parameters controlling gas hydrate formation are described in the forthcoming section. The transfer functions were applied at each grid-point of the ocean in a spatial resolution of 1 • × 1 • . The total methane inventory was integrated by calculating the surface area for each cell of the grid, using the 1984 World Geodetic System ellipsoid as approximation of the shape of the Earth.

Sedimentation rate
The sedimentation rate generally decreases with increasing distance from the coastline and increasing water depth (e.g. . However, the lack of global data sets on sedimentation rate over the Holocene is the main limitation for any global calculation. Moreover, model simulations show that timescales of some Ma are required to run a typical gas hydrate system into steady-state, so that it is also important to consider the long-term sedimentary history (e.g. Dickens, 2011;Wallmann et al., 2012). In order to follow these considerations, we have applied three different approaches to calculate the SR: Approach #1 derives data grids for recent sedimentation rates; Approach #2 aims to consider the sedimentation history by using crustal ages as time constraint; and Approach #3 estimates average sedimentation rates over the Quaternary.

Approach #1
Burwicz et al. (2011) recently published a function for the calculation of Holocene sedimentation rates as a function of water depth: where w is the burial velocity in cm yr −1 , and z is the water depth in m. The best fit to the measured data was obtained for w 1 = 0.117; w 2 = 0.006; z 1 = 200; z 2 = 4000; y 1 = 3 and y 2 = 10. The distribution of burial velocities was calculated by applying Eq. (1) to a 1 • × 1 • global data set of bathymetry (Amante and Eakins, 2009).

Approach #2
Because Approach #1 only provides estimates for modern SR values, which may not hold true for the longer geological timescales on which GH are formed, we provide a comparative data set where we calculated a long-term averaged burial velocity by combining data of the total thickness of marine sediments (Divins, 2010) with the age of the oceanic crust (NOAA, 2010) (Eq. 4).
In order to avoid inconsistencies in young ages, and to prevent infinite values in further calculations, ages younger than 1Myr were discarded. These ages only correspond to the areas adjacent to the mid-oceanic ridges where GH are not formed anyhow because heat flow is too high. Besides these 2 calculations, an additional data set of average SR during the Quaternary was kindly provided by E. Burwicz. In this grid, the main deposition areas have been shifted from the continental shelves into the deeper continental slopes and rises, as a consequence of the low-stand sea level and the enhanced down-slope transport of sediment during glacial conditions typical of the Quaternary (e.g. Hay, 1994). The resulting maps are shown in Fig. 2.

Particulate organic carbon
The POC that accumulates in the sediment on longer timescales depends on the total organic carbon flux to the seafloor (rain rate), and the remineralization processes within the top centimetres of the sediment column (e.g. Suess, 1980;Martens et al., 1992;Seiter et al., 2005;Romankevich et al., 2009;Zonneveld et al., 2010). In order to find good approximations to the distribution of the POC incorporated in the sediment in the world ocean, the following two options were considered.

Approaches #1 and #3
The global data set of particulate organic carbon (POC) in surface sediments published by Seiter et al. (2004) was applied. In order to have a more complete spatial coverage, the value of 1 wt. % was used for data gaps in the Southern Ocean and the south western Pacific. Since a global grid for averaged POC surface concentrations over the Quaternary is currently not available, we apply the Holocene POC data also for the mean Quaternary Approach #3. We, thus, neglect changes in ocean productivity, POC rain rate and POC burial efficiency over glacial-interglacial cycles in our mean Quaternary estimate, which may bias the final outcome, but at present seems to be one of the few practicable ways to account for variable POC accumulation rates in the geologic past.

Approach #2
The above problematic is even more pronounced, if we consider the long-term accumulation of POC. To approach this problem we used a recently published equation by Marquardt et al. (2010) that relates POC concentration to burial velocity to calculate POC contents over geological times: wherew is the averaged sedimentation rate calculated using Eq.
(3). The resulting global distributions of the POC are shown in Fig. 3. The global POCar was calculated applying Eq. (6), as an intermediate result of each scenario. with POCar in gC m −2 yr −1 and POC in wt. %; w is the burial velocity in cm yr −1 ; d S is the density of dry sediment, i.e., 2.5 g cm −3 , and ϕ is the porosity. A porosity value of 0.8 for surface sediments is used for the Holocene and Quaternary scenarios, while 0.5 is used as averaged porosity over the sediment column for the average sedimentation rate scenario.

Thickness of the gas hydrate stability zone
The global thickness of the GHSZ was calculated applying the thermodynamic model provided by Tishchenko et al. (2005) to global grids of bathymetry (Amante and Eakins, 2009), bottom-water temperature, and geothermal gradient. A global grid of averaged bottom-water temperatures over the last 20 yr using the Ocean General Circulation Model (OGCM) (Barnier et al., 2006) was kindly provided by A. Biastoch. Geothermal gradients were calculated based on heat flow models by Pollack (1993) and Hamza et al. (2008). In general, both distribution maps of the heat flow show a similar pattern. The absolute values of heat flow anomalies of Hamza et al. (2008) are generally lower, with values up to 150 mW m −2 compared to values up to 350 mW m −2 presented by Pollack et al. (1993). Regions of major discrepancies between both approaches are located along major midoceanic ridges (e.g. south-east Indian Ocean; South Mid-Atlantic Ridge), where GH are not forming anyhow because heat flow is typically too high and POC very low. Although the estimation of thermal conductivities has been a focus of major interest for the petroleum industry, there is still not a unique way to calculate them. Among the number of means that can be used to calculate the weight of each component (minerals and fluids) in the total thermal conductivity of sediments, the geometric mean is the most accepted one (e.g. Clauser and Huenges, 1995).
where λ is the thermal conductivity in W m −1 K −1 , "min" is mineral, "pw" is pore water, and ϕ is porosity. Thermal conductivity depends on mineralogy and increases with sediment depth due to compaction and cementation (Wallmann et al., 2012). Considering a thermal conductivity of 0.9-3 W m −1 K −1 for clays (Fjeldskaar et al., 2009), 0.6 W m −1 K −1 for pore water in marine sediments (e.g. Deming et al., 1989), and an average porosity of 0.5 over the sediment column, the thermal conductivity of sediments was estimated at 0.7-1.34 W m −1 K −1 . For simplicity, and following previous works, a value of 1 W m −1 K −1 was used in our calculations (e.g. Wood and Yung, 2008;Reitz et al., 2007;Jung and Vogt, 2004;Wallmann et al., 2012). A density of 1000 kg m −3 for water was considered to estimate the hydrostatic pressure at each grid point. The L GHSZ was reduced to the total oceanic sediment thickness (Divins, 2010) where this value was exceeded.

Results and discussion
Approach #1 of POC and SR was applied to calculate the GHI for a Holocene scenario, while Approach #2 of POC and SR was applied to calculate an averaged GH inventory over geological times (Table 1). Additionally, the Holocene data set of POC and the SR grid by  were applied to calculate the GHI for the Quaternary (Approach #3).

Sedimentation rate
The results obtained for the sedimentation rate during the Holocene (Fig. 2) show that the highest values (> 100 cm kyr −1 ) are restricted to the shelf of the Circum-Arctic realm and other shallow water areas, where marine GH are not stable. Low sedimentation rates are characteristic for deep oceanic basins, while intermediate values are common in the margins of Africa, America, and southern Asia, as well as in interior basins such as the Mediterranean Sea.
In contrast, the average of the sedimentation rate during geological times shows lower values along most of the margins, reaching high values (> 100 cm kyr −1 ) in restricted basins off South-eastern Asia, such as offshore Philippines and Java. Relatively high sedimentation rates are also calculated for the Oregon, the southern Argentinean, and the Indian margins, where gas hydrates have been predicted or drilled (e.g. Tréhu et al., 2004;Collett et al., 2008). The average SR calculation does not consider any hiatuses or erosion events that are not included in the total thickness used for the calculation. Specifically, this may lead to a decrease in the SR values along the slopes that cannot be resolved on a global scale.
The average Quaternary scenario is characterized by very high sedimentation rates at the continental slope and rise, since the sediments at the shelf were eroded and transported Table 1. Methane hydrate inventory calculated by two different approaches and equations (see the text for details). GHI is the GH formed by methane generated through degradation of organic matter in the GHSZ at normal compaction conditions, while GHI full comp is the global result in complete compacted sediments. The global integrated POCar for each scenario is also shown for comparison.

Scenario
Global

Particulate organic carbon
In general, POC-rich sediments accumulate in areas of high marine productivity. POC remineralization and preservation strongly depends on environmental conditions such as sedimentation rate, terrigeneous input, bottom-water oxygenation, etc (e.g. Cai and Reimers, 1995;Hartnett et al., 1998;Hartnett and Devol, 2003;Seiter et al., 2005;Burdige, 2007). The two approaches to calculate POC presented here resulted in different distribution patterns (Fig. 3) and show a total range of about one order of magnitude for global POCar (Table 1). Absolute values of POC range from 1 × 10 −2 wt. % in the abyssal plains of the Pacific to ∼ 8 wt. % in rapidly accumulating continental margin settings.
In Approach #1, high POC concentrations of more than 5 wt. % are shown in the Equatorial Pacific, with elevated values of ∼ 2 wt. % towards the central Pacific Ocean. Other relatively high POC concentration areas are offshore Namibia-Angola, Chile, Oman, as well as in the Arctic Ocean. The estimate of Holocene POC accumulation obtained by Approach #1 (1.37 × 10 14 g C yr −1 , see Table 1) is very close to a previous estimate derived from POC concentration data and Holocene sedimentation rates collected by Russian scientists (1.40 × 10 14 g C yr −1 , Baturin, 2007). This conformity strongly supports the validity of Approach #1 to calculate sedimentation rate and POC input during the Holocene.
The POC map calculated after the averaged sedimentation rate over time (Approach #2, Eq. 5) shows a very homogeneous distribution, with values of up to 3 wt. % only. However, it suffers from a lack of data for large continental margin areas such as the Arctic Ocean, the Mediterranean Sea and other marginal seas, where data of the crustal age are not available. These areas coincide with some of the main deposition areas during the Holocene (and likely during the late Neogene), resulting in much lower global sediment inputs as compared to those of the Holocene or Quaternary scenarios (Approaches #1 and #3; Table 1). Due to these restrictions, Approach #2 is clearly not the preferred scenario to calculate the present GH abundance, but can serve as lower boundary estimate.

Thickness of the gas hydrate stability zone
Gridded data of the geothermal gradient used in this study are based on the spherical harmonic approximations of heat flow by Pollack et al. (1993) and Hamza et al. (2008). The spherical harmonic approach does not always fit to measured data, especially in a number of places where crustal flow highly affects the temperature distribution in the sediments (e.g. Costa Rica margin, Grevemeyer et al., 2004). However, both approaches reveal very similar results in terms of GHSZ thickness for most parts of the global ocean. In general, they are also in agreement with recently published data of higher resolution by Wood and Yung (2008), with thinner L GHSZ along the margins and shallower ocean areas, and thicker over the deeper most abyssal plains.
The greatest L GHSZ is found along continental margins with thick sediment cover and where low bottom-water temperatures and low geothermal gradients prevail (such as circumpolar regions and passive margins in the Atlantic and Indian Oceans; Fig. 4). In general, both approaches are also in reasonable agreement with recent reports on BSR depth, which usually have an error of ∼ 5 % (Table 2, Tinivella et al., 2011;Hornbach et al., 2012) However, the comparison of regional BSR depths at various ODP sites (Marquardt et al., 2010, references in Table 2) generally revealed a better match to Hamza's approach than to Pollack's. Therefore, Hamza's approach for heat flow calculation will be applied to calculate global gas hydrate abundances in the following sections.
The global volume of sediment within the GHSZ is calculated to be 59.8 and 66.9 × 10 15 m 3 , for Pollack's and Hamza's data set of heat flow, respectively; from this, only 10.6 to 11.7 × 10 15 m 3 are found at continental margins (8.8 to 9.6 × 10 15 m 3 at passive margins, 2.3 to 2.6 × 10 15 m 3 at active margins), where organic matter accumulation mainly occurs. Applying Hamza's heat flow values and the normal compaction law at sediments published by Wallmann et al. (2012), we estimate the total available pore space within the GHSZ at the margins at 7 to 7.8 × 10 15 m 3 (6 × 10 15 m 3 at passive margins, 1.8 × 10 15 m 3 at active margins). This  value is corroborated by the study of Dickens (2001b), who derived a value of 7 × 10 15 m 3 by multiplying a generic cross section area of a continental margin by the approximate global length of the margins.
A major uncertainty in our approach is the value of the thermal conductivity of marine sediments. The assumed value of 1 W m −1 K −1 is in good agreement with direct thermal measurements of surface and shallow marine sediments at a number of continental-margin sites (e.g. Kaul et al., 2006). For comparison, a maximum thickness of the GHSZ of up to 900 m is calculated applying a thermal conductivity of 1.5 W m −1 K −1 , which can be considered as a maximum value. The distribution maps of the thickness of the GHSZ in this case (not included here) are in good agreement with previous published results by .

Global gas hydrate inventory I
In general, resulting estimates of the global gas hydrate inventory for Approaches #1 and #2, applying two different heat flow models (see Sects. 3.1, 3.2 and 3.3) and compaction settings, are very low (< 10 GtC; Table 1) compared to most previous estimates (cf. Fig. 1). The results are only similar to those of Burwicz et al. (2011) which were based on the same kinetic approach, and similar gridded data sets as Approach #1 for POC and SR. Differences with their result are attributed to offsets between the numerical model and the transfer functions, a different compaction law (Wallmann et al., 2012), and a different thermal conductivity (1.0 W m −1 K −1 instead of 1.5 W m −1 K −1 ) that increases considerably the total thickness of the GHSZ. However, the generally close match between our results (< 10 GtC) and those of 4 GtC) confirms the validity of the transfer functions as useful and simple tools to calculate GH inventories, based on geochemical modelling.
The global gas hydrate inventories predicted for Approach #1 are lower than the values for Approach #2, even though the global POC concentrations and sedimentation rates over the Holocene are higher than the corresponding values averaged over geological times (Table 1). This result is due to the fact that high Holocene sedimentation rates are restricted to the continental shelf areas, i.e. outside the pressure and temperature conditions where gas hydrates can form. As expected, the choice of the heat flow grid (Hamza's or Pollack's) does not have a major effect on the calculation.
Small differences in the POC and SR distributions are responsible for the differences in the GH concentration calculated with Approaches #1 and #2 (Fig. 5, complete compaction, Hamza's heat flow). For the corresponding distribution map for Quaternary conditions, the reader is referred to Fig. 31 in Wallmann et al. (2012). The functions only predict the occurrence of significant amounts of GH at continental margins with high POC and SR (e.g. Blake Ridge, Colombia, Angola-Namibia, Arctic Ocean), including both passive and active margins. In most of these regions GH have been previously recovered or inferred by geophysical detection of a BSR (e.g. Lorenson and Kvenvolden, 2007). In addition, the distribution map considering average sedimentation rates over geological times, shows several grid points in the Indonesia area with values of up to 75 kg m −2 , which ultimately produce the difference in global amounts between both approaches.
The occurrence of GH amounts > 10 kgC m −2 is restricted to a threshold value for the thickness of the GHSZ of ∼ 200 m, POC of ∼ 1 wt. %, and SR of ∼ 30 cm kyr −1 (Eqs. 1 and 2). These values define both the minimum residence time of POC within, as well as the minimum POC flux into the GHSZ required for the generation of sufficient amounts of methane, as well as the residence time of methane in the GHSZ to form GH. A POC content of 0.4-0.5 wt. % has been previously identified as the minimum POC content required for GH formation (e.g. Collett, 1995;MacDonald, 1990;Klauda and Sandler, 2005). Most of the values calculated for GHSZ thickness, POC and SR over the world oceans are below these thresholds, which is ultimately the reason for obtaining only a few isolated patches of significant GH concentrations in Fig. 5.
Obviously, this result corresponds in no way to observations in natural systems (cf. compilations by Lorenson, 2001 or Maslin et al., 2010). The major reason for this mismatch is that the model approach oversimplifies the conditions in complex geological systems, at least with respect to transport processes. The numerical model, which formed the basis for the transfer functions developed by Wallmann et al. (2012) only considered sedimentation, diffusion, and compaction-driven pore water advection (not producing a net upward flow). However, major GH accumulations are related to enhanced fluid and gas advection (e.g. Torres et al., 2004;Buffett and Archer, 2004;Haeckel et al., 2004;Milkov, 2005), which clearly shows the need for an improved transfer function in order to produce more realistic estimates of GH inventories.

Derivation of an extended transfer function considering fluid flow
In order to systematically analyse the effect of upward directed fluid flow on gas hydrate formation, we performed a new series of model runs with the original numerical model (Wallmann et al., 2006). The model explicitly considers steady state compaction of the sediment, diffusive and advective transport of dissolved constituents, input and degradation of POC and particulate organic nitrogen (PON) via sulfate reduction and methanogenesis, anaerobic oxidation of methane (AOM), as well as the formation (and adsorption) of NH 4 , dissolved inorganic carbon (DIC) and methane (CH 4 ). The model calculates the solubility of methane in pore water with respect to the stability field of gas hydrates and considers the formation and dissociation of GH, as well as the formation and dissolution of free methane gas in pore water. When methane concentration exceeds the solubility conditions for gas hydrates or free gas, gas hydrates precipitate or free methane gas accumulates in the pore volume. Gas hydrates dissociate when they are buried below the three phase equilibrium curve described by Tishchenko et al. (2005). The upward transport of gas by rising bubbles is not considered in the model. The POC-degradation rate is calculated as a function of POC input considering the decrease of the reactivity of POC with depth, as a combined effect of POC ageing and the inhibition effect of the accumulation of the degradation products (e.g. DIC and CH 4 ) in pore water (Wallmann et al., 2006). It should be noted in this respect that various degradation kinetics have been suggested in previous studies, which may lead to considerable differences in term of gas hydrate formation (e.g. Burdige, 2011). A systematic comparison of some available POC-degradation kinetics (Wallmann et al., 2012), however, showed that our approach resulted in the best fit to available ODP data (see also Marquardt et al., 2010). Input parameters and boundary conditions for the standard model are shown in Table 2. For a complete description of the model, the reader is referred to the publications by Wallmann et al. (2006Wallmann et al. ( , 2012. Model-runs were performed systematically imposing upward fluid flow velocities of 0.005, 0.01, 0.015, 0.02 cm yr −1 , and covering the range of L GHSZ , POC and SR values typically found in continental margin environments (L GHSZ of 155 to 567 m; POC of 0.1 to 5 wt. % and SR of 10 to 135 cm kyr −1 ). In total 363 model runs were performed. For all model runs the concentrations at the upper boundary of the dissolved species SO 4 , DIC, CH 4 and NH 4 were fixed accordingly to values of the standard seawater composition (Dirichlet conditions), while the concentration gradients of dissolved species were set to zero at the lower boundary of the model domain situated 20-30 m below the GHSZ. A series of model results for fluid flow velocities (FF) ranging from 0 to 0.02 cm yr −1 and constant POC, SR and thickness of the GHSZ (3 wt. %, 100 cm kyr −1 and 567 m, respectively) are shown in Fig. 6. Dissolved methane concentrations increase in response to fluid flow. However, the increase of reaction products in pore water (DIC, CH 4 ) causes a slight decrease of POC degradation with increasing fluid advection, but due to the methane transported from below the GHSZ, the total rates of GH formation increase with increasing fluid flow.
The systematic analysis shows that variations in the four key parameters, thickness of the GHSZ, POC, SR and upward fluid advection velocity (FF) significantly affect the amount of GH accumulation (Fig. 7). For constant FF (0.02 cm yr −1 ), GH amount increases with POC and L GHSZ (Fig. 7a, b), while an increase in SR tends to decrease the gas hydrate inventory (Fig. 7c) because gas hydrate formation is favoured by a high residence time of POC and methane in the GHSZ, it is under high L GHSZ and low SR. However, SR plays a dual role because it also governs the input of POC into the sediment in the GHSZ. Thus, an increase in SR may enhance the gas hydrate inventory if the POC input is the major limiting factor for the gas hydrate formation. Since additional methane is transported into the GHSZ, the gas hydrate inventory is generally enhanced by upward fluid flow (FF). However, FF also intensifies the transport of dissolved methane into the sulphate-bearing surface sediments where methane is consumed by microbial oxidation processes. Moreover, methane production rates within the GHSZ are reduced by FF due to the accumulation of DIC and CH 4 in ambient pore fluids.

E. Piñero et al.: Estimation of the global inventory of methane hydrates
For slow SR, low POC and thin GHSZ, fluid advection has a positive effect on the GH accumulation, whereas for relatively thick GHSZ, high POC and fast SR, it has a negative effect. The latter aspect is further evaluated in the plots in the second row of Fig. 7, which show the results of a series of model runs considering the variation of fluid velocity, for constant SR (50 cm kyr −1 ) and POC (3 wt. %), L GHSZ (567 m) and POC (3 wt. %), and SR (50 cm kyr −1 ) and L GHSZ (567 m). For constant SR and POC, GH accumulation increases with L GHSZ (Fig. 7d). For a thin GHSZ (e.g. 266 m, orange dots in Fig. 7d), the accumulation of GH is enhanced by a faster fluid flow that increases the methane input into the GHSZ. However, GH formation decreases with FF for thicker GHSZ (e.g. 567 m, green dots in Fig. 7d) because the extensive upward migration of fluid produces an increase in DIC and dissolved CH 4 promoting methane loss to the shallow sulphate-bearing sediments. The sedimentation rate shows a similar behaviour with increasing fluid flow, for constant L GHSZ (567 m) and POC (3 wt. %) (Fig. 7e). For low sedimentation rates (10-20 cm kyr −1 , yellow and orange dots), GH accumulation increases with fluid advection, whereas it decreases for fast sedimentation rates (red, green and blue dots). For constant SR and thickness of the GHSZ, FF has a dual effect as well (Fig. 7f). For low POC contents, FF enhances the GH formation, whereas for organic matterrich sediments, it has the opposite effect. Increasing FF high enough would inhibit GH accumulation, due to the decrease in residence time within the GHSZ. Thus, for each combination of SR, POC and L GHSZ there exists a limit of the fluid velocity where GH do not accumulate any more.
The results of the parameter analysis (363 model-runs) could be fitted by empirical functions (Eqs. 8a-8e) indicated by lines in Fig. 7a-f. As a result of this multi-parametric nonlinear least-squares fitting procedure, the critical points that limit the two behaviour domains for upward fluid flow are described as For the domain FF ≥ 0.0001 SR (2 + ln [POC]), the following function was derived and for the domain FF < 0.0001 SR (2 + ln [POC]), provides an appropriate representation. In addition, gas hydrate accumulation is inhibited for: and overall The 9 coefficients were determined in a combined non-linear least-squares fit of Eqs. (8) and (1)  To further test the accuracy of the extended transfer functions (Eq. 8), we calculated the GH content using the same parameters as for the numerical model runs. The extended transfer function reproduces the modelled data very well. Most data points plot along the 1 : 1 correlation line (Fig. 8). The general scatter is low, with a correlation coefficient (σ 2 ) of 0.962 and a variance of residuals of 1167.4 kg m −2 . The linear correlation matrix of fit parameters can be found in Table A3. Significant deviations occur for a combination of low SR, thin L GHSZ , and low POC values (SR < 35 cm kyr −1 , L GHSZ < 300 m and POC < 1.5 wt. %).

Global gas hydrate inventory II
While global data sets for L GHSZ , SR and POC are available (although subject to considerable uncertainties, cf. Sects. 3.1, 3.2 and 3.3) such information is not available for FF. The most important mechanisms generating fluid migration at continental margins are the mechanical reduction of porosity (compaction, compression), and diagenetic and metamorphic processes (release of water from breakdown of hydrous minerals, gas production; e.g. Kastner, 1999;Kopf, 2002;Wallmann et al., 2012). As mentioned above, normal compaction does not directly lead to a net upward fluid flow, and hence values for FF cannot be easily derived from empirical equations of porosity reduction. In fact, fluid velocities have only been calculated for a limited number of areas where focused fluid flow preferentially occurs (e.g. mud volcanoes, cold seeps, etc.), ranging from few mm yr −1 to several m yr −1 (e.g. Wallmann et al., 1997;Linke et al., 2005). However, because these fluid pathways are neither homogenously distributed nor representative for larger areas of the seafloor, those values for FF cannot be extrapolated.
Recently published papers  used the tracers SO 2− 4 and 129 I as proxies to constrain fluid flow in ideal passive and active continental margins. However, the concentration of these species has been only measured at a few locations along the margins, and therefore, their global distribution is not available. A rough but valid method is to approach the global upward FF at continental margins by mass-balance considerations based on a steady-state pore water budget . Due to compaction, porosity of sediments LGHSZ=155 m LGHSZ=266 m LGHSZ=408 m LGHSZ=567 m LGHSZ=155 m LGHSZ=266 m LGHSZ=408 m LGHSZ=567 decreases with increasing sediment depth, which can be described by empirical relationships (e.g. Einsele, 2000). However, a significant volume of pore water remains within the sediment after normal steady-state compaction. The global rate of pore water burial after compaction (PW f in cm 3 yr −1 ) may be estimated from the global sedimentation rate (SR glob ∼ 20 × 10 15 g yr −1 , derived from the bathymetric grid), the mean density of dry solids (d S = 2.5 g cm −3 ) and the porosity after compaction (ϕ f = 0.2-0.3): The resulting value for PW f is 2.0-3.4 km 3 yr −1 . Independent estimates of upward fluid flow at active continental margins yield a range of 1 to 2 km 3 yr −1 , (Von Huene and Scholl, 1991;Moore and Vrolijk, 1992;Jarrard et al., 2003) indicating that a major part of the buried pore fluids is mobilized by tectonic over-pressuring.
Taking into account the extension of active margins (Fig. B1), and a porosity at the upper limit of our model of 0.75 (excluding the uppermost bioturbated sedimentary column, Wallmann et al., 2012), fluid velocities of 0.01 and 0.02 cm yr −1 were obtained at active margins for 1 and 2 km 3 yr −1 , respectively. These values are in the same range as previously published averaged values (e.g. Buffett and Archer, 2004), and fluid flow values determined for diffusive systems in the Cascadia margin (0.017 cm yr −1 ; Malinverno et al., 2008). Thus, the flux of dissolved methane into the GHSZ may be estimated in 1.7-11.5 Tg CH 4 yr −1 , assuming a concentration in the pore fluids of 55-200 mmol CH 4 L −1 , and taking into account the porosity reduction due to burial. Fluid flow at active margins is not only driven by vertical compaction, but also by the subduction and de-watering of sediments transported into the subduction zones in the incoming oceanic plates. The new transfer function (Eq. 8) was thus applied to calculate the inventory of gas hydrates at active margins assuming normal compaction with additional imposed fluid flow velocities of 0.01 cm yr −1 (1 km 3 yr −1 ) and 0.02 (2 km 3 yr −1 ). The transfer functions previously derived by Wallmann et al. (2012) were applied to calculate the gas hydrate inventory at passive margins considering both normal compaction (Eq. 1) and complete compaction (Eq. 2) with no additional fluid inflow from below. The resulting ranges of estimates of the global inventory of GH are listed in Table 4, for the different scenarios to calculate POC and SR presented above, as well as normal or full compaction at passive margins, and two different fluid velocities at active margins.
There are considerable differences in the predicted global mass of GH with respect to the chosen method of calculation. Generally, Approaches #1 and 2 result in much lower estimates (3 to 59 GtC) than Quaternary sedimentation rates (Approach #3). For Holocene conditions major concentrations of GH (up to 515 kg m −2 ) are only predicted in the Peruvian margin, while for the average SR enrichments of up to 105 kg m −2 are predicted in the Chilean margin, as well as in several areas of the Indonesian margin. However, for most of the margins where low SR occur, no GH are predicted evidencing the high impact of SR conditions on the gas hydrate accumulation.
Enhanced fluid flow at active margins has a major effect on the total GH inventory. Under Quaternary boundary conditions (Approach #3), the global gas hydrate inventory continuously increases in response to upward fluid flow at active continental margins. However, the other two approaches (#1 Table 4. Global gas hydrate estimations (GH FF ) in GtC, for the different approaches of SR and POC input (Holocene, Quaternary and the average over geological times). The results of different degrees of compaction at passive margins (normal or full), as well as of different fluid velocities at active margins (0.01 and 0.02 cm yr −1 ) are shown. For all calculations, Hamza's heat flow values were used to estimate the thickness of the GHSZ. At active margins, a volume of 1 km 3 yr −1 of expelled water corresponds to a fluid velocity of 0.01 cm yr −1 , and 2 km 3 yr −1 to 0.02 cm yr −1 . and #2) show a more complex response to fluid flow. FF has a positive effect up to 0.01 cm yr −1 while beyond this averaged fluid velocity, the total GH accumulation decreases (Table 4). This interesting response can be attributed to the fact that upward fluid flow can either favour or inhibit gas hydrate formation depending on local conditions (see Sect. 3.5). Complete compaction of sediments deposited at passive margins produces a significant increase in the gas hydrate inventory as previously shown in Wallmann et al. (2012). As already discussed above, the use of long-term averages of sedimentation rates (Approach #2) can only be regarded as minimum estimates, because large areas underlain by continental crust are neglected. Hence, we consider the simulations for Quaternary conditions (Approach #3) as the most realistic and estimate the global sea floor inventory of GH in ∼ 550 GtC. Figure 9 shows the resulting GH distribution map for Quaternary sedimentation rates (Approach #3) using Hamza's heat flow model, complete compaction at passive margins, and a fluid flow velocity at active margins of 0.02 cm yr −1 . The map shows that minimum concentrations of 10 kg m −2 stored in GH are typical for extended margin areas. Maximum concentrations above 100 kg m −2 are predicted for isolated patches at active continental margins off Chile, Peru, Argentina, Alaska, as well as Japan, Taiwan and Indonesia, with up to 157 kg m −2 predicted for the Gulf of Oman, which is characterized by high POC and SR and a thick GHSZ.
Considering the total available pore volume in the GHSZ, we estimate a global averaged GH saturation of 0.074 %. This saturation is much lower than GH measured in several GH-rich settings, with saturations averaged over the GHSZ ranging from 1 to 10 % (e.g. 1-2 % at Hydrate Ridge, Tréhu et al., 2004;2-4 % at Blake Ridge, Borowski, 2004), but, as outlined before, higher GH occupancy is typically due to specific geo-tectonic settings with enhanced upward flow of methane-rich fluids and rise of free gas.
Although our estimate is in the range of some previous estimates, and hence may be simply regarded as another method to come to the same result, we believe that this approach implies some clear improvements to previous studies: (i) The systematic analysis of the model-runs results reveals a clear distinction of the control parameters of GH formation. (ii) We used improved, state-of-the-art gridded data sets for the global calculation of control parameters. (iii) In contrast to some previous studies we used a calibrated kinetic model of subsurface methane generation and a valid mass balance approach for calculating FF at active continental margins. (iv) We evaluate the effect of porosity reduction in different compaction settings. (v) We provide a simple and realistic tool to calculate GH concentrations applicable to any marine-geological environment. However, although we were able to provide useful estimates of averaged values for FF at active margins, fluid advection still remains a critical parameter. Fluid flow is not homogeneous, but preferentially follows geological pathways, such as faults and other high-permeability conduits, leading to seepage at the seafloor and generating local sub-seafloor gas hydrate enrichments (e.g. . Furthermore, it plays a major role controlling distinct GH accumulations at passive margins (e.g. Blake Ridge, Dickens et al., 2001c). As a result of using an averaged fluid flow distribution, our global map of gas hydrate occurrence shows a smooth distribution along the continental margins. Due to the chosen method, isolated spots of significant GH accumulations such as found at Blake Ridge or Hydrate Ridge cannot correctly be reflected/considered (Fig. 9). We are aware about the fact that the "real" gas hydrate distribution may have a more patchy appearance, with large areas with only minor concentrations of disperse gas hydrate and patches with rich accumulations. However, most of the control parameters influencing the migration pathways cannot be resolved in a 1 × 1-degree resolution. For example, it is known from accretionary margins that fluid flow changes with the distance from the wedge (e.g. Breen and Orange, 1992;Saffer and Tobin, 2011), and hence may cause significant regional variations in the GH inventory. Consequently, any type of regional heterogeneity not resolved in the database will introduce a fairly unknown degree of uncertainty, which clearly shows the need for improved fluid budgets along continental margins. Furthermore, errors in any database may have different weight depending on their latitudinal position, because the areas of each grid cell (resolution 1 • × 1 • ) vary as a function of latitude. As example, an error in the Arctic Ocean would have an impact in the global inventory ∼ 110 times smaller than an error near the Equator.
Similarly, another critical issue is that all model runs were performed assuming steady-state conditions on timescales of several millions of years (Wallmann et al., 2012) assuming major control parameters (heat flow, POC, SR, FF) were constant. Hence, the model runs as well as the transfer functions do not consider any temporal variations of the sedimentary system. However, improved simulations in this regard are likely impossible as they by far exceed available data sources (at least on a global scale).
Overall, we believe that our study provides an improved and well constrained estimate for the global inventory of GH in marine sediments of ∼ 550 GtC. Since the ascent of free methane gas, as another important, but completely unconstrained parameter could not be included one may still consider this value as a minimum estimate.

Conclusions
Our study shows that the global inventory of marine GH can be estimated by the application of transfer functions that were derived from systematic runs of a numerical transportreaction model. Essentially, GH formation in marine sediments is determined by the thickness of the GHSZ, the POC, SR, and upward fluid advection. Without FF, only 3 to 455 GtC would be stored in marine GH, depending on the considered SR conditions. GH would only form in areas of extremely high organic matter input, such as the margins of Chile, Peru and Central America, as well as in the Arabian Sea and the Arctic Ocean. The new transfer function developed in this study explicitly considering FF (0.01-0.02 cm yr −1 at active margins) predicts minimum GH concentrations of ∼ 10 kg m −2 for extended continental margin areas around the globe, with significantly increased concentrations along the American Pacific margin, as well as in Japan, Taiwan, Indonesia and the Arabian Sea. The global inventory is estimated to be in the range of ∼ 550 GtC. Our results suggest that the ascent of methane-rich fluids is a dominant process controlling GH formation and that this is not restricted to sites of focused fluid advection. Due to missing constraints on the upward migration of free methane gas into the GHSZ we still consider the value above as a minimum estimate.
Appendix A Table A1. Symmetric linear correlation matrix of fit parameters of Eq. (1) (normal compaction, no fluid flow). The number of fitted simulations is n = 104; χ 2 = 11315.7; variance of residuals is 115.5 kg m −2 ; linear correlation coefficient for GHI from model simulations and fitted equation is r 2 = 0.9971.