Volume 29, Issue 5 p. 610-625
Research Article
Free Access

Modeling the fate of methane hydrates under global warming

Kerstin Kretschmer

Corresponding Author

Kerstin Kretschmer

GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

Now at MARUM - Center for Marine Environmental Sciences and Faculty of Geosciences, University of Bremen, Bremen, Germany

Correspondence to: K. Kretschmer,

[email protected]

Search for more papers by this author
Arne Biastoch

Arne Biastoch

GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

Search for more papers by this author
Lars Rüpke

Lars Rüpke

GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

Search for more papers by this author
Ewa Burwicz

Ewa Burwicz

GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

Search for more papers by this author
First published: 30 March 2015
Citations: 98

Abstract

Large amounts of methane hydrate locked up within marine sediments are vulnerable to climate change. Changes in bottom water temperatures may lead to their destabilization and the release of methane into the water column or even the atmosphere. In a multimodel approach, the possible impact of destabilizing methane hydrates onto global climate within the next century is evaluated. The focus is set on changing bottom water temperatures to infer the response of the global methane hydrate inventory to future climate change. Present and future bottom water temperatures are evaluated by the combined use of hindcast high-resolution ocean circulation simulations and climate modeling for the next century. The changing global hydrate inventory is computed using the parameterized transfer function recently proposed by Wallmann et al. (2012). We find that the present-day world's total marine methane hydrate inventory is estimated to be 1146 Gt of methane carbon. Within the next 100 years this global inventory may be reduced by ∼0.03% (releasing ∼473 Mt methane from the seafloor). Compared to the present-day annual emissions of anthropogenic methane, the amount of methane released from melting hydrates by 2100 is small and will not have a major impact on the global climate. On a regional scale, ocean bottom warming over the next 100 years will result in a relatively large decrease in the methane hydrate deposits, with the Arctic and Blake Ridge region, offshore South Carolina, being most affected.

Key Points

  • Methane hydrates are vulnerable to climate change
  • The release of methane into the water column is limited within the next century

1 Introduction

The Earth's climate system has experienced significant changes over the last century. Although strongly masked by natural variability on synoptic, interannual, and decadal time scales, it is clear that the long-term temperature trend is due to anthropogenic influences [Huber and Knutti, 2012; Santer et al., 2013]. Observational records reveal that not only the atmosphere but also the global ocean is subject to a warming trend. According to the fifth Assessment Report of the Intergovernmental Panel on Climate Change [Intergovernmental Panel on Climate Change (IPCC), 2013], the upper ocean (0–75 m) temperature has risen from 1971 to 2010 by 0.11°C/10 years, the temperature at 700 m by 0.015°C/10 years. Few observations exist on the deep ocean (here referred to as the regions below 700 m water depth) and warming estimates largely rely on climate models. Owing to its vast volume, the deep ocean is likely to play a key role in modulating the evolution of global warming since it effectively removes heat from the atmosphere and surface ocean [Meehl et al., 2011; Lamarque, 2008].

With observations and climate projections pointing to significant ocean warming over the next 100 years, rising seafloor temperatures may lead to the destabilization of gas hydrates within marine sediments.

Methane (CH4) is the most common gas found in naturally occurring gas hydrates, accompanied by hydrogen sulphide, carbon dioxide (CO2), and, less frequently, by higher hydrocarbons [Sloan, 1990]. The pressure- and temperature-controlled stability limit of methane hydrates is such that the base of the gas hydrate stability zone (GHSZ) outcrops on the upper continental slope at around 300–500 m water depth (Figure 1). Given that most ocean warming will occur within the first few hundred meters water depth, it will be this region where hydrates are most susceptible to dissociation from rising bottom water temperatures.

Details are in the caption following the image
Schematic of a typical continental shelf margin depicting the present-day area of gas hydrate stability (light + dark orange). Within this zone, the accumulations most vulnerable to climate change are marked in dark orange (adapted from Boswell and Collett [2011]). In boxes A and B, a detailed depiction of how the thickness of the GHSZ may change under climate change is shown. The blue dashed lines resemble the phase boundary for methane hydrate, calculated for an exemplary sediment deposited at 500 m water depth and sulfate-free pore water with a salinity of 35. The solid and dashed dark red lines depict the steady state temperature profiles under present-day conditions and a warming scenario of either 3.0°C (box A) or 1.5°C (box B) (red arrows). On the upper continental slope (300–500 m water depth), the point where the GHSZ outcrops on the seafloor is shifted down the slope, which leads to complete dissociation (box A), whereas the lower continental slope and continental rise may face a reduction from the bottom (box B). This all happens on longer time scales, the transient temperature profiles in the sediment are indicated by green dotted lines (after 50, 100, 150, 200, and 500 years from left to right). The base of the GHSZ for present-day is denoted by the Bottom Simulating Reflector (BSR).

Figure 1 illustrates two possible scenarios resulting from a gradual increase in bottom water temperature over the next 100 years, followed by a several hundred years long thermal equilibration phase. If the temperature rise is sufficient to eventually result in the complete dissociation of the hydrate column, progressive hydrate melting will occur from both the top as well as bottom of the GHSZ until the entire stability zone vanishes (box A in Figure 1). This is the main process by which hydrates are likely to destabilize within the next century. However, due to the sluggish thermal diffusion rate in marine sediments, complete dissociation will only occur after several hundred years. If seafloor warming will only reduce the thickness of the stability zone, hydrate melting occurs at the bottom of the stability zone (box B in Figure 1). This, in turn, is the likely process by which the global hydrate inventory is reduced on centennial to millennium time scales in response to global warming. Since the stability of gas hydrates is more vulnerable to small changes in temperature than to pressure [Tishchenko et al., 2005], we concentrate on the influence of varying bottom water temperatures only. Hunter et al. [2013] recently found that even large changes in the sea level (of up to +20 mm/yr) do not lead to significant variations in the global gas hydrate inventory.

During the dissociation of methane hydrates, a free methane gas phase forms that may escape to the seafloor or even the atmosphere, which has important regional and global consequences. Through microbial aerobic oxidation [Valentine et al., 2001; Murrell, 2010], the water column can be affected by oxygen depletion and ocean acidification [Biastoch et al., 2011]. Only if enough methane passes through this microbial filter and reaches the atmosphere, its global warming potential (about 25 times higher than CO2 [IPCC, 2007]) could lead to a further acceleration of climate change. As a consequence, it has been suggested that hydrates played an important role in past and future climate change.

The most prominent and best studied past warming event is the Paleocene-Eocene Thermal Maximum. Enormous amounts of carbon (>2000 Gt) were released into the atmosphere causing a surface temperature rise of ∼6°C [Cohen et al., 2007; Higgins and Schrag, 2006]. The associated carbon isotope excursion points to an isotopically light source, which is consistent with methane released from dissociating gas hydrates [Dickens et al., 1995]. Doubts whether sufficient amount of hydrate were stable under warmer Paleocene marine conditions and the discrepancy between the amount of warming and the carbon excursion pose significant challenges to this hypothesis [Higgins and Schrag, 2006]. Methane hydrate dissociation, according to the clathrate gun hypothesis, may have also played a role during Quaternary climate change [Kennett et al., 2003] but also this view has been convincingly challenged based on ice core data [Sowers, 2006].

Fueled by field studies that indicate an increase in methane fluxes from submarine Arctic permafrost and the seafloor [Westbrook et al., 2009; Shakhova et al., 2010], methane hydrates are also debated in the context of contemporaneous climate change. The sluggish rate of heat diffusion into marine sediments renders catastrophic and widespread hydrate dissociation triggered by warming bottom waters unlikely [Archer, 2007; Maslin et al., 2010; Ruppel, 2011; Biastoch et al., 2011] and recently discovered gas flares offshore Svalbard appear to be related to natural seepage processes rather than to global warming induced hydrate dissociation [Berndt et al., 2014].

A robust assessment of the role of methane hydrates in contemporaneous climate change comes down to the question of how large the current global marine gas hydrate inventory is and how much of it will be affected by warming seafloor temperatures. Despite the clear need to answer these questions, the global abundance of methane hydrates in marine sediments remains poorly constrained. Global estimates range over several orders of magnitude from 500 to 55,800 Gt C (cf. Figure 2) without clear convergence [cf. Kvenvolden and Lorenson, 2001; Buffett and Archer, 2004; Milkov, 2004; Klauda and Sandler, 2005; Archer et al., 2009; Boswell and Collett, 2011; Burwicz et al., 2011; Wallmann et al., 2012; Piñero et al., 2013]. Estimates based on the global extrapolation of known hydrate deposits [e.g., Milkov, 2004] suffer from the sparse global data coverage, while the uncertainty in model-based estimates [e.g., Klauda and Sandler, 2005; Archer et al., 2009; Burwicz et al., 2011; Wallmann et al., 2012] results from differences in model assumptions and variable quality of global input data sets.

Details are in the caption following the image
Global estimates of methane hydrate inventories (Gt C) with error bars indicating estimated ranges of the Gas Hydrate Inventory (GHI).

Similar uncertainty exists for the projected warming of bottom waters. The observation-based average trends in warming given in the latest IPCC report are global averages, which are not necessarily representative for specific regions and problematic to extend into the future. We therefore have to rely on projections based on climate models to assess the spatial and temporal evolution of bottom water temperatures. However, climate models are currently run at a lateral resolution of 1–2° and a vertical resolution of tens to hundreds of meters, which may introduce interpolation errors into the prediction of bottom water temperatures, especially along the steep continental slopes.

The above considerations illustrate the need for improved tools for the assessment of the current global hydrate inventory and its fate under global warming. In this study, we present a new self-consistent methodology, which is based on the combination of a high-resolution ocean hindcast simulation, future climate projections, and a parameterized reaction-transport model for the prediction of hydrates in marine sediments. One advantage of this approach is that the current distribution of marine gas hydrates, which represents the starting condition for the climate simulations, is fully consistent with both the ocean as well as the climate model and is derived from a process-based analysis and not extrapolated from sparse regional data. The high-resolution ocean hindcast simulation provides the baseline bottom water temperature onto which the warming trends predicted by the climate model are projected. This approach increases the effective resolution and assures the spatial accuracy and consistency with the initial temperature field and global input data sets.

2 Data and Methods

2.1 Model Configurations

For analyzing present-day global climate conditions, data from a global ocean/sea ice model, performed under atmospheric boundary conditions of the past decades (hindcast), were used. The model is based on the Nucleus for European Modelling of the Ocean (NEMO) [Madec, 2008] and consists of an eddy-active, z coordinate global grid at 0.25° nominal resolution (ORCA025) [Barnier et al., 2006]. The effective resolution amounts to ∼13.8 km at 60° latitude and gets coarser with decreasing latitude (being ∼27.75 km at the equator).

With a total number of 46 levels in the vertical, the grid spacing is fine at the surface (10 levels in the upper 100 m) and then increases with depth to 250 m at the bottom. Bottom topography is described by partially filled bottom cells to better represent topographic slopes and providing a reasonable representation of the global oceanic circulation [Barnier et al., 2006]. Vertical mixing is parameterized using the turbulent kinetic energy scheme of Blanke and Delecluse [1993], lateral mixing is aligned along isopycnal surfaces.

The model run is initialized with temperature and salinity fields combined from the Levitus climatology [Levitus et al., 1998] for mid and low latitudes, the Polar Science Center Hydrographic Climatology (PHC 2.1) for high latitudes [Steele et al., 2001] and from the Medatlas climatology for the Mediterranean Sea [Jourdan et al., 1998]. The model integration is spun-up to a quasi steady state using atmospheric fields for the years 1978–2007 then integrated over the hindcast period 1948–2007 [Behrens et al., 2013]. Momentum, heat, and freshwater fluxes at the sea surface are based on the atmospheric data sets developed by Large and Yeager [2004, 2009] and implemented through bulk formulae according to the CORE-II-protocol (Coordinated Ocean-ice Reference Experiments), suggested by Griffies et al. [2009]. These forcing data sets are originally based on the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis products and corrected and adjusted using various observational (e.g., satellite) data sets to correct for known biases [Large and Yeager, 2009]. The forcing fields are provided at 6-hourly (windspeed, temperature, and humidity), daily (short- and long-wave radiation), and monthly (precipitation) resolution. Annual mean river runoff was provided by Dai and Trenberth [2002].

For analyzing future climate change, simulations of the Kiel Climate Model (KCM) [Park et al., 2009], a coupled ocean/atmosphere/sea ice general circulation model, were used. KCM consists of the atmospheric European Centre Hamburg Model (ECHAM5) [Roeckner et al., 2003] and the ORCA2 model configuration, the latter also based on the same numerical framework as the ORCA025 model configuration, but at lower resolution (2° horizontally, 31 vertical levels). Here two experiment series were used: A twentieth century equivalent control experiment and a series of twenty-two 100 yearlong global warming simulations. The control simulation describing “present-day” climate conditions was forced by constant greenhouse gas concentrations (CO2=348 ppm) and was integrated for 1100 years in total. For the global warming experiments (21st century equivalent), the CO2 concentration was increased by 1%/yr (compound) until CO2 doubling is reached after about 70 years and stabilized thereafter for another 30 years. Those experiments were started from different initial conditions chosen semiregularly at 30 to 40 year intervals from the control run. According to Hawkins and Sutton [2009], future predictions of climate change are dominated by model uncertainty on decadal time scales. In order to reduce the model uncertainty of KCM, an ensemble mean is used in this study to provide a robust assessment of temperature changes under future climate change.

2.2 Stability Analysis of Gas Hydrates

The stability of gas hydrates within marine sediments depends on temperature, pressure, and to a smaller extent on salinity. The boundary between the hydrate and free gas phase fields for structure I methane hydrates is obtained via the modified Pitzer approach of Tishchenko et al. [2005]. Pore fluid pressures are assumed hydrostatic, and salinity within the sediment column is assumed to be the same as at the seafloor. The temperature gradient within the sediment column is computed from global heat flow data:
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0001(1)
where qH is the conductive heat flow obtained from Hamza et al. [2008] and κB is the bulk thermal conductivity of sediments. Throughout this study, a constant thermal conductivity of κB=1.5 W/m/K is used, which is based on the characteristic properties of shaly sandstones [Hantschel and Kauerauf, 2010]. With these assumptions, the steady state geotherm can be calculated:
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0002(2)
where TBW is the bottom water temperature predicted by the ocean circulation model. To assess possible gas hydrate dissociation as a consequence of warming bottom waters, a transient thermal solution is necessary. The sedimentary temperature profile is then given by the solution of the heat conduction equation:
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0003(3)
where αB is the bulk thermal diffusivity (in m2/s), which is defined as
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0004(4)
with ρ being the density (in kg/m3) and cp the specific heat capacity at constant pressure (in J/kg/K) of bulk sediments. Furthermore, the volume-specific heat capacity, (ρ·cp)B, can be calculated by
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0005(5)
where Φ denotes porosity (ranges between 0 and 1) and the subscripts s and f represent solids and fluids, respectively. For the calculation of αB the following parameter values have been chosen: κB=1.5 W/m/K, ρs=2650 kg/m3, ρf=1025 kg/m3, urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0006 J/kg/K, and urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0007 J/kg/K. Although porosity changes with depth, a simple constant thermal diffusivity is assumed based on an average porosity of Φ = 0.5. Equation 3 is solved using finite differences and thus the thickness of the GHSZ (LGHSZ) in marine sediments for transient state conditions can be evaluated. All transient calculations start from the steady state solution given by equation 2 and use the predicted changing bottom water temperature as top boundary condition and heat flow qH as bottom boundary condition.
We follow Wallmann et al. [2012] and compute the hydrate inventory, GHI, in kg C/m2 at every grid point with a transfer function that takes gas hydrate stability zone thickness, LGHSZ (m), the particulate organic carbon (POC) concentration at the seafloor (wt %), and sedimentation rate (cm/kyr) as input:
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0008(6)
The following fitting parameters (parameter value ± standard deviation) are used: 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, and f = 0.3925 ± 0.032. The derivation of the transfer function in Wallmann et al. [2012] was done by parameterizing the predicted hydrate inventory computed by a transport reaction model [Wallmann et al., 2006; Burwicz et al., 2011], which resolves kinetics-controlled microbial degradation of POC via sulfate reduction as well as methanogenesis, hydrate formation, and pore fluid chemistry (CH4, DIC, SO4). Equation 6 implicitly assumes a compaction law that describes exponential compaction with depth:
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0009(7)
where Φ0 is the initial porosity at the sediment surface and c0 is the compaction length scale (in m−1).
As sedimentation rates are not available on a global scale, water depth is used as a proxy to estimate them. The evaluation of sedimentation rates of Holocene surface sediments measured at different stations (in total more than 500) showed a decrease in the latter with increasing water depth [Burwicz et al., 2011]. Due to this depth dependency, Burwicz et al. [2011] were able to derive a parameterization that provides sedimentation rate w (in cm/yr) as a function of water depth (z in m):
urn:x-wiley:gbc:media:gbc20265:gbc20265-math-0010(8)
with w1=0.117 cm/yr, w2=0.006 cm/yr, z1=200 m, z2=4000 m, c1=3, and c2=10. Evaluating equation 8 shows that most sedimentation (∼72%) takes place on the continental shelf (z = 0–200 m), where w ranges between 0.01 cm/yr and 10 cm/yr. While this is correct throughout the Holocene, direct observations revealed that the thickness of Holocene sediments is only a fraction of the actual GHSZ thickness. It is therefore more appropriate to consider mean sedimentation rates averaged over a period of several million years for the prediction of hydrate accumulation.

Hay [1994] found that under glacial conditions, as they have repeatedly occurred throughout the Quaternary, the anomalous Holocene shelf accumulation rate (as predicted by equation 8) might have been diminished by an order of magnitude. It is therefore reasonable to shift the main deposition areas from the shelf to the deeper continental slope and rise. Hence, the sedimentation rate was increased over a ∼500 km wide range around the continental margins by moving deposited material from shelf regions (z = 0–200 m) to continental slopes [Burwicz et al., 2011]. At distance greater than 500 km to the coast sedimentation rates are computed via equation 8. It should be noted that the total amount of sediment supplied to the global ocean stays constant. For estimating the GHI of present climate conditions, we use this modified sedimentation mode as today's marine hydrate accumulation rates are likely controlled by Quaternary mean values [Wallmann et al., 2012].

The total global methane inventory for present climate conditions can be obtained by integrating the GHI (given by equation 6) over the surface area for each grid cell of the ocean model. To evaluate the possible impact of future climate change onto marine hydrates, the change in the total global methane inventory is calculated by first using present-day data, then inferring bottom temperature warming due to future projections, whereby the change in the global gas hydrate inventory under global warming is calculated by scaling the present GHI with the variations in the LGHSZ. Note that this approach of scaling the GHI by the relative change in LGHSZ will slightly underestimate the change in the global gas hydrate inventory, as GHI scales exponentially with the exponent b =∼ 1.7 and not linearly. The reason for this lies in the nonlinear pressure- and temperature-dependent solubility limit of methane in seawater. Example simulations with the full transient reaction transport model upon which the employed transfer function is built [Burwicz et al., 2011] show that the time the system takes to reequilibrate to a new steady state is much longer than the here considered time scales of seafloor warming. It is therefore legitimate to assume that the GHI scales linearly with the reduction in the stability zone thickness in the global warming experiments.

For the present-day quantification, global water temperatures and salinities, averaged over the last 20 years of model integration (1988–2007), have been extracted from the lowest water grid box of the ORCA025 model. To evaluate future climate change, the temperature and salinity distributions extracted from the KCM have been three dimensionally interpolated onto the tripolar grid of ORCA025. This allows a more accurate determination of the bottom slopes and values due to the higher horizontal and vertical resolution of ORCA025. Temperature changes and steady state GHSZ changes are calculated from the differences between the 22 individual ensemble members (last 30 years) and their corresponding periods of the reference model. For the calculation of transient GHSZ changes, the ensemble mean was derived after the application of the individual gas hydrate stability analyses.

In many deep sea areas the calculated thickness of the GHSZ is larger than the actual sediment thickness; in these cases the vertical extent was restricted to the sediment thicknesses. For this purpose, a combined data set of global sediment thicknesses based on the data sets provided by the National Oceanic and Atmospheric Administration National Geophysical Data Center for the World's oceans and marginal seas [Divins, 2003] and by Laske and Masters [1997] for the Arctic region was used.

Global POC concentrations in marine surface sediments are based on the combination of the data sets provided by Seiter et al. [2004] and Romankevich et al. [2009]. For some areas (e.g., Southern Ocean and south western Pacific), POC concentrations are not available and a default value of 1 wt % POC is assumed.

3 Results

3.1 Present-Day Climate Conditions

The global map of the thickness of the GHSZ under steady state conditions for present-day climate (Figure 3) shows values of up to 600–900 m in polar regions due to low bottom water temperatures (<0°C) and along continental margins where bottom water temperatures are low and a thick sediment cover exists. Thinner GHSZs (<150 m) are predicted for most of the deeper ocean, where the vertical extent is restricted by thin sediment thicknesses. The shelf regions (above 300 m water depth) feature no stable conditions for structure I gas hydrates because of too high bottom water temperatures (except for permafrost regions but those are not considered in this study). The global volume of sediment within the GHSZ under present-day climate conditions is estimated to be 87·1015 m3. Reassuring is the good agreement with observational values of the depth of the Bottom Simulating Reflector (BSR), which is often used as an indicator for the base of the GHSZ [Haacke et al., 2007], for different regions (Table 1).

Details are in the caption following the image
Global map of the gas hydrate stability zone thickness (LGHSZ in m) under present-day climate conditions (mean 1988–2007). The white shadings indicate areas where LGHSZ=0 m. The dashed contour line marks the 300 m isobath.
Table 1. Thickness of the GHSZ and Depth of the Bottom Simulating Reflectora
Region LGHSZ BSR Depth Reference
Arctic (Beaufort Sea) 452 m 300–700 m Andreassen et al. [1997]
Blake Ridge 447 m 440–600 m Kvenvolden and Lorenson [2000]
Cascadia Margin 202 m 230 m Riedel et al. [2006]
Sea of Okhotsk 331 m 50–900 m Lüdmann and Wong [2003]
Norwegian Continental Margin 456 m 270–300 m Mienert et al. [2005]
  • a Comparison of the mean present-day thickness of the GHSZ (LGHSZ in m), including all values >0 m, estimated according to the method described in 3 with recently published results obtained after calculating the depth of the Bottom Simulating Reflector (BSR in m) for different regions (adapted from Piñero et al. [2013]). Regions are depicted in Figure 6.

Using the predicted LGHSZ (Figure 3), the Quaternary sedimentation rates, and the POC concentration in surface sediments, the abundance and distribution of methane hydrates in marine sediments are calculated using equation 6. Figure 4 shows widespread accumulations of hydrate along the continental slopes. The shallower (<300 m water depth) and most of the deeper parts of the ocean do not contain any methane hydrate because the former features too warm bottom water temperatures and the latter is subject to too slow sedimentation as well as POC accumulation rates. In fact, the highest amounts of methane hydrates with values exceeding 260 kg C/m2 are found in areas with high amounts of organic carbon and high sedimentation rates, i.e., off Colombia, Uruguay, the eastern coast of Africa, and the Arabian Sea (Figure 4). The present-day global inventory of marine methane hydrates calculated for Quaternary sedimentation rates amounts to 1146 Gt of methane carbon (Figure 2). In contrast, using Holocene sedimentation rates would provide a minimum estimate of methane hydrate accumulation in marine sediments of just 3 Gt C.

Details are in the caption following the image
Present-day global distribution of methane hydrates (kg C/m2) in marine sediments calculated for Quaternary sedimentation rates. For the performed calculations, the sedimentation rates were increased over a 500 km wide range close to the coast line; at distance (>500 km) sedimentation rates have been obtained according to equation 8. Note that our sedimentation model explains the predicted hydrate occurrence around, e.g., Hawaii and the Galapagos Islands, although no hydrate is observed there. The dashed contour line marks the 300 m isobath.

3.2 Future Climate Conditions

Considering the ensemble mean trend of bottom water temperatures for future climate change (Figure 5), it is clearly visible that the most dominant changes will occur in the shelf regions, with shallow regions of the Arctic as well as off the coast of Australia and Nova Scotia, around Indonesia, and in the Yellow Sea rising by more than 3°C within the next century. However, the deep ocean appears to remain largely unaffected. Li et al. [2013] found that the deep ocean will respond to global warming with some delay on centennial to millennial time scales (refer to Figure 3 therein). These time scales are not considered in our 100 year approach. In the context of destabilizing methane hydrates, the changes in the deeper ocean play a minor role. The Labrador Sea is even predicted to cool by about 0.5° to 1.5°C at a depth below 1000 m because of a decreased northward heat transport related to a weakened Meridional Overturning Circulation (MOC) [Weaver et al., 2012; IPCC, 2013].

Details are in the caption following the image
Global distribution of the ensemble mean trend in bottom water temperatures (°C) over the next 100 years. The dashed contour line marks the 300 m isobath.

3.2.1 Steady State Analysis

Under global warming conditions, the thickness of the GHSZ will decrease along the global continental slopes by up to 60 m if steady state conditions are considered (Figure 6), i.e., the rate of heat diffusion is ignored and a steady state geotherm for the changed bottom water temperature is assumed. This analysis provides a clearer picture of what may happen to marine gas hydrates under global warming conditions but over predicts the total volume of possible dissociating hydrates. More realistic volume estimates are provided later in the transient analysis. Under steady state conditions, most changes occur in the depth range 300 to ∼1500 m; the deeper ocean and the shelf regions will probably not be subject to any considerable changes due to temperatures remaining in the hydrate or in the gas phase because of the prevailing temperature and pressure conditions.

Details are in the caption following the image
Global map of the predicted thickness of the gas hydrate stability zone (LGHSZ in m), calculated under steady state conditions for the ensemble mean trend, where the colored shadings indicate either a reduction (red) or an increase (blue) in the LGHSZ. The colored boxes mark the regions of the Cascadia Margin (CM), Gulf of Mexico (GOM), Blake Ridge (BR), Gulf Stream (GS), Norwegian Continental Margin (NOR), and Sea of Okhotsk (OS). The dashed contour line marks the 300 m isobath.

By comparing the bottom water temperature trend (Figure 5) and predicted GHSZ change (Figure 6), it becomes evident that most warming will occur on the continental shelves and the upper parts of the continental slopes. As a consequence, most hydrate melting occurs on the upper continental slope due to the retreat of the GHSZ (Figure 1, box A and Figure 6). The vertical reduction of the hydrate stability zone (box B in Figure 1) is of secondary importance within hundreds of years, due to the slow diffusion of heat, and only occurs on the lower parts of the continental slope. The deeper ocean is not significantly affected by seafloor warming and hydrate melting over the next century. In some areas, e.g., the Labrador Sea, the deeper ocean is even subject to a cooling trend. We do not, however, consider the (unlikely) formation of new hydrates over the next hundred years in those areas and for the following only the areas featuring a reduction of the LGHSZ under the future global warming scenario (red shadings in Figure 6) are taken into consideration.

In the steady state scenario the global volume of sediment within the hydrate stability zone will reduce by approximately 0.3%. In respect to the methane hydrate inventory, this results in a reduction of the global inventory by about 6.7 Gt C (−0.6%).

We chose seven example regions (Table 2 and Figure 7) with known hydrate occurrences for further detailed analysis [cf. Kvenvolden and Lorenson, 2001; Lüdmann and Wong, 2003; Riedel et al., 2006; Bohrmann and Torres, 2006; Heeschen et al., 2007; Hester and Brewer, 2009]. Assuming steady state, the greatest relative changes in the methane hydrate inventory will occur at the Blake Ridge region offshore South Carolina (−11.5%), followed by the Arctic (−1.4%). In contrast, the Cascadian Margin will feature the smallest relative changes which amount to a reduction of about 5 Mt C (<−0.1%). The Norwegian Continental Margin will also be subject to a rather weak decrease in the GHI of 2 Mt C. Note that although covering only ∼3% of the global seafloor, the Arctic Ocean accounts for ∼24% of the global hydrate inventory loss in the steady state analysis.

Table 2. Volume of Sediment Within the GHSZ and Methane Hydrate Inventories (Present-Day Estimates and Future Changes)a
Sediment Volume Change in Gas
Region Within GHSZ Gas Hydrate Inventoryd Hydrate Inventory
Arctic 3.80·1015 m3 116 Gt C −1.4%b, −0.12%c
Blake Ridge 0.18·1015 m3 1 Gt C −11.5%b, −0.99%c
Cascadia Margin 0.11·1015 m3 9 Gt C <−0.1%b, <-0.01%c
Sea of Okhotsk 0.91·1015 m3 39 Gt C −0.7%b, −0.05%c
Gulf of Mexico 0.54·1015 m3 43 Gt C <−0.9%b, −0.03%c
Norwegian Continental Margin 0.04·1015 m3 1 Gt C <−0.2%b, −0.04%c
Gulf Stream 0.12·1015 m3 5 Gt C −0.9%b, −0.05%c
Global 87·1015 m3 1146 Gt C −0.6%b, −0.03%c
  • a Present-day estimates of volume of sediment within the GHSZ (in m3) and methane hydrate inventories (GHI in Gt C) for the different regions depicted in Figure 6. The changes in the methane hydrate inventories under global warming (in %) are related to present-day climate conditions. Note that no active fluid flow is considered for the calculation of the GHI.
  • b Change in GHI under steady state conditions.
  • c Change in GHI within the next 100 years under transient state conditions.
  • d The GHI values have been rounded to the nearest whole number.
Details are in the caption following the image
The predicted change in the distribution of methane hydrates (kg C/m2) under steady state conditions for the ensemble mean trend for (a) the Arctic and the Norwegian Continental Margin, (b) the Cascadia Margin, Gulf of Mexico, Blake Ridge, and Gulf Stream, and (c) the Sea of Okhotsk. The solid contour lines mark the 300 m isobath.

Figure 8 shows seafloor warming over the next 100 years for the seven selected example sites in the context of the methane hydrate phase diagram. The open symbols indicate bottom water temperatures of the present-day climate and the filled symbols those of the global warming experiment. Progressive (and eventual complete) dissociation of hydrate occurs where the seafloor temperature rises across the phase boundary (cf. box A in Figure 1). At midlatitude regions this happens between 450 m (Sea of Okhotsk in Figure 8c) and 600 m (Blake Ridge in Figure 8b and Gulf of Mexico in Figure 8c). In the Arctic, this phase shift occurs already at shallower depths (∼380 m) due to the lower bottom water temperatures. Regions where the seafloor is under present-day climate within the hydrate stability field and where warming does not result in a phase shift at the seafloor, a progressive vertical reduction in the GHSZ thickness from below occurs (cf. box B in Figure 1) on time scales of hundreds to thousands of years. At all example sites the magnitude of warming diminishes with increasing water depth.

Details are in the caption following the image
Phase diagrams of methane hydrate as a function of pressure (here denoted by water depth in m) and temperature (°C) considering a constant salinity of 35 for (a) the Arctic and Norwegian Continental Margin, (b) the Blake Ridge and Cascadian Margin and for (c) the Sea of Okhotsk, Gulf Stream, and Gulf of Mexico. The open symbols indicate bottom water temperatures of the present-day climate (time mean 1988 to 2007) and the closed symbols those of the global warming experiment. The error bars mark the vertical resolution of the ORCA025 model configuration. Note the different scales of the x axes.

3.2.2 Transient State Analysis

A steady state will not be reached within the next century since heat penetration into marine sediments is relatively slow (see green dotted lines in boxes A and B in Figure 1) [Ruppel, 2011; Parmentier et al., 2013]. The more realistic calculation thus accounts for transient heat diffusion into the sediments (equation 3).

Figure 9a shows that within the next 100 years the global volume of sediment within the GHSZ will decrease on average by only 2.4·1013 m3 (∼ 0.03% of present-day volume). The here not considered consumption of latent heat during hydrate melting might reduce this value even further. Since the global warming scenario with KCM is limited to the next 100 years, it is assumed that bottom water temperatures will remain constant beyond that time frame and volume changes of the GHSZ were estimated for the next 200 to 500 years (see Figure 9a) by considering a longer time of heat diffusion into the sediment. However, this also means that we calculate a lower estimate, since it is unlikely that bottom water temperatures remain constant after the stabilization of the greenhouse gas concentrations [Park et al., 2009].

Details are in the caption following the image
(a) Global mean volume change in the GHSZ thickness (1013 m3) and (b) global mean change in the methane hydrate inventory (Mt C) as a function of time. The error bars indicate the standard deviations of the 22 ensemble members. The changes are related to present-day climate conditions (Figure 2 and Table 2).

Similar to the global values, the regional breakdown also features relative small reductions. Since comparisons of absolute changes make no sense for individual (somehow arbitrary defined) regions, we concentrate on the relative changes. As already indicated by the analysis of the bottom water temperature (Figure 8), largest relative changes are expected in the Arctic and at Blake Ridge offshore South Carolina. Accordingly, the amount of methane released from the seafloor varies among the investigation areas by 3 orders of magnitude. The largest relative changes will occur at Blake Ridge, where the inventory is on average going to be reduced by almost 6% (−54 ± 5 Mt C) after 500 years of heat penetration into the sediment owing to the strong bottom water temperature increase at the end of the 21st century (cf. Figure 5). The largest absolute release of methane gas will occur in the Arctic, with between 140 ± 10 Mt C and 410 ± 30 Mt C entering the water column of the Arctic Ocean within the next 100 to 500 years. In contrast, changes along the Norwegian Continental Margin are relatively small (∼0.9 ± 0.2 Mt C after 500 years), despite the strong expose to the poleward flowing continuation of the North Atlantic Current. In total, the seven chosen investigation areas account on average for about 52% of the global hydrate inventory change when considering a transient warming of the sediment column over 100 years, whereby the global inventory reduces by ∼355 ± 23.5 Mt C (Figure 9b).

4 Discussion

For the present-day climate we estimate that about 1146 Gt C are stored within marine methane hydrate deposits. This value fits fairly well into the range of previously published global estimates (cf. Figure 2), which were either based on field observations or on model studies. Our study improves previous model-based estimates in that we use a high-resolution ocean hindcast model to predict current bottom water temperature and salinity values, global input data sets compiled at the resolution of the ocean model, and a transfer function that is based on a 1-D reaction transport model, which takes the full complexity of the problem into account and has been successfully verified at various test sites [Wallmann et al., 2006, 2012].

It should be noted that direct observations might point to very different local hydrate inventories compared to the predicted values shown in Figure 4. This mismatch arises as local pathways for fluid and gas flow are not considered in our approach. Both fluid and gas flow have a strong effect on the accumulation of gas hydrates in marine sediments and could enhance the hydrate inventory [Zatsepina and Buffett, 1997; Xu and Ruppel, 1999; Buffett and Archer, 2004; Piñero et al., 2013]. Therefore, the value of 1146 Gt C should be considered as a minimum estimate based on the in situ microbial degradation of POC. In order to match specific local observations on hydrate occurrences with numerical models, detailed 3-D model-data integration is necessary. Such an approach is currently not feasible on a global scale due to limitations in both current models and available data sets.

In respect to future climate change, we have shown that the stability of methane hydrates as well as the total hydrate inventory is substantially affected by changing bottom water temperatures. Within the next century, hydrate dissociation will especially occur where the hydrate stability limit is shifted down the continental slope. Hydrate deposits that are located at greater depth will not be affected considerably within the next century owing to the long ventilation times of the deep ocean of up to 1000 years [Krey et al., 2009] and the slow rate of heat diffusion within marine sediments [Ruppel, 2011; Parmentier et al., 2013]. Only a small fraction of the global hydrate inventory is therefore affected by global warming over the next century. A possible effect on global climate thus only seems plausible on millennial and longer time scales [cf. Archer, 2007].

Nevertheless, warming bottom waters will destabilize hydrate deposits on the continental slopes at water depths of ∼300−700 m over the next century (see Figure 8), which will result in the release of noticeable amounts of methane into the ocean and potentially into the atmosphere. Especially, the deposits found in the Arctic, at relatively shallow water depth, could undergo rapid dissociation. Observed gas flares and seepage sites [Westbrook et al., 2009; Shakhova et al., 2010] suggest that this is already happening but the evidence remains inconclusive, and the observed venting may well be unrelated to global warming [Berndt et al., 2014]. What is clear, however, is that the absolute numbers of methane release will be low. Compared to the 354 Mt CH4/yr of present-day annual anthropogenic methane emissions [IPCC, 2013], the here predicted ∼473 Mt CH4 of globally released methane covering the full time span of 100 years in the transient scenario is small. This is even an unrealistic upper estimate, since most of the methane released from the seafloor is utilized by microbial processes. Although highly uncertain, it is possible that at least 50% of the dissolved methane is retained inside marine sediments due to anaerobic methane oxidation [Treude et al., 2003; Reeburgh, 2007; Knittel and Boetius, 2009]. Mau et al. [2007] recently investigated a massive seep field and found that only ∼1% of the methane leaving the seafloor reaches the atmosphere.

Compared to Biastoch et al. [2011], we note that the amount of methane and carbon inventories in the Arctic presented here is much smaller. Specific details in the calculation of the GHSZ thickness (higher-resolved ocean model, more global warming ensemble members, and a gradual temperature increase instead of a step function) are of minor importance, given the fact that the GHSZ volume north of 60°N is quite similar. Most important is certainly the inventory calculation. The constant mean hydrate pore filling of 2.4% to 6.1% in Biastoch et al. [2011] led to a much higher inventory of 9000 Gt C for the present climate (note that the inventory in the original version of Biastoch et al. [2011] was incorrectly given by 900 Gt C). The almost 2 orders of magnitude lower value of ∼116 Gt C calculated here demonstrates that the main uncertainty in model-based hydrate estimates is not the volume of sediments within the hydrate stability zone but the computation of the amount of hydrate stored within it. Here we use a transfer function instead of simple pore-filling estimates. Future models that further improve on this will help reducing the large uncertainty in current global estimates (cf. Figure 2). Since the amount of released methane is much lower compared to Biastoch et al. [2011], the large-scale impact on acidification and deoxygenation is also smaller. Regional refinement (which is beyond the scope of this study) is needed to identify the local imprint of released methane on the near-bottom chemistry.

Apart from the Arctic, hydrate deposits in relatively shallow waters (300–700 m) in the Gulf of Mexico or along the western North Atlantic Margin will, most likely, be affected by global warming as well [Reagan and Moridis, 2007; Phrampus and Hornbach, 2012; Skarke et al., 2014]. In addition, variations in the location of ocean currents and an additional warming on decadal time scales might lead to the local dissociation of marine gas hydrates. Phrampus and Hornbach [2012] found that the onset of methane hydrate destabilization can already be observed along the western North Atlantic Margin and is attributable to a northwestward shift in the Gulf Stream flow direction as warmer waters are introduced in regions, which have previously been exposed to colder bottom water currents. Both changes in the Gulf Stream flow direction and in ocean temperature could trigger the melting of methane hydrate deposits. It is likely that changing ocean currents, i.e., variations in the flow direction and/or temperature distribution, might have a significant impact on the stability of methane hydrates in specific regions of the world ocean.

5 Conclusion

To improve our understanding of methane hydrates and their fate under global warming, detailed estimates of hydrate inventories are required. Here we provide improved present-day estimates, both for the global-scale and specific regions in which hydrates were confirmed by observations.

Despite the significant ocean bottom warming over the next 100 years and the resulting potential phase shifts from methane hydrate to free gas at middepth along the continental margins, we find that the release of methane into the water column is very limited throughout the next century. On a global scale, it is negligible compared to the current anthropogenic releases of methane and other greenhouse gases. However, due to the long time scales of deep oceanic processes and the slow penetration of heat into the sediment, the decrease of methane hydrates will continue even beyond a potential recovery of the atmospheric warming.

On a regional scale, the relative decrease of methane hydrate is larger, with the Arctic and Blake Ridge region being most affected by climate change. In this regard, our maps could provide some guidance in which regions the interplay between ocean circulation changes, ocean bottom warming, sediment structure, and localization and quantification of methane hydrate deposits needs to be further investigated.

Apart from the physical and geological approaches, advances are also needed to better quantify the amount of released methane from the sediment and into the atmosphere. Uncertainties of methane consumption within the water column due to microbial consumption are too large to provide reliable estimates on the impact of decreasing methane hydrates onto the climate system.

Acknowledgments

This paper has benefited from the constructive comments and suggestions of the two anonymous reviewers, which greatly helped to improve our manuscript. We gratefully thank E. Behrens and W. Park for providing the data of ORCA025 and KCM, respectively. We further like to thank K. Wallmann for his helpful advice and recommendations. The model integrations have been performed at the North German Supercomputing Alliance (HLRN) and the computing center at Kiel University. All data can be requested from the corresponding author (Kerstin Kretschmer, [email protected]).