Climate of the Past Inferred gas hydrate and permafrost stability history models linked to climate change in the Beaufort-Mackenzie Basin , Arctic Canada

Atmospheric methane from episodic gas hydrate (GH) destabilization, the “clathrate gun” hypothesis, is proposed to affect past climates, possibly since the Phanerozoic began or earlier. In the terrestrial Beaufort-Mackenzie Basin (BMB), GHs occur commonly below thick ice-bearing permafrost (IBP), but they are rare within it. Two end-member GH models, where gas is either trapped conventionally (Case 1) or where it is trapped dynamically by GH formation (Case 2), were simulated using profile (1-D) models and a 14 Myr ground surface temperature (GST) history based on marine isotopic data, adjusted to the study setting, constrained by deep heat flow, sedimentary succession conductivity, and observed IBP and Type I GH contacts in Mallik wells. Models consider latent heat effects throughout the IBP and GH intervals. Case 1 GHs formed at ∼0.9 km depth only∼1 Myr ago by in situ transformation of conventionally trapped natural gas. Case 2 GHs begin to form at ∼290–300 m∼6 Myr ago in the absence of lithological migration barriers. During glacial intervals Case 2 GH layers expand both downward and upward as the permafrost grows downward through and intercalated with GHs. The distinctive model results suggest that most BMB GHs resemble Case 1 models, based on the observed distinct and separate occurrences of GHs and IBP and the lack of observed GH intercalations in IBP. Case 2 GHs formed>255 m, below a persistent ice-filled permafrost layer that is as effective a seal to upward methane migration as are Case 1 lithological seals. All models respond to GST variations, but in a delayed and muted manner such that GH layers continue to grow even as the GST begins to increase. The models show that the GH stability zone history is buffered strongly by IBP during the interglacials. Thick IBP and GHs could have persisted since ∼1.0 Myr ago and∼4.0 Myr ago for Cases 1 and 2, respectively. Offshore BMB IBP and GHs formed terrestrially during Pleistocene sea level low stands. Where IBP is sufficiently thick, both IBP and GHs persist even where inundated by a Holocene sea level rise and both are also expected to persist into the next glacial even if atmospheric CO 2 doubles. We do not address the “clathrate gun” hypothesis directly, but our modls show that sub-IBP GHs respond to, rather than cause GST changes, due to both how GST changes propagates with depth and latent heat effects. Models show that many thick GH accumulations are prevented from contributing methane to the atmosphere, because they are almost certainly trapped below either ice-filled IBP or lithological barriers. Where permafrost is sufficiently thick, combinations of geological structure, thermal processes and material properties make sub-IBP GHs unlikely sources for significant atmospheric methane fluxes. Our sub-IBP GH model histories suggest that similar models applied to other GH settings could improve the understanding of GHs and their potential to affect climate.


Proposed gas hydrate impacts on past climate and their future climate change potential
Gas Hydrates (GHs) are proposed to have exerted a major control on past climates, both for the glacial-interglacial periods (Nisbet, 2002;Kennett et al., 2003), and also for other important geological events extending back to the start of the Phanerozoic or earlier (e.g., Benton and Twitchett, 2003;Kirschvink and Raub, 2004).Some infer that episodic atmospheric methane flux, attributed to massive GH destabilization, called commonly the "clathrate gun" hypothesis (Kennett et al., 2003), is the trigger for, or a significant cause of, important global climate changes including glacialinterglacial transitions.
The "clathrate gun" hypothesis is counterintuitive.It requires that methane from GHs enters the atmosphere easily once GHs disassociate (i.e., GHs are dynamically trapped primarily by clathration, rather than being formed in situ from conventional natural gas trapped previously).Methane released from GHs entrapped previously as conventional natural gas accumulations can neither enter the atmosphere nor affect climate.Ground surface or seafloor temperature changes can be inferred to precede and be greater than temperature changes of GHs at depth, since ground surface temperature (GST) variations propagate downward like a damped wave.Therefore, it seems intuitively more likely that GH destabilization follows rather than triggers climate change.In addition, ice melting and GH dissociation are endothermic and the latent heat effects should delay, additionally, any phase changes at depth, such that even where the GH destabilization mechanism is pressure reduction (Kennett et al., 2003) GH destabilization would be an unlikely major contributor to climate change.
As a contribution for evaluating potential GH climate effects generally and for understanding better the specific formation history, processes, characteristics and fate of Beaufort-Mackenzie Basin (BMB) IBP and GHs, we simulated linked profile histories (time dependant onedimensional models) of terrestrial climate (i.e., GST), subsurface temperatures, IBP and GHs in the BMB.The BMB is an insignificantly glaciated Canadian Arctic region lacking isostatic uplift.We also employed models to characterize the more recent glacial-interglacial history and the future of GHs that formed terrestrially during the Pleistocene sea level low stand, some of which were inundated and are now offshore due to Holocene marine transgression.
Our previous analysis, using a 1-D model (Majorowicz et al., 2008), of Holocene-Pleistocene temperature history implications for IBP and GH layers in onshore BMB during the last 600 ka showed that GHs can be stable through interglacial intervals despite large GST variations due to the buffering effect of overlying IBP and latent heat effects.This work, using similar models, provides an analysis of the origin, growth and persistence of BMB GH and overlying IBP layers.Models for the formation and persistence of sub-IBP GHs provide us insights into the geological history of GHs, their probable formation age, and their potential to be a source of atmospheric methane flux that might affect climates past and future.Our analysis will use model IBP and GH responses to GST histories to constrain simultaneously the models of past environments that led to IBP and GH formation and persistence.The work will also evaluate the climate risks posed by future global and regional temperature change to IBP and GH stability and their potential to serve as a future methane source to the atmosphere.

Geological setting of BMB ice bearing permafrost and gas hydrates
The geological setting is a clue to terrestrial GH formation, with important implications for past climates.Osadetz and Chen illustrated BMB GH accumulation characteristics, noting particularly that most GH resources occur in a small number of large, rich accumulations, generally overlying conventional petroleum pools (Osadetz and Chen, 2010).Although this could be, in part, due to the preferential subsurface sampling at locations favourable for conventional petroleum accumulations, it remains consistent with the apparent thermogenic petroleum system origin of BMB GH methane (Lorenson et al., 2005).In northern Alaska, where rich GHs occur partly within the IBP, there is strong evidence for stratigraphic natural gas entrapment, suggesting that a conventional petroleum accumulation was transformed in situ in response to past climate change (Collett, 1997).In the Canadian Arctic, the BMB GHs occur typically below thick permafrost, but are rare within the IBP interval.Typically the GH top occurs well below the permafrost base, commonly below a lithological barrier.The distinct and separate occurrence of BMB IBP and GHs separated by a "barren" interval is industrially very important, as the "barren" interval provides a reliably occurring and geotechnically suitable substrate for cementing petroleum well surface casing.Together, both the Alaskan and BMB GH geological settings suggest that the large, rich GH accumulations may be conventional gas pools that will remain inaccessible to the atmosphere, regardless of their physical phase.
The Mallik well sites in the Mackenzie Delta provide an excellent example where the position of the IBP and GHs contacts are well described (Dallimore and Collett, 2005) and can be used to evaluate model results.At the same location the geothermal environment and GH stability conditions are well described (Dallimore and Collett, 2005;Benton and Twitchett, 2003;Henninges et al., 2005a), both onshore and offshore (Fig. 1).Offshore there are relict terrestrially formed IBP and GH that formed during the Pleistocene sea level low stand.

Thermal environment and history reconstruction
We required a reliable GST history model to simulate IBP and GH formation and dissipation histories that include latent heat effects.Aspects of the current and more recent BMB GST forcing environment and history can be obtained from shallow temperature profiles (Henninges et al., 2005a;Rath and Henninges, 2008).However, a complete and continuous glacial-interglacial GST forcing history must be derived from sources distant from the BMB study area.Our numerical GST forcing models consider both the general cooling trend that started in the Late Miocene as well as a more detailed and variable glacial-interglacial history, as discussed below.In addition, we projected our models into the future,  considering possible changes in GST expected to accompany a doubling of atmospheric carbon dioxide.Using these GST models we simulated IBP and GH formation and their timedependent thickness variations forced by BMB climate variations (Fig. 1).

Ground surface temperature history
We employ a GST history based on several sources (Peltier and Andrews, 1976;Osadetz and Chen, 2010;Lorenson et al., 2005;Collett, 1997;Majorowicz et al., 2008;Dallimore and Collett, 2005;Judge and Majorowicz, 1992;Henninges et al., 2005a;Rath and Henninges, 2008) adapted to our model study location.Our evidence for past temperatures comes mainly from global temperature reconstructions based on isotopic data, especially δ 18 O, which provide a consistent and continuous temperature record of the last few tens of millions years (Myr).These sources show that during the interval 100 to 65 Myr, average global temperatures reached the highest values attained during the last ∼200 Myr (Lisiecki and Raymo, 2005).During that interval large-scale ice sheets were absent and forests extended to polar latitudes.Temperature also peaked 52-50 Myr ago, when atmospheric CO 2 concentrations were 1125-3000 parts per million, but by 20 Myr ago, atmospheric CO 2 concentration dropped to about 400 ppm (Polsson, 2007;Henninges et al., 2005a;Rath and Henninges, 2008).Following the Paleocene-early Eocene the climate cooled variably towards the Pleistocene glacial environment (Judge and Majorowicz, 1992;Henninges et al., 2005;Rath and Henninges, 2008).The Paleocene-Eocene Thermal Maximum (PETM) has itself been attributed to abrupt methane releases from oceanic GHs (Dickens et al., 1997).
Between 6 and 3 Myr ago, both globally averaged temperatures (Lisiecki and Raymo, 2005) (Fig. 2) and polar surface temperatures were higher than current temperatures and the Northern Hemisphere likely lacked continental glaciers (Muller and MacDonald, 2000;Collett, 1997;Henninges et al., 2005a).Climate changed dramatically, subsequently from ∼3 Myr to the present in concert with astronomical effects, specifically Milankovitch cycles, causing major GST variations that gave rise to the glacial and interglacial cyclicity within a progressively and gradually deepening Ice Age.The growth and retreat of Northern Hemisphere glacial epochs and continental ice sheets occurred at 2 frequencies.Forty-one kyr cycles were characteristic between ∼2.58 Myr ago and ∼900 kyr ago, and 100 kyr cycles were characteristic subsequently to the present (Fig. 2).The gradual intensification, or deepening, of the Ice Age over approximately the last 3 Myr is attributed to both declining concentrations of atmospheric CO 2 , which partly influences changes in temperatures, and to changes in atmospheric circulation, as El Ni (n)o is inferred to have been persistent rather than intermittent prior to ∼3 Myr ago.Significant amplification of the climate response to orbital forcing also began 3 Myr ago, resulting in the oscillations between glacials and interglacials (Fig. 2).
Our model study focuses on the BMB terrestrial onshore and shallow offshore environments at the Canadian margin   of the Amerasian Basin of the Arctic Ocean.The contiguous study area, a region inferred neither persistently nor thickly glaciated, lacks any significant post-glacial isostatic uplift (Peltier and Andrews, 1976).The study region is unlike the Queen Elizabeth Island Group in the Arctic Archipelago to the north and east, which was persistently and thickly glaciated, and where there is a strong post-glacial isostatic uplift record (Peltier and Andrews, 1976;Osadetz and Chen, 2010).
In order to simulate IBP and GH formation, we constructed a general surface temperature history for the last 14 Myr (Figs. 4-5) based on available temperature histories (Muller and MacDonald, 2000;Collett, 1997;Rath and Henninges, 2008).The most recent surface climate history for the end of the Wisconsinan and Holocene is from Taylor et al. (2005).If the current Vostok T temperature of 0 • C (Fig. 2) is adjusted to the present ground surface temperature in the study area, −6 • C (Fig. 1), then GST oscillations at the Mallik site are inferred to be between −6 • C and −4 • C during the period 6-3.5 Myr.A preliminary steady-state and the dynamic subsurface models of temperature, IBP and GHs, discussed below, use this GST model.Considering the latest Oligocene to Middle Miocene warm period (25 to 15 Myr ago), the GST model suggests that GST at the Mallik site dropped below 0 • C 14 Myr ago, initiating IBP formation.Constrained by the GST model, our simulations used a stepwise approximation of the surface temperature with a step length of 1 Myr between 14 and 6 Myr ago and 0.5 Myr between 6 and 2.5 Myr.Beginning ∼2.5 Myr ago we considered 41 ka glacial-interglacial cycles and from ∼0.9 Myr we use a 100 ka glacial-interglacial cyclity.In comparison to our previous work (Majorowicz et al., 2008), the length of the most recent glacial-interglacial cycles was reduced from ∼115 ka to ∼100 ka at the expense of the glacial part of the cycle, which is now 75 ka instead of 90 ka.
3 Preliminary steady-state pressure and temperature models of permafrost and gas hydrate formation Preliminary steady-state models indicate that GSTs boundary conditions between −6 • C and −4 • C are critical for subsurface GH formation (Fig. 3).The steady-state geotherm, corresponding to a GST of −5 • C (Fig. 3), approaches the GH phase curve from above at about 300 m, such that GH would not have formed if the GST were higher than −5 • C. Before 6 Myr ago the average climate was even warmer, which means that the first GH deposits could have formed no sooner than ∼6 Myr ago, compared to the first IBP formation which began 14 Myr ago (Fig. 4).This means that the IBP was a potentially impermeable barrier to any potential upward gas migration, long before pressure and temperature conditions permitted GH formation.The steady-state preliminary model IBP thickness is ∼250 m when GST is −5 • C.This means where gas migration to shallow depths is possible, the GHs could have formed in contact with the IBP base ∼6 Myr ago. Figure 3 shows that the geotherm and the GH stability phase curves are nearly parallel at these pressures and temperatures, such that even a small downward shift in the subsurface temperatures implies a large broadening interval, over which GH could form.Further surface cooling below −5 • C causes the IBP base to grow downward into and through the depths already occupied by GH, such that residual pore water saturation would be incorporated into the IBP after GHs formed resulting in intercalated GHs and IBP > 250 m.Thus, the preliminary steady-state model suggests that GHs would be expected to occur commonly intercalated within the IBP interval if there were no geological impediment or seal that prevented upward natural gas migration.
For 30 % porosity and 60 % GH saturation, the remaining water available occupies only 12 % of the total volume and the latent heat released by its freezing is 2.5 times smaller than that in the shallower, <250 m, intervals above the top of GHs.As BMB GHs are extremely rare within the IBP interval, even the preliminary steady-state model suggests that the observed separation and segregation of IBP above and GHs below, with an intervening barren zone, results from a lithological impediment to upward natural gas migration, nor did this gas have any potential to enter the atmosphere.
4 Dynamic pressure and temperature models of permafrost and gas hydrate formation We predicted model intervals of IBP and GH stability as functions of depth and time using an adapted computer model that simulates temporal subsurface temperature changes in response to surface forcing following the model used by Safanda et al. (2004), which itself employs an implicit finite-difference method similar to that described by Galushkin (1997).
The solution of the transient heat conduction Eq. (1) gives the temporally dependent subsurface temperature change in response to surface forcing: where T is the temperature, λ is the thermal conductivity, C v is the volumetric heat capacity, A is the rate of heat generation per unit volume, z is the depth, and t is the time in  a one-dimensional, layered geothermal model.Within the model, Eq. ( 1) is solved numerically using an implicit finitedifference method like that described by Galushkin (1997).
The upper boundary condition is the temporally varying GST and the lower boundary condition is a constant heat flow density at 15 km depth.The depth grid steps are 2, 5, 10, 50, 100, 250 and 500 m in the depth layers at 0-100, 100-1500, 1500-2000, 2000-2500, 2500-5000, 5000-10 000, 10 000-15 000 m, respectively.Time steps vary between 0.5 to 50 yr, depending on the amplitude of GST changes.The finite-difference scheme of Eq. (1) on the depth and time grids ... z k−1 , z k , z k+1 , ... and ... t n , t n+1 , ..., respectively, has a form: where the subscript k and the superscript n denote values at the k-th depth step and the n-th time step, respectively; This difference scheme, together with the upper and lower boundary conditions, leads to a system of difference equations for unknown values k+1 within a tridiagonal matrix, which was solved by the Peaceman and Rachford's forward method (Peaceman and Rachford, 1955).
To estimate effective thermal conductivities and volumetric heat capacities, it was necessary to consider the respective geometric and arithmetic averages of the constituent values for the rock matrix, water, ice and GH in proportion to their volumetric fractions (Galushkin, 1997;Taylor et al., 2005).A consumption or release of the latent heat, L, in water/ice (334 kJ kg −1 ) and GH (430 kJ kg −1 ) accompanying either thawing or freezing was included.The interstitial ice and GH effects were accounted for using apparent heat capacity, C L , according to Carslaw and Jaeger (1959), when the volumetric heat capacity increases in the model depth sections, where the thawing and freezing occurs (i.e., where the temperature is within the thawing range between the temperature of solidus T S , and liquidus, T L , at the actual simulation time step) by the term C L = ρ L/(T L − T S ).Here, ρ is the density of either ice or GH and is a fraction of the total volume occupied by these phases.The model assumes a simple hydrostatic pressure gradient in the water-wet sediment pore space.The liquidus and solidus temperatures of water/ice and GH are depth and pressure dependent (Galushkin, 1997) and solidus temperatures were 0.2 • C lower than liquidus temperatures.
In the IBP zone we infer 30 % of rock matrix porosity to be fully occupied by water at temperatures above T L , and by ice at temperatures below T S .Within the GH stability zone, the GH saturation in matrix porosity was inferred to be 60 %.A salt concentration of 9 g l −1 was considered constant with depth and the p-T phase curves were so adjusted.The salinity of 9 g l −1 (Galushkin, 1997) was used so that the liquidus temperature at the permafrost base in Mallik (−1 • C at 600 m) corresponds with the value given by formula: If there is a fresh water within Mallik sediments, the liquidus temperature would be 0.58 • C higher.If there is sea water with 40 g l −1 , the liquidus would be 1.98 • C lower.For temperature gradient 20 K km −1 , it means a shift of the permafrost base 30 m downward and 100 m upward, respectively.
Model calculations were tested and verified by comparing results against the analytical solidification problem solution (Carslaw and Jaeger, 1959).The error estimate for the IBP and GH numerical simulations is of the order of tenths of a degree C. A similar error range was estimated by halving the model time and depth steps.
Our model uses deep heat flow, thermal conductivity, present IBP and Type I GH depth and thicknesses, a surface melting temperature (−0.576 • C) and observed groundwater salinities (9 g l −1 ).It employs latent heat effects throughout the IBP and GH layers, which is an improvement upon previous models (Taylor et al., 2005).Deep heat flow is constrained by bottom hole temperatures from deep wells (Majorowicz et al., 1990;Taylor et al., 1982), and wellconstrained estimates of thermal conductivity, latent heat, present IBP thickness and observed GH thicknesses (Henninges et al., 2005b;Wright et al., 2005).The models consider the hydrostatic pressure-depth dependence of ice and GH thawing points over the entire expected extent of the model IBP and GH layers (Sloan, 1998)(Fig.3), whereas previously published models (Taylor et al., 2005) considered only a thin layer using a constant dissociation temperature.
In simulating the models temporal subsurface changes with vertically separated occurrences of IBP and GH, a consumption or release of the latent heat during decay or formation of IBP or GH was considered separately.This means that in the IBP zone, only the latent heat and p-T phase curve of water ice were taken into account, and in the GH zone only latent heat and p-T phase curve of the GH is considered.The division into IBP and GH zones was prescribed explicitly in the model (see also Taylor et al., 2005).Such a separate treatment of ice and GH latent heats may not be appropriate in models of intercalated IBP and GH occurrence.Where ice and GH are intercalated, our models used the relative fractions of ice and GH in the pore space, which are prescribed specifically at each depth.The code then checks the P -T conditions for IBP and GH at each depth independently, during each time step, and it adjusts the volumetric heat capacity to the actual temperature at that depth (Carslaw and Jaeger, 1959).At Mallik, where 30 % porosity was assumed, the prescribed pore GH fraction in the uppermost 250 m was zero and that of IBP 100 %.Where the GH pore space saturation is 60 %, it means that the pore fraction ratios occupied by GH and IBP are 60:40.

Model 14 Myr surface and subsurface temperature histories and implications for IBP and GHs
Subsurface temperature history depends on the GST forcing model, the downward propagation of GST variations and the latent heat effects of phase changes in both the IBP and GH layers.In our GST forcing model, the surface temperature at Mallik drops below 0 • C no earlier than 14 Myr ago, initiating IBP formation.Subsequently, we use a stepwise surface temperature approximation with a step length of 1 Myr between 14 and 6 Myr ago and 0.5 Myr between 6 and 2.5 Myr; from 2.5 Myr ago we consider the 41 ka cycles and from the 0.9 Myr, the 100 ka cycles.
Using this GST model we simulated the downward propagation of GST variations on eastern Richards Island, the Mallik gas hydrate well site, and its implications for IBP and GH formation and dissipation using a GST model (Fig. 1).Individual computational models use the characteristics of IBP and GH formation and dissipation as functions of temperature history, constrained by present temperature observations and current IBP and GH layer thicknesses (Dallimore and Collett, 2005).Thermal conductivity dependence on rock section specific heat as a function of both porosity and interstitial water and ice proportions is important.It is also crucial to account for latent heat effects attending the freezing and thawing of interstitial water or ice in the IBP layer to match observations successfully.All models account for latent heat using the apparent specific heat, which is a standard treatment.The models also consider diffusive heat flow  related to surface-subsurface coupling.All models permit IBP and GH to occur simultaneously, either separately or intercalated.Subsurface temperature changes would proceed much faster, generally without the heat sink provided by IBP latent heat.
We simulated two model "cases" of IBP and GH accumulations using the Mallik GST model.In Case 1 (Figs. 4 and 5) we allowed GHs to form due to an in situ transformation of a pre-existing conventional natural gas accumulation trapped below a geologically impermeable seal.This case can be viewed as describing a GH occurrence end-member for the in situ transformation of a pre-existing conventional gas pool.In Case 2 (Figs.6-9) we simulated a situation lacking a lithological seal that inhibits upward natural gas migration.This case can be viewed as a GH occurrence end-member for natural gas that is dynamically trapped by GH formation itself.Both these end-member models have a geological basis, depending on the local lithological composition of the sedimentary succession along petroleum migration pathways.There are reasonable grounds to infer that both cases are possible in the heterolithic and lenticularly bedded, clastic deltaic deposition setting of the BMB (Majorowicz et al., 2008) even though IBP and GH are generally not intercalated and GH accumulations are not distributed uniformly either geographically or stratigraphically.Rather, GHs are strongly associated with underlying conventional petroleum pools (Osadetz and Chen, 2010).

Case 1 -gas hydrate forms in a conventionally entrapped natural gas accumulation
Case 1 is the most probable situation for large, rich BMB GH accumulations (Dallimore and Collett, 2005).Both at Mallik  and at other BMB GH accumulations, the GHs occur distinctly and separately below the IBP with an intervening "barren zone".GHs are hosted in coarse-grained sandstone beds overlain and separated by thin non-hydrate-bearing, fine-grained siltstones and mudstones within an anticline trap (Dallimore and Collett, 2005;Osadetz and Chen, 2010).A fault bounding the structure is the likely conduit for interstratal gas migration into the reservoir, which now lies within the GH stability zone (GHSZ).GH methane isotopic compositions, ∼ −42.7 %, are indicative of a thermogenic petroleum system (Sloan, 1998).The upward migrating gas was likely trapped in the anticline long before (Kroeger et al., 2008) the geologically recent cooling, which subsequently formed the GHs.This conventional accumulation hypothesis is also consistent with the coincidence of GH occurrences above deeper conventional petroleum reserves (Osadetz and Chen, 2010;Lorenson et al., 2004).To calculate Case 1 the vertical division and separation of the IBP from GHs was prescribed explicitly in the model, as expected for gas trapped conventionally below a lithological seal.The consumption or release of latent heat during the decay or formation of IBP or GH is likewise considered separately, such that, in the IBP zone, only the latent heat and p-T phase curve of water/ice were taken into account, while in the GH zone only the latent heat and p-T phase curve of GH were used.
Case 1 results (Fig. 4) indicate that GH started to form only after the onset of the 100 ka cycles 0.9 Myr ago, when the glacial interval surface temperature dropped to −15  1060 m (Dallimore and Collett, 2004), are near but slightly above their observed positions at Mallik, which are 600 m and 1107 m, respectively.In previous simulations (Majorowicz et al., 2008), these bases occurred deeper, at 600 m and 1165 m.The difference is attributed to a change in the GST model, specifically the shorter duration of the cold glacial intervals, which are 75 ka in the current simulations, instead of the 90 ka used previously (Majorowicz et al., 2008).
Considering the quite generalized GST history (Fig. 4), the agreement between model and observed GH and IBP depths is very good, suggesting that the GST model is consistent with the only observable constraints on IBP and GHs.Details of the entire 14 Myr history of GST forcing effects for IBP and GH (Fig. 5) show that the IBP has persisted since its initial formation shortly after 3 Myr and that the GH layer has persisted since its first formation below the lithological seal at 900 m shortly after the 100 ka cycles began.

Case 2 -gas hydrate forms where upward gas migration is unimpeded by lithological barriers
Case 2 simulates a dynamic GH trap, where natural gas can migrate freely toward the surface without a lithological seal.
The model permits the simultaneous occurrence of intercalated IBP and GH within the pore space.This case resembles the fate of the majority of petroleum generated from thermogenic petroleum systems, as typically only a small fraction of the petroleum generated is actually entrapped (e.g.36).The Case 2 simulation considers the last 14 Myr (Figs. 6-9) using the same GST history used in Case 1 above.
Prior to 14 Myr ago, the Mallik site GST was +1 • C, and subsurface temperatures remained above p-T phase curves of both IBP and GH at all depths, where the basal heat flow is 60 mW m −2 and conductivity is 2.1 W/(m K −1 ).With the onset of a gradual cooling 14 Myr ago, the GST dropped below  the solidus IBP temperature and permafrost started to propagate downward from the surface.However, for our model the frozen thermal conductivity is 3.4 W/(m K −1 ) when all the pore water is frozen (30 % porosity).Using a latent heat of water freezing = 0.334 MJ/kg, the subsurface temperatures stay above the GH p-T phase curve until the GST drops to −5.5 • C, as a GST of −5 • C was not low enough to permit GH formation.It is uncertain what might have happened to methane within the freezing rock, which should have been there.If all the pore water had been frozen, methane would have been prevented from migrating to the surface and the natural gas might have been displaced downward, below the downward migrating IBP base or either escaped at the margins of the IBP to resume its migration upward, as this might occur currently where IBP "weak" thermokarst occurs around shallow Mackenzie River Delta lakes (Bowen et al., 2008).
Figure 6 shows transient temperature-depth profiles corresponding to a model time 42 500 years after 6 Mya ago, when the GST cooled from −4.5 to −5.5 • C, and at a model time 20 years after the subsequent warming back to −4.5 • C 5.5 Myr ago.After the cooling 6 Myr ago, the T -z profile started to move downward and touched the p-T phase curve for GHs at 290-300 m depth 42.5 ka after the cooling (i.e., the light blue line in Fig. 6 at model time 5.9575 Myr ago).When GHs first began to form, the IBP base was at 250 m.Thus, at that time, there was an unfrozen intervale between 250 m and 290 m thick that separated IBP from GHs.With further cooling, the GH layer propagated not only downward, but also upward, such that after an additional model time interval of 100 ka, the model upper GH boundary touched the IBP base at 255 m.Upward growth of the GH zone then stoped, because it contacted the ice in the already frozen IBP pore space.Nevertheless, the permafrost continued migrating downward through the permafrost controlled GH top to freeze the residual pore space water (i.e. the model considers 60 % GH saturation in the pore space and the remaining 40 % of the pore space, or 12 % of the entire volume, remains water-filled).It was that interstitial water that froze with continued cooling and downward propagation of the IBP layer.
This scenario and its consequences for IBP and GH occurrence were accounted for by considering separately water/ice-filled and gas/GH-filled porosity.For ice, this was 30 % by volume above the IBP-GH contact at ∼250 m and 12 % below it due to the prior consumption of pore space water as GHs formed.Similarly, the porosity available for the GH above the IBP-GH contact at 250 m was 0 % and 18 % (60 % GH saturation in 30 % porosity rock) below it.This presented the possibility to treat latent heat changes in the solidus-liquidus zones of IBP (0.334 MJ kg −1 ) and GH (0.430 MJ kg −1 ) both independently and appropriately.The thermal conductivity in the zone of intercalated IBP and GH was estimated at 2.7 W/(m K −1 ), which is lower than the value in pure IBP (3.4 W/(m K −1 ) because of the presence of the intercalated GH fraction that has a low conductivity = 0.45 W/(m K −1 ) compared to waterice = 2.2 W/(m K −1 ).
The Case 2 simulation results indicate that the first GH, which began forming 5.9575 Myr ago, reached its thickness of 125 m (between 255-380 m) shortly after the GST warmed from −5.5 • C to −4.5 • C, 5.5 Myr ago (Fig. 7).Note specifically that for Case 2, the GH layer continued to grow downward even after the GST warming began because of the time lag between GST changes and temperature changes at depth.The same is true for the IBP, which penetrated the GH layer only 4 m (255-259 m). Figure 7 shows that, following GST warming, the last GH, in the depth range 310-340 m, disappeared at 5.472 Myr ago when it was separated from the IBP base at 240 m by a 70 m thick melted zone.GH formed again, shortly after the next GST decrease from −5.0 • C to −5.5 • C at 4.5 Myr ago.Since that time, some GH has persisted through all subsequent times and GST variations.The top of the persistent GH layer at 250 m, is always within the GH stability zone, while it extended downward to a maximum depth of 990 m during the most recent glacial interval.Depth variations of the simulated IBP and GH upper and lower boundaries during the entire modelled interval of 14 Myr are depicted in Fig. 8. Details of Fig. 8 for the last 3 Myr are depicted in Fig. 9. Case 1 results are shown for comparison (Figs. 8 and 9).
The shallower than currently observed position of both the model IBP and GH bases in the Case 2 model results from the lower total conductivity and consequently, higher temperature gradient within the IBP, below 250 m, where IBP is intercalated with GH. Figure 9 also shows that the depth variations of both the IBP and GH bases are slightly larger for Case 2 than for Case 1 during both the 41 ka and 100 ka glacial-interglacial cycles.This is attributed to the larger variations of the IBP base exerting a relatively smaller latent heat dampening effect due to the GHs intercalated within the IBP below 250 m.The smaller damping effect at the IBP base also results in the faster propagation of subsurface temperature to the GH base that results in the larger Case 2 depth variations of GHs compared to the Case 1 results.
In Case 2, methane is present in the pore space from at least below 250 m, where it transforms to GH prior to the arrival and downward propagation of the IBP base such that: 1. IBP thickness did not exceed 200 m before 6 Myr ago.
2. The subsurface temperatures accompanying GST temperatures was > −5.5 • C prior to 6 Myr ago precluding GH formation at all depths.
3. The first GH formed at 290-300 m shortly after 6 Myr ago when the GST dropped from −4.5 • C to −5.5 • C. The GH layer subsequently expanded both downward and upward until 100 ka later when it contacted the downward propagating IBP base at 255 m.
4. The GH accumulation formed initally disappeared completely during the subsequent warm interval, but GHs reformed shortly after GST again dropped from −5 • C to −5.5 • C at ∼4.5 Myr ago to persist to the present.The upper contact of the GH stability zone at 250 m, which was intitally constrained by contact with the IBP has remained within the GH stability zone to the present and the maximum downward GH layer extent, 990 m, was reached during the most recent glacial invterval.
5. Shortly after the formation of a GH layer, the downward migrating IBP base and upward migrating upper GH boundary met at 250 m.Since that time, IBP forming subsequently below 250 m has been intercalated with GHs and pore space ice could form from only the residual pore water that was not already bound in GHs.
Although based on the assumption that there is no lithological barrier to upward methane migration, the Case 2 model also provides no opportunity for methane to escape to the surface after the IBP "lid" forms a dynamic and persistent impediment beginning at 14 Myr ago.-9, which indicates the upper potential limit of GH stability.However, it does not imply that these previously ice-filled sediments contain model GHs, since the model precludes pore space GH formation, where complete freezing of pore space water preceeded GH stability at a specific depth (i.e., model GH does not form where the pore space is already occupied by ice).
6 Sensitivity of the model to uncertainty in the history of surface temperature change Our calculations require a continuous numerical GST model.
In the absence of other comprehensive and continuous local GST models, we used a global air surface temperature reconstruction constrained by the recent surface temperature history of the study area as a "plausible" surface forcing model.This model allowed us to simulate how glacial-interglacial temperature variations affect permafrost stability and how the latent heat effects preserve underlying GHs, the phase changes within which have their own intrinsic latent heats.
Our GST forcing model is consistent with global GST variations; it is forced to match current study area GST; it is consistent with recent thermal history derived from shallow borehole temperature anomalies (Henninges et al., 2005a;Rath and Henninges, 2008), and it predicts both the currently observed IBP and GH bases (Dallimore and Collett, 2005).Still, we acknowledge that our GST forcing model is also simple and that it is inadequately calibrated for the studied area, generally qualitative indications of previous temperature histories, especially for the older parts of the model, prior to the 100 ka cycles.Unfortunately, there is no wellestablished, or consensual local GST model that covers sufficiently our simulated interval.Therefore, we conclude that it is both necessary and entirely appropriate to adapt the well studied and more globally relevant temperature histories and to employ them, suitably adjusted for location and current GST state.The use of regionally and globally consistent temperature histories provided testable predictions consistent with currently observable conditions (Kennett et al., 2003;Henninges et al., 2005a;Rath and Henninges, 2008).Still, the extent to which local GST history followed regional and global GST variations, especially at earlier times (14.0 to 1.0 Myr), remains uncertain, especially for the more distant past.
To explore GH formation and stability sensitivities to uncertainties in the early GST forcing model, we made additional simulations for Case 1 (Fig. 5) prior to the 100 ka glacial cycles (Fig. 10).The sensitivity analysis simulation begins 3 Myr ago with a steady-state T -z profile corresponding to GST of 0 • C. In the sensitivity analysis test, the GST decreases linearly from 0 • C 3 Myr ago to −10 • C at the onset of the 100 kyr glacial cycles 0.9 Myr ago.The resulting model is up to 6 K higher for most of the period between 3 Myr and 0.9 M difference between these two models is tens of meters, at most.It indicat youngest part of the model is not very sensitive to the early parts of the G The "warmer" model is up to 6 K higher for most of the period between 3 Myr and 0.9 Myr ago.The current difference between these two models is tens of meters, at most.It indicates that the youngest part of the model is not very sensitive to the early parts of the GST history.
differences between the linearly declining or "warm" GST model during the last 3 Myr and "cold" GST models that include a 41 ka cyclicity is observed primarily to affect the IBP base prior to ∼1 Myr ago.The effect for Case 1 GH differs by only tens of metres during the past 0.9 Myr.Therefore, we conclude that uncertainties in earlier parts of the GST forcing model have an insignificant impact on the Case 1 model results during the more recent past and present.Despite our use of a general and regional GST history, our simulation results agree well with the currently observed IBP and GH bases.Neither do the variations in the earlier parts of the GST model produce significant differences between model results subsequent to the beginning of the 100 ka cycles, nor do they require that we modify significantly any of the conclusions drawn from simulations that employ the regional GST history.
The fact that the Ice Age has deepened progressively means that the current IBP base is near its deepest extent today, effectively obliterating any potential to find and possibly date previous IBP bases that would distinguish the "warm" or the "cold" GST histories of the sensitivity analysis.However, the sensitivity analysis does suggest that earlier temperature histories result in different IBP bases.Since the IBP base can provide a "lid" Case 2 GH occurence, it might be possible to use the top of Case 2 GH occurrences to additionally calibrate the GST history in specific GST scenarios.However, as mentioned above, the Case 2 GH occurrences are rare and poorly substantiated (i.e., GH as opposed to free gas in the IBP) in the BMB.

The impacts of Holocene marine transgression and anticipated future warming accompanying increased atmospheric CO 2 content
The two cases simulated above illustrate that IBP and GH models, predicted by a global GST history, model successfully current observations.Still, the models explain the main features of terrestrial BMB GH occurrences including the general lack of intercalation of GH within the IBP and generally deeper GH top contact that would be expected from the convolution of equilibrium phase thermodynamics without the constraints provided by either lithological or permafrost seals.
Since the end of the last glacial maximum, two processes have impacted or could impact, past and future forcing temperature boundary conditions with implications for IBP and GHs.The first effect is Holocene marine transgression over terrestrial GHs formed during the Pleistocene low stand.The second effect is the anticipated warming accompanying increased atmospheric CO 2 content.We illustrated both these effects sequentially using a single simulation that examines, sequentially, the results of changing the temperature 34 Fig. 12. Changes in the onshore and offshore (Holocene marine transgression) GH and IBP during the present glacial cycle including the possible future warming (see Fig. 11) in its end induced by a doubling of atmospheric CO2 content.The Holocene starts at time 550 ka of the simulations and the present time is 563.5 ka.

Fig. 12.
Changes in the onshore and offshore (Holocene marine transgression) GH and IBP during the present glacial cycle including the possible future warming (see Fig. 11) in its end induced by a doubling of atmospheric CO 2 content.The Holocene starts at time 550 ka of the simulations and the present time is 563.5 ka.
boundary conditions accompanying Holocene marine transgression beginning 13.5 ka and followed by future warming anticipated due to changes in atmospheric CO 2 content (Intergovernmental Panel on Climate Change, 2007).The effects on IBP and GH model histories for these two surface forcing effects are illustrated by a detailed model of the last ∼0.5 Myr (Figs. 11 and 12).
This detailed simulation (Figs.11 and 12) considers a proposed temperature history that includes both marine transgression and future warming following Taylor et al.'s previous discussion (Taylor et al., 2005, their Fig. 3b).The GST model is similar to our previous model 3 (Majorowicz et al., 2008) and it uses a glacial-interglacial periodicity of 115 ka and a frozen versus thawed conductivity varying from 3.4 (frozen) to 2.1 (unfrozen) W/(m K −1 ).
The simulation considers Holocene marine transgression (Fig. 11) as it affected the previously exposed, Pleistocene continental shelf that is now inundated by the shallow Beaufort Sea.To consider both these effects, we have used the surface temperature model for Mallik (Taylor et al., 2005;Nixon, 1986) and our previous model (Majorowicz et al., 2008), at the end of the last glacial.The model warms from −2 • C to +1 • C (Fig. 11) due to marine inundation instead of the terrestrial cooling to −5 • C that occurred 8 ka ago (Taylor et al., 2005;Nixon, 1986).In addition, we also consider a future warming, starting at present, from +1 • C to +7 • C during the next 300 years consistent with a suggested warming at a rate of 2 • C per century accompanying the projected doubling of atmospheric CO 2 during the same interval (Intergovernmental Panel on Climate Change, 2007).We are aware that marine transgression also occurred during previous interglacials, however, the record of their extent is neither well preserved, nor well documented and the results of www.clim-past.net/8/667/2012/Clim.Past, 8, 667-682, 2012 the sensitivity analysis suggest that previous departures from the GST model affect the most recent history insignificantly.Therefore, we consider only the Holocene marine transgression, although with suitable GST modifications the simulations could include all such effects.
Figure 13 shows the temperature-depth profiles resulting from a model based on Taylor's surface temperature history until marine transgression was followed by anticipated future warming effects (Taylor, 1999).Because the transgression results in a shallow water column, the change of the hydrate pressure and temperature phase curve, due to hydrostatic pressure increase attending the transgression being neglected.The phase curve corresponding to 50 m water depth is shown by a dashed line (Fig. 13).After the ∼50 m sea level increase, the hydrate phase boundary shifts downward ∼25 m and the heat released by GH formation increases the temperature in that zone, although not significantly.
When this simulation is projected into the future (Fig. 12) to the hypothetical "natural" end of the current interglacial, about 11.5 ka in the future, the simulation predicts that only a thin vestige of the IBP will remain due to thawing by ∼150 m from below and 70-80 m from above.However, most of the thinning occurs in the IBP, which buffers the effect on the GH layer, such that the accompanying GH layer thinning is very small and within the range attributed to previous interglacial variations, in spite of the accelerated surface warming that is anticipated to accompany a doubling of atmospheric CO 2 .

Discussion
We used four models that illustrate BMB, terrestrial IBP and GH formation and dissipation as a function of GST history, which includes a realistic treatment of phase change latent heat.Initial models examined the IBP and GH formation and dissipation employing a locally adjusted, globally consistent, GST applied to two geological end-member GH accumulations model cases: Case 1, in situ GH formation from natural gas previously entrapped under a geological seal; and Case 2, GH formation from a freely upwardly migrating natural gas entrapped dynamically by clathration.These two models have distinctively different GH and IBP accumulation characteristics despite their common GST model.In both cases GH forms only after IBP.The most noticeable differences are the depth and age at which GHs form.Case 1 GHs formed at 0.9 km depth no earlier than about 0.9 Myr ago, roughly coincident with the onset of the 100 ka glacial-interglacial cycles.Case 2 GHs first formed in the depth range 290-300 m shortly after 6 Myr ago when the GST dropped to −5.5 • C.
Case 1 GHs grew downwards subsequently, limited by the lithological top seal of the conventional petroleum accumulation.Case 1 GHs persist to the present, although the GH layer base changes depth following, but out of phase with GST variations, such that the GH layers are still expanding downward even after the GST warming begins.This is due   to the delay in GST changes affecting the deeper subsurface temperature field.Case 1 suggests that IBP and GH will occur separately and distinctly, possibly separated by a barren zone containing neither IBP nor GH, if the depth of the lithological seal occurs below the IBP base.
Case GHs initially grow both downward and upward, as they are unconstrained by lithological seals, however, their upward growth continues only until they contact the downward thickening IBP layer, which serves as an impermeable barrier that mimics the effect of the Case 1 lithological seal at a much shallower depth.Case 2 GH that formed about 6 Myr ago did not persist past 5.5 Myr ago, due to a warmer interval, but Case 2 GHs were re-established and have persisted since the progressive cooling trend that began about 4.5 Myr ago.Case 2 GHs also respond out of phase to the changing GST, especially after the onset of the 41 ka glacial interglacial cycles to the present.Likewise, Case 2 IBP grows downward and it similarly responds by thickening and thinning as a delayed response to GST fluctuations.Case 2 IBP grows downward to freeze the residual water in pores already partly filled by GHs.Case 2 models are characterized by an absence of GHs at depths shallower than 255 m and that intercalated GHs and IBP should be common in the interval between 255 m and below.The intercalation of GHs and IBP changes the thermal conductivity and latent heat of that depth interval compared to Case 1 models such that at a given time the Case 1 IBP and GH bases are generally deeper than for Case 2.
The differences between Case 1 and Case 2 current predictions (Fig. 9) indicate strongly that Mallik and most other BMB gas hydrate accumulations are Case 1-like GH that formed by in situ transformation of previously entrapped conventional gas accumulations below a lithological seal.If Case 2 models prevailed, it would be common both to find GH at much shallower depths as high as but no shallower than 255 m and that both IBP and GHs would be commonly intercalated through the depth interval between 255 m and the IBP base.Our models lead us to conclude that Case 1 GH occurrences predominate in the BMB since GHs generally occur separately and distinctly from the IBP.Regardless of its phase, Case 1 natural gas was never and will not be available to influence climate until the conventional petroleum trap fails.However, the surprising result of our simulations is that neither is the Case 2 methane available to migrate into the atmosphere and influence climate because of the IBP that formed first at 14 Myr ago and which completely froze pore space fluid prior to the onset of GH stability probably.This provides as effective a seal to upward gas migration as that provided by lithological barriers (Fig. 8).
During the ∼100 kyr cycles, GH layer thickness generally increases during colder intervals (i.e., glacial) and decreases during warmer intervals (i.e., interglacial) and these variations are ∼0.2 km for the IBP and about 0.1 km for the GH layers, although the delayed effect of GST at depth combined with latent heat effects means that these changes follow rather than lead GST changes.Where the IBP layer is sufficiently thick, it is unlikely that any sub-IBP GH accumulations disappeared entirely during previous interglacial intervals.
Our third model examined the effects of a warmer GST history on model results.A warm distant past has little effect on present Case 1 IBP and GH occurrences, although the effect prior to ∼1 Myr ago is significant.It is possible that the rare Case 2 type of GH accumulations could exhibit greater sensitivity to prior GST variations, although such occurrences remain unsubstantiated.
Our fourth model combines and examines both Holocene marine transgression and anticipated GST warming accompanying a doubling of atmospheric CO 2 .Both marine transgression and anticipated GST warming do affect the models, but neither effect is significant.Neither are these effects expected to result in complete IBP disappearance, where it was sufficiently (∼500-600 m) thick.The attending GH thickness variations are "protected" by the latent heat effects of IBP thawing, such that prior to the "natural" end of the current interglacial, GH thicknesses vary no more than during any of the previous 100 ka cycles.
Our models suggest that GH layers can act as a persistent sink for methane from deep thermogenic petroleum systems, where terrestrial IBP is thick and that stratal barriers to upward gas migration can be augmented by persistent thick IBP, which is also an effective barrier to the migration of methane into the atmosphere.Our models also illustrate how IBP and GH layers have a delayed phase transformation response to GST variations, both because of how GST changes are filtered by the Earth at depth and due to the latent heat requirements of thawing and GH dissociation.
The hypothesis that GHs destabilize rapidly in response to environmental change, late in glacial intervals, and that they serve, at other times, as a sink for and barrier to the migration of methane into the atmosphere has been attributed primarily to marine non-IBP GHs (Nisbet, 1990(Nisbet, , 2002;;Kennett et al., 2003;Benton and Twitchett, 2003).We do not now address that setting with these models and we acknowledge that marine GHs may be more easily destabilized than the terrestrial sub-IBP GHs we have analyzed here.
Our study shows that sub-IBP GHs below vary in thickness in response to surface temperature history changes, but that terrestrial thermal inertia conserves both IBP and GHs, delaying and reducing methane formation, which is also detained below persistent lithological and IBP seals, such that GH methane from neither Case 1 nor 2 sources can enter the atmosphere, since BMB permafrost first formed perhaps ∼14 Myr ago, although other mechanisms do result in active methane seepages (Bowen et al., 2008).Terrestrial thermal inertia also imposes a phase-delay between surface temperature warming and the subsequent onset of GH dissociation, making it unlikely that terrestrial GHs below thick IBP would either trigger, or rapidly reinforce, climate warming events.The implications of latent heat and thermal diffusivity effects for submarine GHs on the continental shelves remain to be determined, however, our model results appear consistent with observations of end-Pleistocene methane anomaly isotopic compositions from ice cores that inferred inconsistently with the "clathrate gun" hypothesis (Sowers, 2006).
We recognize the limitations of a one-dimensional modelling approach, despite its apparent success prediction of ibp and gh stabilities, using a very restricted range of used physical properties and implied simplifications in upper and lower thermal boundary conditions.It is possible that a one dimensional (1-D) model does not adequately represent situations where thin permafrost will be much more sensitive to surface warming and where the possibility of ibp buffering of surface warming to preserve underlying ghs is reduced, as is present in some parts of the bmb especially near the margins of the thick ibp area.
We also recognize opportunities to improve both the model and the GST history, as well as limitations attending model uncertainties, especially as they relate to our GST model.There is no well established or consensus model for a protracted and quantitative local temperature history.However, the current prediction obtained using an adjusted global GST history provided an encouraging model that is in general agreement with current observation.Another forcing model limitation could be due to the assumption of the ground surface condition.In some cases, surface temperature history can differ from GST history, such as the difference between the temperature at the base of an ice sheet compared to its upper surface (Marshall and Clark, 2002).Similarly, the geothermal heat flow and the rate of vertical ice advection can affect basal ice sheet temperatures (Demezhko et al., 2007).We have focused on terrestrially formed BMB GH www.clim-past.net/8/667/2012/Clim.Past, 8, 667-682, 2012 occurrences in an area that was neither persistently nor thickly glaciated as it lacks significant post-glacial isostatic uplift.This permitted us to make a simpler model that is free of the effects of both the pressure and temperature of more active ice dynamics.
As a result, we restrict our conclusions regarding the significance of sub-permafrost GH models to the terrestrial BMB and in detail to the Mallik well area of the Richards Island.Neither do we generalize our models to the entire Arctic region, where variety of paleo-environmental conditions and geological settings, including ice sheet loading pressure effects may provide different results.One such specific condition can be the rapid change of temperature from cold terrestrial to warm sea bottom conditions that have been attributed to coastal IBP and GH destabilization in the Arctic Islands (Majorowicz, 1998;Wright et al., 2005).

Conclusions
This study shows that thermal history modelling as a function of GST history contributes significantly to our understanding of the origin, growth, persistence, and characteristics, of sub-IBP GH accumulations.It helps to determine and distinguish the features of GH geological natural gas resource model, with implications for both climate and resources.The key conclusions are 1.A GST history based on marine isotopic data adjusted for study location predicts both IBP and GH current bases as observed.
2. Two end-member GH models, where gas is either trapped conventionally or where it is trapped dynamically by GH formation have distinctive characteristics that lead us to conclude that most BMB GH accumulations formed coincident at the onset of the 100 ka cycles by the in situ transformation conventional natural gas accumulations, based on the general distinct and separate occurrence of GHs and IBP layers and the IBP lacks GH intercalations.Uncertainties in the early GST history do not affect significantly either current results or their interpretation.
3. Case 2 GH models formed ∼6 Myr ago, but these occurred only >255 m, below ice-filled permafrost.The persistent IBP layer is as effective a seal to upward methane migration as are lithological seals, such that methane from sub-IBP GHs has been prevented from entering the atmosphere since ∼14 Myr ago.
4. All models respond to GST variations, but in a delayed and muted manner such that GH layers continue to grow even as the GST begins to increase.Both thick IBP and GHs have persisted since ∼1.0 Myr ago and ∼4.0 Myr ago for Cases 1 and 2, respectively and they persisted 5.Where permafrost is sufficiently thick, combinations of geological structure, thermal processes and material properties make sub-IBP GHs unlikely sources for significant atmospheric methane fluxes.
6.Although we do not address the marine GH accumulations attributed to driving the "clathrate gun", our models and the details of sub-IBP GH histories suggest that similar models could provide a robust understanding of GHs and their potential to affect climate.

Fig. 1 .
Fig.1.Depth to GH base of Type I stability in the Mackenzie Delta.Model studies were performed for northern Richards Island and its adjacent offshore (near the centre of the figure).

Fig. 1 .
Fig. 1.Depth to GH base of Type I stability in the Mackenzie Delta.Model studies were performed for northern Richards Island and its adjacent offshore (near the centre of the figure).

Fig. 3 .
Fig. 3. Steady-state profiles for Mallik wells on Richards Island area geothermal model.Conductivity of the frozen and thawed rock is assumed to be 3.4 and 2.1 W m −1 K −1 , respectively.Note: the figure shows that GH is stable only for steady state surface temperatures < −5 • C for the terrestrial heat flow and thermal properties of the Beaufort Mackenzie Basin.

Fig. 4 .Fig. 4 .
Fig.4.Case 1 model results illustrating the surface temperature and IBP and GH bases during the last 3 Myr according to a numerical simulation of subsurface temperatures for the last 14 Myr of surface temperature forcing Fig. 4. Case 1 model results illustrating the surface temperature and IBP and GH bases during the last 3 Myr according to a numerical simulation of subsurface temperatures for the last 14 Myr of surface temperature forcing.

Fig. 5 .
Fig. 5. Case 1 model results for the entire 14 Myr history of surface temperature forcing showing both the time dependence IBP and GH base depth variations.

Fig. 5 .
Fig. 5. Case 1 model results for the entire 14 Myr history of surface temperature forcing showing both the time dependence IBP and GH base depth variations.
Fig.6.Case 2 model transient temperature-depth profiles illustrating the response to cooling the GST from -4.5 °C to -5.5 °C at 6 Myr ago (and to all previous GST history since 14 Myr ago) followed by a warming back to -4.5 °C 5.5 Myr ago.

Fig. 6 .
Fig. 6.Case 2 model transient temperature-depth profiles illustrating the response to cooling the GST from −4.5 • C to −5.5 • C at 6 Myr ago (and to all previous GST history since 14 Myr ago) followed by a warming back to −4.5 • C 5.5 Myr ago.

Fig. 7 .
Fig. 7.Case 2 transient temperature-depth profiles illustrating the response to the sudden GST warming from -5.5 °C to -4.5 °C at 5.5 Myr ago and to all previous GST changes since 14 Myr ago.The profiles correspond to times 20, 10 000, 28 000 and 40 000 years after the warming begins.The curve for time 28 000 years (5.472Myr ago) shows the moment when the last GH disappears in the depth interval 310 -340 m, when it is separated from the IBP base (240 m) by a 70 m thick melted zone.

Fig. 7 .
Fig. 7. Case 2 transient temperature-depth profiles illustrating the response to the sudden GST warming from −5.5 • C to −4.5 • C at 5.5 Myr ago and to all previous GST changes since 14 Myr ago.The profiles correspond to times 20, 10 000, 28 000 and 40 000 years after the warming begins.The curve for time 28 000 years (5.472Myr ago) shows the moment when the last GH disappears in the depth interval 310-340 m, when it is separated from the IBP base (240 m) by a 70 m thick melted zone.

Fig. 8 .
Fig. 8.A comparison of Case 1 and Case 2 results illustrating the respective depth variations of IBP and GH upper and lower boundaries during the last 14 Myr.

Fig. 9 .
Fig. 9. Detail, from Figure 8, of the IBP and GH upper and lower boundaries for the last 3 Myr before the Present.

Fig. 9 .
Fig. 9. Details from Fig. 8 of the IBP and GH upper and lower boundaries for the last 3 Myr before the present.

Fig. 10 .
Fig. 10.An illustration of the sensitivity of current Case 1 model results t uncertainties in the earlier GST history, prior to the 100 ka cycles.The m comparison of a -warmer‖ GST model (light green) to a -colder‖ GST m 14 Myr record of global marine isotopic data adjusted to the Mallik site.model is up to 6 K higher for most of the period between 3 Myr and 0.9 M difference between these two models is tens of meters, at most.It indicat youngest part of the model is not very sensitive to the early parts of the G

Fig. 10 .
Fig. 10.An illustration of the sensitivity of current Case 1 model results to potential uncertainties in the earlier GST history, prior to the 100 ka cycles.The model comparison of a "warmer" GST model (light green) to a "colder" GST model based on 14 Myr record of global marine isotopic data adjusted to the Mallik site.The "warmer" model is up to 6 K higher for most of the period between 3 Myr and 0.9 Myr ago.The current difference between these two models is tens of meters, at most.It indicates that the youngest part of the model is not very sensitive to the early parts of the GST history.

Fig. 11 .Fig. 11 .
Fig. 11.The onshore and offshore (the Holocene marine transgression) GST models of the present glacial cycle and the future warming in its end.The Holocene starts at time 550 ka of the simulations and the present time is 563.5 ka. .Fig. 11.The onshore and offshore (the Holocene marine transgression) GST models of the present glacial cycle and the future warming in its end.The Holocene starts at time 550 ka of the simulations and the present time is 563.5 ka. 35

Fig. 13 .
Fig. 13 .Tempertaure (T)-depth (z) profiles corresponding to a model based on Taylor's surface temperature history affected by both the Holocene marine transgression and the possible future warming due to a doubling of atmospheric CO2 content.

Fig. 13 .
Fig. 13.Tempertaure (T )-depth (z) profiles corresponding to a model based on Taylor's surface temperature history affected by both the Holocene marine transgression and the possible future warming due to a doubling of atmospheric CO 2 content.
Fig. 8.A comparison of Case 1 and Case 2 results illustrating the respective depth variations of IBP and GH upper and lower boundaries during the last 14 Myr.

Table A1 .
Nomenclature.A Rate of heat generation [W m −3 ] C v Volumetric heat capacity [J/(m 3 K −1 )]where transgressed by Holocene sea level rise.As well, where IBP is sufficiently thick, both are expected to persist into the next glacial even if atmospheric CO 2 doubling increases GST warming. even