Estimates of future warming‐induced methane emissions from hydrate offshore west Svalbard for a range of climate models

Methane hydrate close to the hydrate stability limit in seafloor sediment could represent an important source of methane to the oceans and atmosphere as the oceans warm. We investigate the extent to which patterns of past and future ocean‐temperature fluctuations influence hydrate stability in a region offshore West Svalbard where active gas venting has been observed. We model the transient behavior of the gas hydrate stability zone at 400–500 m water depth (mwd) in response to past temperature changes inferred from historical measurements and proxy data and we model future changes predicted by seven climate models and two climate‐forcing scenarios (Representative Concentration Pathways RCPs 2.6 and 8.5). We show that over the past 2000 year, a combination of annual and decadal temperature fluctuations could have triggered multiple hydrate‐sourced methane emissions from seabed shallower than 400 mwd during episodes when the multidecadal average temperature was similar to that over the last century (∼2.6°C). These temperature fluctuations can explain current methane emissions at 400 mwd, but decades to centuries of ocean warming are required to generate emissions in water deeper than 420 m. In the venting area, future methane emissions are relatively insensitive to the choice of climate model and RCP scenario until 2050 year, but are more sensitive to the RCP scenario after 2050 year. By 2100 CE, we estimate an ocean uptake of 97–1050 TgC from marine Arctic hydrate‐sourced methane emissions, which is 0.06–0.67% of the ocean uptake from anthropogenic CO2 emissions for the period 1750–2011.


Introduction
Methane hydrate in marine sediments may contain ;500-2500 Gt [e.g., Piñero et al., 2013] of carbon, of which ;100-600 Gt may be stored in the Arctic [Archer et al., 2009]. Hydrates form at low temperature-high pressure conditions where the dissolved methane concentration in the pore water within the gas hydrate stability zone (GHSZ) is at or above saturation value [Sloan and Koh, 2007]. The most intense future climate warming is predicted in the Arctic [Hassol, 2004] and hydrates are most sensitive to ocean warming at high latitudes and in shallow water depths [e.g., Hunter et al., 2013]. Past warm periods have been linked to hydrate dissociation [e.g., Dickens, 2011] and various authors propose that hydrate may dissociate again in the future in response to plausible future temperature scenarios [MacDonald, 1990;Nisbet, 1989;Reagan and Moridis, 2008;Reagan et al., 2011;Biastoch et al., 2011;Mar ın-Moreno et al., 2013].
Perhaps the best-documented example of an apparent response of Arctic hydrate to ongoing ocean warming is on the west Svalbard continental margin. Here, widespread methane venting has been observed near the landward limit of hydrate stability in 350-400 m water depth [Westbrook et al., 2009;Sahling et al., 2014]; hydrate-related bottom-simulating reflectors are widespread [Sarkar et al., 2012]; and hydrate has been sampled from beneath the seabed in pockmarks in water depths of 900-1200 m Panieri et al., 2014]. Thatcher et al. [2013] showed that the observed methane venting could be attributed to warming-induced hydrate dissociation over the last century. However, extensive methane venting also occurs further landward [Westbrook et al., 2009;Rajan et al., 2012;Sahling et al., 2014]. In addition, based on the ages of carbonate samples from the 350 to 400 mwd venting area and their isotopic composition, Berndt et al. [2014] showed that methane seepage has been active there for more than 500 year. From their modeling of the response of the subseabed temperature field to seabed temperature variations that were measured continuously over a period of nearly 2 years, Berndt et al. [2014] also suggested that the observed present-day emissions may be controlled by seasonal changes in temperature. However, they did not model To assess the range of possible future ocean temperature behavior and its impact on seafloor methane emissions west of Svalbard during the 21st century, we use seven different climate models for which outputs are available over this period. We focus on water depths of 400-500 m, which is the depth range of potential gas hydrate dissociation [Giustiniani et al., 2013;Mar ın-Moreno et al., 2013] during the 21st century. The temperature increase predicted by the climate models in deeper water is not sufficient to cause hydrate dissociation [Mar ın-Moreno et al., 2013]. Our approach differs from that of Hunter et al. [2013], who model hydrate dissociation globally in response to ocean temperature changes predicted by a range of climate models, in that we focus on a single geographical area with known emission of methane, has been well characterized geophysically which means that we make fewer assumptions about the physical parameters involved. We account for the important effects of latent heat, and we track the consequences of such dissociation in terms of seabed methane emissions. We also assess the influence of decadal and seasonal fluctuations on the past and future response of marine hydrate-bearing sediments offshore west Svalbard, and thus constrain the associated future Arctic methane emissions.

Processes Modeled
We ran one-dimensional (1-D) simulations for water depths of 400, 420, 450, and 500 m using the TOUGH 1 HYDRATE (T 1 H) code [Moridis et al., 2012]. We adopted a 1-D approach because: (1) a twodimensional (2-D) analysis of our study area [Reagan et al., 2011] suggested that for a slope of ;38, which is greater than the observed slope of 1-1.58 between 350 and 550 mwd [Thatcher et al., 2013], the 1-D and 2-D models generate similar results for low values of hydrate saturation, such as the 5% that we use here (Table 1); (2) an aquifer with a 18 slope generates a horizontal component of the head gradient of only 2% of the vertical component [Thatcher et al., 2013], so that the horizontal component of fluid flux is likely to be negligible; and (3) the run times for 2-D models are too long for the long time series and large number of temperature scenarios required in this analysis. Each 1-D model represents the response of the geological system beneath an area of about 1 km 2 and, at its lower limit, simulates the net behavior of the system feeding an individual gas flare, which occupies an area of about 15,000 m 2 [Thatcher et al., 2013].
T 1 H is a thermohydraulic code for the simulation of the behavior of gas-hydrate-bearing geologic systems under nonisothermal conditions. Here we describe the most important physical processes modeled; more detailed information is available from the online T 1 H manual (http://esd.lbl.gov/files/research/projects/ tough/documentation/ TplusH_Manual_v1. pdf). We imposed equilibrium conditions for hydrate formation and dissociation and considered three possible mass components (water, methane, and salt) and heat. These mass components were partitioned between four possible phases: hydrate (for water and methane), aqueous (for water, methane, and salt), gas (for water and methane), and ice (for water only). Heat exchanges due to conduction, convection, hydrate formation and dissociation, and methane and salt dissolution were modeled. Water and gas flows driven by pressure changes were modeled using Darcy's law. Darcy's law and a Fick's type law were adopted for the advective and molecular diffusive transport, respectively, of methane and salt within the aqueous phase. Calculated seabed methane emissions include contributions from both methane bubble flow and dissolved methane. When methane bubble flow occurs, its Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 flux is 2-3 orders of magnitude larger than the dissolved methane flux. Over the timespan of our runs (2100 year), the dissolved methane contribution to the total methane released to the ocean is more than 2 orders of magnitude smaller than that from methane bubble flow.
Following previous studies [e.g., Reagan et al., 2008Reagan et al., , 2011Thatcher et al., 2013], methane and water relative permeabilities were computed according to a modified version of Stone's [1970] first three-phase relative permeability method and the capillary pressure was calculated using van Genuchten's [1980] law (Table 1). Both models were adjusted using the Evolving Porous Medium (EPM) model #2 [Moridis et al., 2012].  / 0 is the porosiy at surface conditions; k i and k io are the initial intrinsic permeability and initial intrinsic permeability at surface conditions; n k is a fitting parameter; k rA and k rG are the relative permeabilities for aqueous and gas phases; S A , S G , and S H are the saturations for aqueous, gas, and hydrate phases; S irA and S irG are the irreducible aqueous and gas saturations; S mxA is the maximum water saturation; P is the pore pressure; P cap is the capillary pressure; P max is the maximum value of capillary pressure; P 0 is the capillary entry pressure; T is the temperature; n, n G and k are the fitting parameters.

10.1002/2015GC005737
The relative permeability model was adjusted for changes in the saturation of solid phases (ice or hydrate) occupying the pore space. The capillary pressure model was adjusted for both changes in porosity, resulting from changes in pressure and in saturation of solid phases in the pores, and changes in intrinsic and relative permeabilities, resulting from changes in porosity (Table 1). In our runs, changes in porosity due to changes in pressure (Table 1) are almost negligible.

Model Parameters
The physical properties of the gas and hydrate system and constraints derived from seismic data are summarized in Table 1. The 1-D models have a total height of 1.1 km and a variable cell height (Table 2). We applied the past and future temperatures as a top boundary condition, changing them annually or every 6 months, and initialized the models with seabed temperatures at 1 CE [Mar ın- Moreno et al., 2013], hydrostatic pressure, constant heat flow in the entire column, and 7 m of hydrate-free sediment at the top of the model [Thatcher et al., 2013] to simulate the presence of a sulfate reduction zone (SRZ). In marine anoxic sediments, the SRZ is the zone below the sediment-water interface, within which the interstitial sulfate concentration decreases with depth, by microbial sulfate depletion, and aqueous methane is consumed by anaerobic oxidation [Boetius and Wenzh€ ofer, 2013;Borowski et al., 1996]; we did not attempt to model this process. We assumed an irreducible gas saturation (defined as the concentration of gas above which gas flows) of 2%, consistent with other modeling studies [e.g., Liu and Flemings, 2007;Thatcher et al., 2013]. We adopted a two-layer model comprising glaciogenic sediments overlying hemipelagic marine sediments [Thatcher et al., 2013] with initial porosities of 0.3 and 0.5, respectively. We estimated the initial hydrate saturation in the GHSZ and the free-gas saturation below the GHSZ to be 5% and 3-4%, respectively, from the analysis of P and S wave velocities in our study area [Chabert et al., 2011] (Table 1), using the approach applied by Westbrook et al. [2008] to similar seismic data from farther downslope, in water depths of 1400 m. For both sediment types, we assumed a thermal conductivity for fully water-saturated sediment of 1.4 W m 21 K 21 [Mar ın- Moreno et al., 2013], and an initial intrinsic permeability of 10 213 m 2 [Thatcher et al., 2013]. At the bottom of the model, there is a constant heat flow equal to heat flow used as initial condition, and a constant methane flow that approximately matches the rate of buoyancy-driven methane flow from below the GHSZ to the GHSZ. A discussion on intrinsic permeability is presented in the next section. Other details of our model assumptions and parameter uncertainties (hydrate concentration, heat flow, and thermal conductivity) are as given by Mar ın-Moreno et al. [2013].

Seabed Temperatures
The method used to construct the seabed temperature series for the period 1-2005 CE is described in Mar ın-Moreno et al. 's [2013] electronic material. We focus on the period 2005-2100 CE, for which the thermal effects of a range of climate models can be explored. In addition to the CCSM4 and HADGEM2 models used by Mar ın-Moreno et al. [2013], the models we used were: the NOOA Geophysical Fluid Dynamics Laboratory model (NOAA GFDL) [Gordon and Stern, 1982]; the Institute Pierre-Simon Laplace model (IPSL) [Marti et al., 2010]; the Max-Plank-Institute f€ ur Meteorologie model (MPI) [Marsland et al., 2003]; the Meteorological Research Institute model (MRI) [Yukimoto et al., 2001]; and the Norwegian Climate Centre model (NCC Nor-ESM1) [Bentsen et al., 2012;Iversen et al., 2012]. These models are leading, well-documented climate models that have been used widely in high-profile published model intercomparison studies [e.g., Stroeve et al., M-1 is the default mesh. At depths deeper than 100 m the mesh size increases by a factor of 1.1 with respect to the previous depth z-interval.
c Model runs using a processor of 2 GHz Intel Core i7 and memory of 4 GB 1333 MHz DDR3.
Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 2012]. The mean annual seabed temperatures given by each of the seven global climate models and two scenarios were interpolated to our study location.
Some of the climate models predict present-day temperatures poorly for our study area (supporting information Text S1, Table S1), so we used the climate model output to predict the future temperature changes rather than using the actual temperature values. To avoid an abrupt present-day step in temperature, ideally we would subtract from future temperatures the mean difference in temperature between historical observations and model predictions over a period of overlap. However, CTD measurements in the study area were compiled for the period 1975-2008 [Westbrook et al., 2009] and our model-derived temperature series starts at 2005, so there are only 4 years of overlap, which is insufficient to estimate the mean difference robustly. Therefore, instead, we corrected the model-derived series by the difference between the model temperature at 2005 and the mean observed temperature over the period 1975-2005 (supporting information Text S1, Table S1). A further limitation of our approach is that the CTD data are limited to the commonly ice-free period from May to October [Westbrook et al., 2009], whereas our model-derived series are mean annual seabed temperatures. For five of the climate models, the corrections were less than 618C (supporting information Text S1, Table S1), which is similar to annual model temperature variations around the tie point in 2005 (Figures 1a, 1c, 1e, and 1g), so the temperatures predicted by these models are similar to those measured in our study area and the corrections have no significance. In contrast, adjustments of more than 38C were required for the HadGEM2 and MPI models (supporting information Text S1, Table S1). For these two models, the uncorrected temperatures greatly overestimate those measured.

Seabed Temperature Variations
In our study area, predicted temperatures over the 21st century are significantly more sensitive to RCP scenario than to the climate model. The temperature variation predicted for the 21st century can be approximated well by a linear increase for scenario RCP 2.6, and by a quadratic increase for RCP 8.5, as indicated by the norms of the residuals between the temperatures from the models and those from the regression curves (Table 3 and supporting information Text S1, Tables S2, and S3). A linear regression of the mean model temperatures for scenario RCP 2.6 and for all water depths gives an increase of 0.008 6 0.0038C yr 21 . For RCP 8.5, the corresponding linear increase is 0.032 6 0.0058C yr 21 at 400 mwd and 0.031 6 0.0058C yr 21 at 500 mwd (Table 3 and supporting information Text S1, Tables S2, and S3). For all climate models and water depths, temperatures over the first quarter of the century are similar and independent of the scenario RCP used (Table 3 and Figures 1a, 1c, 1e, and 1g). The temperatures for the two scenarios start to increase at different rates at about 2050 year (Figures 1a, 1c, 1e, and 1g), and in the last quarter of the century they differ by about 1.7-1.88C (Table 3).

Sensitivity to Climate Model
At 400 mwd, the future response of the hydrate system varies little for five of the models, irrespective of RCP scenario, but shows significant differences for climate models IPSL and MPI (Figures 1b, 1f Figure 2) will disappear completely and that therefore the uncertainty in its future response to ocean warming is small. Here, for the climate models that give a subseabed temperature profile (supporting information Text S1, Figures S3, and S4) in which the temperatures shallower than ;5 m are generally a little lower than those needed to produce dissociation, perturbations in temperature move the system in and out the hydrate stability field producing seabed methane pulses. The results for the MPI climate model under RCP 2.6 show the most episodic seabed methaneemission behavior of all our models (Figure 1b), because for this climate model the long-term average seabed temperature is slightly lower than the temperature of dissociation throughout the century (Figure 3).
To test whether some of the observed peaks in methane emissions may be due to numerical instabilities derived from the relatively coarse mesh (M-1) used in the T 1 H models, we reran the T 1 H model corresponding to the MPI climate model at 400 mwd with two finer meshes (M-2 and M-3; Table 2 and Figure 3). The mesh was refined mainly within the top 20 m where hydrate is currently stable. The time step for the three models was varied automatically within the range 10-10 6 s depending on the stability of the convergence procedure. While there are differences between the resulting model outputs, the principal features Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 of the predicted variation in methane flow with time are well correlated ( Figure 3). The default mesh size chosen generates results that yield a normalized root mean square (NRMS) difference (equation (1)) in predicted methane fluxes during the period 2000-2100 of 10% and 14% from M-2 and M-3, respectively.

Sensitivity to Intrinsic Permeability
Our parameters represent the average properties of the system beneath an area of about a 1 km 2 if the minimum increment of water depth is 20 m and slope of the seabed is 1-1.58. Seismic data [Sarkar et al., 2012] show that the system is very heterogeneous and incorporates fractures. The pattern of occurrence of gas seeps shows that gas migration is focused along pathways of locally higher permeability, almost certainly provided by fractures. Consequently, the effective permeability of this heterogeneous system will be higher than that of pristine samples of the lithologies present. Our intrinsic permeability value of 10 213 m 2 is 2-4 orders of magnitude greater than that for hemipelagic sediments. It is unlikely that very low permeability can be maintained in the hemipelagic sediments when gas is being produced rapidly from hydrate dissociation, because, for an intrinsic permeability of about 10 216 m 2 , the pore pressure exceeds the lithostatic load only a few years after the dissociation of hydrate commences [Thatcher et al., 2013]. Consequently, fractures are likely to form and fluid flow to become dominated by the greater permeability of the fractures, rather than the original intergranular permeability of the sediments in which they occur. Similarly, Smith et al.
[2014] in modeling a methane vent system in the Ursa Basin (Gulf of Mexico) found that they needed to use an intrinsic permeability of 10 212 m 2 , up to 3 orders of magnitude higher than the values of permeability measured from samples of the lithologies present. Accordingly, the value of permeability we used implicitly embodies the assumption that shallow fractures increase the permeability of the system. Our modeling approach assumes that the separation of fractures is small compared with the area of the seabed represented by each 1-D model, and so the fracture permeability can be considered to contribute to the bulk permeability used for each cell of the model [cf. Freeze and Cherry, 1979;Hiscock, 2009].
We assessed the influence of intrinsic permeability in the time for hydrate-sourced methane to reach the seabed and in the magnitude of seabed methane flow. We also ran the T 1 H models corresponding to the HadGEM2 climate model for both RCP scenarios at 400 and 420 mwd using an intrinsic permeability of 10 215 m 2 (Figure 4). Reducing the intrinsic permeability 2 orders of magnitude delays the onset of methane emissions at the seabed between two to four decades. However, the maximum magnitude of seabed methane flow remains similar to that calculated with an intrinsic permeability of 10 213 m 2 . Both results are in agreement with those of Thatcher et al. [2013]. The onset of seabed methane emissions at 400 mwd is mainly controlled by the intrinsic permeability (Figure 4a) because, as explained above, at that depth the system response is similar irrespective of the rate of future seabed temperature increase. Figure 4a supports Thatcher et al.'s [2013] result that to account for emission of gas from the seabed at ;400 mwd at the  Geochemistry, Geophysics, Geosystems

10.1002/2015GC005737
present-day in response to measured warming in the late 20th and early 21st centuries, the models need an intrinsic permeability higher than the typical permeability of the kinds of hemipelagic sediments present. At 420 mwd, a similar time delay occurs both when decreasing the intrinsic permeability 2 orders of magnitude or when decreasing the average rate of future seabed temperature increase by about 4 times, from the 0.0318C yr 21 of RCP 8.5 to the 0.0078C yr 21 of RCP 2.6 ( Figure 4b).   Table 2) at 400 m water depth for the 21st century using temperatures from the MPI climate model RCP 2.6.

Sensitivity to Short-Period Temperature Fluctuations
At 420 mwd, interannual temperature perturbations from the climate models, at about 60.58C, are similar to those at 400 mwd, but significant pulses of methane emission are not observed (Figures 1d and 1h).
Here such temperature perturbations are not large enough to cause hydrate dissociation and hence, for water depths deeper than ;420 m, a long-period of temperature increase (decades to centuries) due to ocean warming is required to cause methane emissions. Temperature variations in sediments containing hydrate and/or gas, driven by seabed temperature fluctuations, depend on porosity, hydrate and gas saturations, and thermal conductivity. These parameters are similar at 400 and 420 mwd and, at the same depth below seafloor and similar temperature, the higher pore pressure of the model at 420 mwd makes the hydrate stable (pressure at 420 mwd above that at the aqueous-gas-hydrate phase boundary for a given temperature) and explains these different can cause hydrate dissociation at 400 mwd but not at 420 mwd. If the supply of gas to the base of the GHSZ has not changed significantly in the past two millennia, our models predict emissions at 400 mwd during periods when temperatures were similar to the mean over the last century, which is ;2.68C (supporting information in Mar ın-Moreno et al. [2013]). The occurrence of such emissions is also sensitive to the assumed thickness of the hydrate-free zone [Thatcher et al., 2013].
In the models of Mar ın-Moreno et al. [2013], 50 year (1-1900 CE) and 15 year (1900-1950 CE) running means of annual temperatures were used, so seasonal to decadal temperature fluctuations were not considered. For periods after 1950 CE, decadal temperature fluctuations were included but seasonal fluctuations were not, because temperature was sampled at 1 year intervals. To explore the models' sensitivity to such fluctuations, they were run again for the period 1-2300 CE, at 400 and 420 mwd, using temperatures predicted by HadGEM2 RCP 2.6, with the addition of random temperature variations with periods of years (from 1 to 2300 CE) to decades (from 1 to 1950 CE). Seasonal variations were represented by a square wave of 21.0 to 1.08C amplitude simulating warm (summer-autumn) to cold (winter-spring) seasons, to which was added a random series with values between 21.0 and 1.0 (as illustrated in the inset in Figure 5a). We represented random decadal variations using a sinusoidal wave with an amplitude of 0.58C peak to peak, to which were added random values between 20.15 and 0.158C in phase with the sinusoidal wave values, resulting in Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 decadal variations with a maximum amplitude of 0.88C peak to peak as shown by Biastoch et al. [2011] in our study area. Note that such random series also introduce longer-period variations [Hasselmann, 1976]. The random series (annual and decadal) give different long-timescale behaviors with respect to the times at which maxima and minima of different periods occur. However, the amplitudes of temperature variation for different periods in each random series are very similar to each other and similar to the natural temperature variations in the latter part of the 20th century (compare Figure 6b with the temperature fluctuations shown in Biastoch et al. [2011, Figure 2b]. Here we discuss the results from one combination of series, as an  (Table 1). (c) Perturbed temperatures used in plots (a, red solid line) and (b, black solid line) and difference between associated seabed methane flows.
Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 example of the possible response of the system, but not a product of actual decadal and annual fluctuations, which are unknown for the period prior to 1950. The models were run for initial hydrate-free thicknesses of 1 and 7 m to test also for sensitivity to this parameter. As methane rises into the hydrate-free zone, the models reduce its thickness (Figures 2a and 2b).
At 400 and 420 mwd, the first methane emissions at the seabed are later for an initial hydrate-free thickness of 7 m than for 1 m thickness ( Figure 5 and supporting information Text S1, Figure S2). The delay is determined by the time taken for the interstitial water in the hydrate-free zone to saturate in methane and the distance the hydrate-sourced methane needs to travel to reach the seabed.
Annual and decadal temperature fluctuations at ;400 mwd may have contributed to gas hydrate dissociation and methane emissions in the past (Figures 5a and 5b). Annual to decadal temperature fluctuations only have a significant effect to a depth of ;5 m below seabed, but the long-period components (few decades to centuries) of these fluctuations (Figure 6b) in the quasiannual series (annual plus random) penetrate more deeply (Figure 6a). These long-wavelength fluctuations lead to the migration of the geotherm toward the hydrate phase boundary in the bottom part of the GHSZ (initially located at ;37 m below seafloor; Figures 2a and 2b). Once the system is at the phase boundary, they lead to hydrate dissociation by supplying the latent heat of fusion, and release of gas at concentrations above the irreducible gas saturation of 2%. However, the upward methane flow is slow, because hydrate reforms in the next cooling interval. This cyclic process may have governed the system until about 1900 CE, except during long episodes of temperatures similar to the mean over the last century of ;2.68C (Figure 5a). During such episodes, the warmer temperatures would have increased the overall rate of dissociation by keeping the temperature close to that of the hydrate phase boundary and making the hydrate system beneath the seabed with a water depth of around 400 m vulnerable to the effects of the longer-period components of the short-period temperature-fluctuation series (Figures 5a and 5b). These fluctuations produce multiple short periods (several years to a few tens of years) of emissions. This type of seabed emission activity can occur in any water depth if the longterm variation in temperature has brought hydrate close enough to the phase boundary. Although the resulting periods of temperature above the phase boundary are relatively short, they are long enough to produce emissions of methane lasting several years or even a few tens of years. These periods of higher temperature are compensated by similar periods of temperature lower than the long-term temperature trend, during which no methane emission occurs. The periods of no emission can be much longer than the periods of emission, but the important effect of long-term variations in temperature is that some methane emission occurs. Consequently, there will be evidence of emissions, in the form of authigenic carbonates [Berndt et al., 2014], for example, even though the modeling of past emissions from historical records and proxy data predicts no emission when the resolution of the thermal time series is too poor to resolve seasonal variation. Similarly, the prediction of future emissions is likely to provide an underestimate if the climate model does not predict seasonal variation adequately.
At present and over the 21st century, the long-term temperature at 400 mwd is predicted to be above that of the phase boundary. The gas hydrate system at 400 mwd is more influenced by decadal temperature fluctuations than by seasonal fluctuations. Therefore, at 400 mwd, the seasonal fluctuations result in only small increases in methane emissions (Figure 7). In our models, emissions before ;1900 CE do not result in significant depletion of the hydrate reservoir and its current thickness of ;18.5 m (Figures 2a and 2b) remains consistent with seismic observations (Table 1). Note that during the 2300 year of simulation, methane entering the base of the model (Table 1) does not rise far enough to contribute directly to the saturation of the hydrate reservoir, and it is the inflow of free methane initially present in the sediments below the GHSZ (Table 1) that prevents a significant depletion of the hydrate reservoir due to emissions before 1900 CE. The methane emissions over the past 500 year or so at 390 mwd [Berndt et al., 2014] could have been caused in large part by the effect of short-period temperature variation. Furthermore, particularly in the more recent past, the disappearance of both the GHSZ and hydrate, because of the relatively shallow depth of the seabed and the changing temperature of the seawater, may have allowed methane migrating up the continental slope (a process not represented in our 1-D models) to emerge from the seabed where the temperature is too high for hydrate to be stable. During cooler episodes, the GHSZ will have been reinstated and more hydrate formed. In sufficiently shallow water, however, the GHSZ has been absent for at least two millennia and methane emission has been perennial.

Predictions of Future Gas Emission
The time taken for methane gas to reach the seabed in our models is mainly controlled by the intrinsic permeability of the sediments and by the rate of change of the seabed temperature. The higher the rate of temperature increase, the sooner the methane emissions occur (Figure 8a). Seabed methane emissions at 420 mwd start after 2080 CE for RCP 2.6 and between 2060 and 2085 CE for RCP 8.5 (Figure 8a). A difference in the rate of temperature increase between 1.2 3 10 24 8C yr 21 (RCP 2.6) and 3.5 3 10 24 8C yr 21 (RCP 8.5), results in a 1 year difference in timing of the first seabed methane emission (Figure 8a), and a centuryaveraged difference in seabed methane emissions of 0.22 mol yr 21 m 22 (Figure 8b). At 450 mwd, seabed methane emissions occur only if using RCP 8.5 and these start between 2074 and 2096 CE for five of the climate models and not at all for MPI and NorESM1 (Figure 8b). At 500 mwd, none of the models predict gas emissions (supporting information Text S1, Figures S1d, and S1h).
Seabed temperatures when the first methane emission occurs vary little between climate models, but differ by about 0.5-0.78C between RCP scenarios (Figure 8a). The ratio between the rate of hydrate dissociation Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 and the rate of pressure dissipation is greater for RCP 8.5 than for RCP 2.6. Therefore, during dissociation, RCP 8.5 models require a greater increment in temperature to compensate for the excess pore pressure and move the system toward the hydrate phase boundary. Generally, greater rates of temperature increase lead to higher dissociation temperatures.
At 400 mwd, all climate models result in a similar maximum rate of methane emissions of 75-95 mol yr 21 m 22 for both RCP scenarios. For scenario RCP 2.6, we obtain a maximum rate of 10-20 mol yr 21 m 22 at 420 mwd and for RCP 8.5, the rate is 25-35 mol yr 21 m 22 at 420-450 mwd (Figures 1b, 1d, 1f, and 1h and supporting information Text S1, Figure S1b, and S1f). The magnitude of methane emissions is limited by the rate of supply of the heat required for dissociation [Thatcher et al., 2013], and this rate varies little with water depth. During the 21st century, at 400 mwd, the subsurface is mainly out of the GHSZ, and methane emissions come from    Figure  1h), which is the model that predicts the earliest seabed methane emissions, does methane from dissociation at the base of the GHSZ at 420 mwd start to contribute to the total methane outflow before 2100 CE.
The totals of hydrate-related methane emissions from the study area are estimated by using 21st century average methane fluxes calculated from our models at 400, 420, 450, and 500 mwd and applying a linear interpolation of the fluxes between these water depths. For the seabed depth range of our study, methane emissions are a monotonic function of hydrostatic pressure. We consider a seabed slope of 1.58 [Mar ın- Moreno et al., 2013], and assume emissions extend 11 km along the margin (Area 3 in Sahling et al. [2014]), margin length that comprises the same geological setting as our study area. Over this century, using scenario RCP 2.6 the active seabed area of methane emissions is then ;11.8 km 2 , from 400 to 430 mwd, releasing 0.4-1.5 Gg yr 21 of methane (2.4-8.3 3 10 3 mol yr 21 ) per meter along the margin. For RCP 8.5, methane emissions occupy a seabed area of ;31.4 km 2 , from 400 to 480 mwd, releasing 1.7-4.5 Gg yr 21 (9.6-25.7 3 10 3 mol yr 21 ) per meter along the margin (Figure 9). Sahling et al. To illustrate the potential significance of Arctic hydrate dissociation, we have extrapolated our 21st century hydrate-sourced methane emissions over the entire Eurasian Margin (738N-858N; 08W-1608W, going eastward). These extrapolations have many limitations: current ocean temperatures are colder further east (so hydrate is stable at shallower ocean depths); ocean temperature changes are likely to vary along the margin; more gentle margin slopes in shallow waters may result in a larger potential area of gas hydrate dissociation; and porosity, permeability, thermal conductivity, hydrate saturation and distribution, and heat flow [e.g., Crane et al., 1991] in the sediments are likely to be heterogeneous along the margin. However, the predicted emissions do allow a comparison with emissions from other global natural [e.g., Dlugokencky et al., 2011;McGuire et al., 2012] and anthropogenic [Rhein et al., 2013] methane sources and with other published estimates of hydrate-sourced methane flux [Westbrook et al., 2009]. If we ignore the above limitations, the future potential dissociation area is ;38,878 km 2 (seabed range from 400 to 430 mwd) for RCP 2.6 and ;98,475 km 2 (seabed range from 400 to 480 mwd) [Jakobsson et al., 2008] for RCP 8.5. Therefore, the scaling factors between the plume area and the Eurasian Margin are 3295 (RCP 2.6) and 3136 (RCP 8.5), and the corresponding methane emissions are 1.3-4.9 Tg yr 21 (0.97-3.65 TgC yr 21 ) and 5.3-14.1 Tg yr 21 (3.95-10.5 TgC yr 21 ), respectively (Figure 9). These numbers are slightly lower than the 20 Tg yr 21 estimated using a similar approach by Westbrook et al. [2009], and more than an order of magnitude lower than the 162 Tg yr 21 estimated for the entire Arctic by Biastoch et al. [2011] using ocean temperatures from a coupled climate model. Hydrate-sourced methane emissions are unlikely to reach the atmosphere directly [e.g., Westbrook et al., 2009], and are likely to increase ocean acidification [e.g., Biastoch et al., 2011]. Our estimated methane emissions to the ocean are 2-3 orders of magnitude smaller than the anthropogenic CO 2 average ocean uptake of 1.0-3.2 PgC yr 21 (from fossil fuel combustion, cement production, and deforestation and other land use change) [Rhein et al., 2013]. Neglecting the contribution of methane transfer to the atmosphere by equilibration, by 2100 CE the total carbon taken up by the ocean from marine Arctic hydrate-sourced methane emissions may reach 97-1050 Tg. This carbon uptake is only 0.06-0.67% of the total 155 Pg of carbon taken up by the ocean from anthropogenic CO 2 emissions over the period 1750-2011 [Rhein et al., 2013].

Conclusions
From our modeling of past and future methane emissions resulting from hydrate dissociation beneath the continental slope west of Svalbard, we conclude the following: Geochemistry, Geophysics, Geosystems 10.1002/2015GC005737 1. Over the past 2000 year and during periods in which the average seabed temperature over a few decades was similar to the mean over the last century (2.68C), the gas hydrate system at water depths of 400 m and shallower would have been vulnerable to the effects of decadal temperature fluctuations and the longer-wavelength components arising from seasonal fluctuations. These fluctuations would have produce methane emissions over periods of several years to a few tens of years, even when the long-term temperature trends, by themselves, were insufficient to generate methane emissions.
2. The shorter-period fluctuations in seabed temperature could explain the presence of authigenic carbonates that are more than 500 years old at the seabed in water depths of around 400 m. In shallower water, this process of stimulation of gas emission becomes increasingly important, until long-term temperature change causes the sediments beneath the seabed to be permanently outside the GHSZ for hundreds of years, allowing gas flowing through the sediments toward the seabed to enter directly into the ocean.
3. The total carbon taken up by the ocean from marine Arctic hydrate-sourced methane emissions over the 21st century may reach between 97 and 1050 TgC, which is 0.06-0.67% of the ocean uptake of 155 PgC from anthropogenic CO 2 emissions for the period 1750-2011 [Rhein et al., 2013].
4. The pattern of occurrence of gas seeps shows that gas migration is almost certainly provided by fractures. Our models require an intrinsic permeability of about 10 213 m 2 to produce the present-day methane emissions from the seabed at ;400 mwd in response to measured warming in the late 20th and early 21st century. This permeability value is higher than the typical permeability of most of the kinds of hemipelagic sediments present and embodies the assumption that the real system is heterogeneous and incorporates closely spaced fractures compared with the area of the seabed represented by each 1-D model.

5.
Predicted methane emissions are mostly insensitive to the choice of climate model and RCP scenario (with the exception of IPSL) over the first half of the 21st century. Therefore, the uncertainty in our estimated methane emissions over that period is small. In contrast, predicted emissions for the period 2050-2100 are significantly sensitive to the choice of both climate model and RCP scenario.