Modeling of stability of gas hydrates under permafrost in an environment of surface climatic change – terrestrial case , Beaufort-Mackenzie basin , Canada

Introduction Conclusions References


Introduction
The attempt to model the onset of ice bearing permafrost (IBP) formation, the formation of petroleum gas hydrate (GH) and modeling for the more detailed Holocene temperature history requires knowledge of the surface temperature forcing environment.This needs to be derived from sources sometimes distant from the study area of the Canadian Arctic and specifically the Beaufort Mackenzie Basin (BMB) area where very thick GHs are located currently.
Recent results of the analysis of the Holocene-Pleistocene temperature history using 600 k yr by Majorowicz et al. (2008) showed that GH can be stable through interglacial intervals despite large variations in surface temperatures due to the buffering effect of the IBP and the retardation of GH dissociation due to latent heat effects.The present paper focuses on the analysis of the origin, growth and persistence of the GH and IBP layers.Understanding of formation and persistence of sub-permafrost GH accumulations should help us to determine if a geological natural gas resource model can be formulated.Our analysis will use the characteristics of IBP and GH response to variations of surface temperature history, to simultaneously constrain the models of past environments that led to the initiation, growth and persistence of the linked occurrences of both IBP and GH accumulations in the subsurface of the Canadian Arctic.The work will also address the risks posed by global and regional temperature change, as well as will provide a tool that will assist the appraisal of risks related to surface imposed changes on the GH layer.In order to model conditions for the IBP formation due to cooling of the surface and following GH formation considering latent heat "delay" effects (due to deeper formation conditions) we need to know surface conditions in the past.
To understand better the formation and history of petroleum GHs in terrestrial IBP regions we have performed numerical modeling of the surface forcing due to general cooling trend started in late Miocene and more detailed glacial-interglacial history.Persistent GH layers in a terrestrial environment of thick IBP in cold regions sequester methane and impede its migration into the atmosphere.The Mallik site in the Mackenzie Delta is an excellent example of such GH deposits situated in the large area of deep GH Type I stability in the onshore and offshore (Fig. 1) under continental and continental relic IBP (accumulated during the Pleistocene marine low-stand prior to the Holocene sea level rise).More detailed descriptions can be found in a collection of works edited by Dallimore and Collett (2005), Judge and Majorowicz (1992), Majorowicz and Hannigan (2000) and Majorowicz and Smith (1999).These works describe the geothermal environment of the BMB hydrates and Mallik site in particular.We model the onset of permafrost formation and subsequent GH formation in the changing surface temperature environment for the Beaufort-Mackenzie (BMB), Fig. 1.We also model terrestrial BMB GH thickness variations below an ice-bearing IBP layer in response to past surface temperature changes onshore and offshore in BMB using a 1-D thermal model that assumes no water or gas flow.

Assumed Surface temperature history
Our evidence for past temperatures comes mainly from isotopic considerations (especially δ 18 O).Our models of surface forcing of temperature are based on several sources (Muller and MacDonald, 2000;Frakes et al., 1992;Taylor et al., 2005;Intergovernmental Panel on Climate Change, 2007).
The analysis of the paleoclimatic data for the Phanerozoic and specifically for the last tens of millions years (Myr) from the above sources and ones in the References show: Sixty-five to one hundred Myr ago average global temperatures were the highest during the last ∼200 Myr (Wikipedia.org, 2008).This obviously prevented the formation of large scale ice sheet.Forests extended all the way to the poles.Temperature also peaked 50-52 Myr ago, when atmospheric CO 2 ranged 1125-3000 parts per million.By 20 Myr ago, CO 2 dropped to about 400 ppm (Polsson, 2007;Lawrence, 2006;Wikipedia.org, 2008).
Following the Paleocene to early Eocene peak warming the climate cooled variably towards the Pleistocene glacial environment (Wikipedia.org, 2008).However, 3 to 6 Myr ago, globally averaged temperatures (Wikipedia.org, 2008) (Fig. 2) were still higher than today with surface temperature at poles higher than present temperatures.The Northern Hemisphere likely had no continental glaciers [Muller and MacDonald, 2000;Frakes et al., 1992;Lawrence, 2006].
Climate during the following 3 Myr before present changed dramatically in response to astronomical effects (Milankovitch cycles) caused changes in surface forcing.These Figures

Back Close
Full resulted in cycles of glacials and interglacials within a gradually deepening ice age.
The growth and retreat of continental ice sheets in the Northern Hemisphere occurred at 2 frequencies; 41 kyr and 100 kyr (Fig. 2).The gradual intensification of this ice age over the last 3 Myr has been associated with declining concentrations of the CO 2 which partly influences this change in temperatures.El Nino was continual rather than intermittent until 3 million yr (Myr) ago.Significant amplification of the response of climate to orbital forcing began 3 Myr ago, resulting in drastic oscillations between ice ages and warmer periods over the past 1 Myr (Fig. 2).

P-T time envelope and formation of permafrost and gas hydrate -preliminary considerations
In order to simulate the IBP and GH formation, we have constructed the general model of the surface temperature forcing history for the last 14 Myr (Figs. 4-5) based on temperature history given in (Muller and MacDonald, 2002;Frakes et al., 1992;Wikipedia.org, 2003).The most recent surface climate history for the end of the Wisconsinan and Holocene is after Taylor et al. (2005).If the equivalent Vostok ∆T temperature of 0  formation.The equilibrium IBP thickness is some 250 m for a GST of −5 • C.This means that GH could have formed immediately below the IBP base.Because the geotherm and the phase curve are nearly parallel in this range of p-T, even a small downward shift of the geotherm means a large broadening of the zone with the p-T conditions favorable for GH formation.
Further surface cooling below −5 • C caused the IBP base to descend into the layer already occupied by GH.For 30 % porosity and 60 % GH saturation, the remaining water available to freeze occupied not 30 %, but only 12 % of the whole volume and the latent heat released by its freezing was 2.5 times smaller than in the IBP interval above the GH.The question is, "what had happened with the methane in the uppermost 300 m, if any gas was initially collected under the IBP?" where T is the temperature, K 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.We employed a computer code to simulated temporal subsurface temperature changes in response to surface forcing by Safanda et al. (2004).Within the model Eq. ( 1) is solved numerically by an implicit finite-difference method similar to 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 kth depth step and the nth time step, respectively; This difference scheme, together with the upper and lower boundary conditions lead to a system of difference equations for unknown values T fractions Galushkin (1997) and Nixon (1986).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 effects of interstitial ice and GH were accounted for using apparent heat capacity according to Carslaw and Jaeger (1959), when the volumetric heat capacity is increased in the depth sections of the model where the thawing and freezing occurs, i.e.where the temperature is within the thawing rangebetween the temperature of solidus T S , and liquidus, T L , at the actual simulation time step.
The liquidus and solidus temperatures of water/ice and GH are depth and hydrostatic pressure dependent (Galushkin, 1997) and solidus temperatures were 0.2 • C lower than liquidus temperatures.A contribution to the heat capacity from the latent heat = ρΦL/(T L − T S ) was considered, where ρ is the density of either ice or GH and Φ is a fraction of the total volume occupied by these phases.In the IBP zone we infer the 30 % 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 %.The salt concentration 9 g L −1 was considered constant with depth and the p-T phase curves were adjusted to this value.
Numerical code performance was tested by comparing model 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 degree C. A similar error range was estimated by halving the time and/or depth steps.according to Sloan (1998), (Fig. 3).Previously published models have considered only a thin layer using a constant dissociation temperature Taylor (2005).
In simulating the temporal subsurface changes in models with vertically separated occurrence of IBP and GH, a consumption/release of the latent heat during decay/formation of IBP or GH was considered separately.This means that in 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 in the model explicitly.
The above conditions may not be appropriate in models of possible simultaneous IBP and GH occurrence.In such cases it is prescribed at each depth point, which fraction of the pore space can be occupied by water ice and which by GH.Then the code checks independently the p-T conditions for IBP and for GH in each depth point, in each time step and adjusts the volumetric heat capacity to the actual temperature at a given depth Carslaw and Jaeger (1959).For example, in the Mallik case, where a porosity of 30 % was assumed, the prescribed pore fraction of GH in the uppermost 250 m was zero and that of IBP 100 %, because as the preliminary calculations had shown, IBP forms prior to hydrate.On the contrary, below 250 m the GH forms first.For the 60 % saturation of the pore space by GH this means that the ratio of pore fractions occupied by GH and IBP is 60: 40 at this depth.

Model for the 14 Myr history of surface temperature change -results
We have simulated the downward propagation of the surface warming and cooling considering the cyclical glacial and interglacial models for the eastern Richards Island location where Malik well was drilled for hydrate (Fig. 1).
The dependence of the thermal conductivity on water/ice content and the specific heat of the rock section on the porosity and the proportion of interstitial water and iceare important.Accounting for the effect of the latent heat necessary to thaw the interstitial ice in the IBP layer, is crucial for matching observations at realistic time rates.In the absence of this heat sink provided by thawing ice in the IBP, the subsurface warming 2871 Figures

Back Close
Full would proceed much faster.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.
All models account for latent heat by means of the apparent specific heat, which is a standard treatment.The model also considers diffusive heat flow related to surfacesubsurface coupling.
There was a warm period (25 to 15 Myr ago) with maybe no IBP in Mallik.According to our model, the surface temperature in Mallik dropped below 0 • C 14 Myr ago, which initiated the IBP formation.We have used a stepwise approximation of the surface temperature with a length of the step 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 cycle and from the 0.9 Myr the 100 ka cycle.

Case 1 -Gas hydrate formation controlled by the geological gas entrapment conditions
In Case 1 (Fig. 4) we allow occurrence of GH only in the geologically confined area under a geologically impermeable seal.This is a probable situation for case of the Beaufort-Mackenzie (BMB) GHs as shown by stratigraphic sections in the Mallik boreholes Majorowicz et al. (2008).At Mallik, the interbedded character of the GH bearing strata indicates a lithologic control on GH occurrence.GH layers occur in coarse-grained sandstone beds separated by thin non-hydrate-bearing, fine-grained siltstone and claystone beds.The gas and GH at Mallik appears to be entrapped in association with the diastrophic structure.The fault bounding the structure is the likely conduit for migration of the gas into the GH stability zone (GHSZ).Methane isotopic compositions of 12 GH samples averaged −42.7 % and clearly indicate a thermogenic source according to Lorenson et al. (2004).
Thermogenic gas likely migrated up along listric-normal growth faults from larger depths where gas generation is possible due to availability of the source rock and 2872 Figures

Back Close
Full thermal conditions (thermogenic gas generation envelope temperatures are likely in 4-5 km in the area of study).The deep upward migrating gas could have been trapped in the anticline and or tilted fault blocks and turned into GH when cooling took place.Such a hypothesis has some grounds in the coincidence of GH occurrences above the deep known hydrocarbon reserves as estimated by Majorowicz and Osadetz (2001).
We relate our model to the above cited findings.The division into IBP and GH zones was prescribed in the model explicitly.We have vertically separated the occurrence of IBP and GH.A consumption/release of latent heat during the decay/formation of IBP or GH is considered separately.It means that in IBP zone, only the latent heat and p-T phase curve of water ice were taken into account, and in the GH zone it was only latent heat and p-T phase curve of GH.
GH formation was considered below 900 m only.The GH did not form even during the coldest 41 ka cycles (from −11.5 • C to −7.Despite our use of quite general surface temperature history (Fig. 4), agreement between the model and the observation ofthe present-day state of GH and IBP is reasonably good.
The plot of the whole 14 Myr history of surface forcing effects on IBP and GH is shown in Fig. 5 Case 2 -considering simultaneous occurrence of permafrost and gas hydrate Case 2 deals with situation when there is no impermeable seal in the sedimentary formation and gas can flow freely.The simultaneous occurrence of IBP and GH is allowed.This has some geological grounds as the geological seals commonly seen in the BMB area are not widely available Majorowicz et al. (2008).
The numerical simulation of the subsurface temperature response to changes of the surface temperature forcing was done for the last 14 Myr (Figs. 6-9).The GST history is the same as in Case 1 model shown in Fig. 5.
In the times of warm climate, like before 14 Myr ago, when GST at the Mallik site was +1 • C, the subsurface temperatures were above p-T phase curves of both IBP and Figure 6 shows transient temperature-depth profiles corresponding to times 20, 10 000, 20 000 and 42 500 yr after the GST cooling from −4.5 • C to −5.5 • C 6 Myr ago, and to time 20 yr after the consequent warming back to −4.5 • C at 5.5 Myr ago,all that preceded by the whole GST history since 14 Myr ago.After the cooling at 6 Myr ago, the T-z profile started to move downward and touched the p-T phase curve of the GH at the depth interval 290-300 m 42.5 ka after the GST drop (light blue profile in Fig. 6 for time 5.9575 Myr ago).At that moment, the IBP base was at the depth of 250 m.It means that there was an unfrozen depth section (250-290 m) separating IBP from hydrate.With further subsurface cooling, the zone occupied by GH propagated not only downward, but also upward and after another 100 ka the upper GH boundary touched the IBP base at the depth of 255 m.Upward spreading of the GH zone stopped at this moment, because all pore space in overlying IBP was occupied by frozen water.Nevertheless, this contact of IBP with GH, which stopped upward growth of the GH layer, does not stop downward migration of the IBP layer, because the model considers 60 % saturation of the pore space by hydrate.The rest of 40 % of the pore space (12 % of the whole volume) was occupied by water, which could have frozen with subsequent cooling.This scenario and consequences following from it for IBP and GH occurrence were taken into account by considering separately the porosity available for water ice and that for the hydrate.For the ice, it was 30 % above the line of the IBP-GH contact at approximately 250 m and 12 % below it.Similarly, the porosity available for the GH above the line of contact at 250 m was considered to be 0 % and below it 18 % (60 % saturation in 30 % porosity rock).This opened the possibility to treat the latent heat released/consumed in the solidus-liquidus zones of IBP (0.334 MJ kg −1 ) and GH Depth variations of the IBP and GH upper and lower boundaries during the last 14 Myr are depicted in Fig. 8 and in detail for the last 3 Myr in Fig. 9. Also shown are variations of these parameters for the Case 1 Model, where the GH formation was restricted to the depth below 900 m only, and simultaneous occurrence of IBP and GH was not possible.The shallower position of both IBP and GH in this model is caused by lower conductivity and consequently higher temperature gradient within the IBP in the zone of its simultaneous occurrence with GH below 250 m.It is also evident from Fig. 9 that depth variations of both the IBP and GH basis are slightly larger in Case 2 than Case 1 during the 41 ka and 100 ka cycles.The most probable reason for the larger variations of the IBP base is a smaller damping effect of the latent heat release/consumption at the base due to a smaller amount of freezing/melting water/ice.Most of the pore water is bounded in the GH, which is stable in this depth range.The smaller damping effect at the IBP base also means that the subsurface temperature changes propagate faster to the GH base and that these cause its larger depth variations compared to Model -Case 1.
The Case 2 model assuming presence of methane in the whole rock, or at least below 250 m, before the onset of cooling 14 Myr ago predicts the following subsurface temperature-GH and IBP changes resulting from the gradual climate cooling since this time: Introduction

Conclusions References
Tables Figures

Back Close
Full the GH survived all the following GST variations.Its upper constraint at 250 m was always within the GH stability zone and its maximum extent to 990 m was reached during the most recent glacials.
5. Shortly after the definitive GH formation, the downward migrating IBP base and upward migrating upper GH boundary met at the depth of about 250 m.Since that moment, IBP below 250 m could have "used" only that part of the pore water, which was not bounded in hydrate.
The Case 2 model provides no mechanism enabling methane to escape to the earth surface from below the IBP lid formed 14 Myr ago.Its forecast of simultaneous occurrence of IBP and GH below 250 m is, however, not generally observed Majorowicz et al. (2008).
The simultaneous occurrence also implies lower conductivity in this zone, which means higher temperature gradient, which shifts the IBP and GH basis upward, above 2877 Introduction

Conclusions References
Tables Figures

Back Close
Full their position observed currently.The indicated upper GH boundary above 250 m shown in Figs.8-9 is a formal one.It was calculated according to the p-T phase curve.In reality, no GH was considered there because all pore space above 250 m is occupied by water ice and GH cannot form in a region already occupied by ice.

Discussion and conclusions
Historical surface temperature forcing as well as geological control have significant implications for both IBP and GHs model results that consider latent heat effects of water/ice and GH formation and dissipation and 14Myr GST history.
Our model, in which we considered a stepwise approximation of the GST with a length of the step 1 Myr between 14 and 6 Myr ago and 0.5 Myr between 6 and 2.5 Myr followed from 2.5 Myr ago by the 41 ka cycles and from 0.9 Myr by the 100 ka cycles, shows that the onset of GH formation always comes after IBP formation begins.IBP, provides a potentially impermeable seal for any upward gas migration and this is inferred to have existed long before the p-T conditions enabled formation of the first GH.
The models based on the two scenarios proposed, i.e; Case 1: formation of GH from natural gas previously entrapped under a deep geological seal Case 2: formation of GH from gas in all of the free pore space simultaneous with IBP formation; indicate two entirely different GH generation onset times for depending on model assumptions.For Case 1: GHs could have formed at 0.9 km only some 0.9 Myr ago.For Case 2: the first GH formed in the depth range 290-300 m shortly after 6 Myr ago when Figures

Back Close
Full the GST dropped from −4.5 • C to −5.5.• C and started to spread both downward and upward.
In case 2, after all the pore water had frozen, methane contained in the freezing rock could be prevented from escaping to the surface by a IBP seal and this methane might have been pushed downward below the downward migrating IBP base.It is also possible that this methane could have escaped to the sides of the IBP seal and then migrated upward.
Detailed models of recent ∼100 kyr cycles in the BMB area shows that GH layer thickness generally increases during colder intervals (i.e.glacial) and decreases during warmer intervals (i.e.interglacial).In the more recent glacial history of the ∼100 ka cycles these variations are ∼0.2 km for the IBP and about 0.1 km for the GH layer.Where the IBP layer is thick it is unlikely that sub-IBP GHs disappeared entirely during previous interglacial intervals, nor are they expected to disappear prior to the 'natural' end of the current interglacial.In regions of thick terrestrial IBP like the Mackenzie Delta, GH layers can act as a persistent sink for methane from deep thermogenic sources and petroleum systems, and as a barrier to the migration of methane into the atmosphere.
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 applies mainly to marine non-IBP GHs (Nisbet, 1990(Nisbet, , 2002)), which may be more easily destabilized than are the terrestrial sub-IBP GHs we modeled.Our study shows that sub-IBP GHs below thick IBP vary in thickness in response to surface temperature history changes, but that terrestrial thermal inertia conserves both IBP and sub-IBP GHs delaying and reducing methane release.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 could rapidly reinforce climate warming events, consistent with the hypothesis.The implications of latent heat effects and thermal inertia for submarine GHs remain to be determined, however, our model results appear Figures consistent with the implications of recent observations of methane isotopic compositions from ice cores for the "clathrate gun" hypothesis of Sowers (2006).
This study shows that thermal modeling is essential for the understanding of the origin, growth and persistence of sub-IBP GH accumulations, as a function of temperature history.It can help to determine and distinguish the features of the GH geological natural gas resource model, which have both climate and resource implications.
a 1-D model that separated IBP and GH layers in an onshore setting during the last Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | C indicated for the present time (Fig. 2) corresponds to the present ground surface temperature of −6 • C at Mallik Richards I. Beaufort-Mackenzie basin, Canada (Fig. 1), then during the period 6-3.5 Myr ago it corresponds to ground surface temperature (GST) oscillations at Mallik between −6 • C and −4 • C. The range of GSTs between −6 • C and −4 • C is critical for GH formation as shown in preliminary test for the steady-state in Fig. 3.The geotherm coresponding to the GST of −5 • C (Fig. 3) approaches the GH phase curve from above at the depth of about 300 m.According to it, the GH could not have formed when the GST was higher than −5 • C. Before 6 Myr ago the climate was even warmer, which means that the first GH deposits might have formed around 6 Myr or later, compared to the onset time of 14 Myr ago for the IBP (Fig. 4).It means that the IBP, a potentially impermeable seal for any potential gas migrating upward, existed for a long time before p-T conditions permitted the first hydrate Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Solving the transient heat conduction equation gives the temporally dependent subsurface temperature change in response to surface forcing: tridiagonal matrix, which was solved by the forward method ofPeaceman and Rachford (1955).To estimate effective thermal conductivity values and volumetric heat capacity, 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 Discussion Paper | Discussion Paper | Discussion Paper | Our model uses deep heat flow, thermal conductivity, present IBP, Type I GH thicknesses and a surface melting temperature (−0.576 • C) that considers groundwater salinities (9 g L −1 ).It employs latent heat effects throughout the IBP and GH layers, which improves upon previous models Taylor (2005).The models are constrained by deep heat flow from bottom hole temperatures determined from deep wells Majorowicz et al. (1990), Taylor et al. (1982), and thermal conductivity, latent heat, present IBP thickness and observed present Type I GH thicknesses Henniges et al. (2005) and 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 IBP and GH layers Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5• C), because the steady-state geotherm corresponding to the mean surface temperature of the cycle, −9.5 • C, crosses the GH stability curve just above 900 m (see Fig.3) and duration of the cold phase of the cycle with surface temperature of −11.5 • C is too short (20.5 ka) for the temperature to cool close to the steady-state curve corresponding to −11.5• C. Results of the simulations (Fig.4) indicate that GH started to form only with the onset of the 100 ka cycles 0.9 Myr ago, when the surface temperature of the glacial phase dropped to −15 • C. The simulated present-day bases of IBP, 541 m, and GH, 1060 m, are slightly above its observed positions in Mallik, 600 m and 1107 m, respectively.In our previous simulations Majorowicz et al. (2008), the bases appeared deeper, at 600 m and 1165 m.The main reason is probably the shorter duration of the cold phase of the glacial, (75 ka instead of 90 ka), considered in the latest simulations.The length of the whole glacial cycle circa 100 ka, is consistent with other sources and the length of the interglacial is 25 ka.The surface temperatures during the glacial, −15 • C, and the interglacial were the same as that of previous model 3 by Majorowicz et al. (2008) modified from Taylor et al. (2005, Fig. 3b) Taylor (2005).Discussion Paper | Discussion Paper | Discussion Paper | GH at all depths (for the geothermal model with basal heat flow of 60 mW m −2 and conductivity of 2.1 W (m K) −1 ).With the onset of a gradual cooling at 14 Myr ago, when the GST dropped below the solidus temperature of IBP, the IBP started to propagate from the surface downward.However, for the model considered (30 % porosity), with all the pore water frozen, thermal conductivity of the frozen rock 3.4 W (m K) −1 and the latent heat of water freezing 0.334 MJ kg −1 ) the temperatures stay above the p-T phase curve of GH until the GST drops to −5.5 • C. A GST of −5 • C was not low enough for GH formation in this geothermal model.It is questionable, what could have happened with the methane contained in the freezing rock.If all pore water had frozen, it could not have escaped to the surface and might have been pushed downward below the downward migrating IBP base or escaped to the sides of the IBP body and then migrated upward.Such IBP "weak" thermokarst areas are under water bodies of the present Mackenzie River Delta.Discussion Paper | Discussion Paper | Discussion Paper |

( 0 .
430 MJ kg −1 ) independently and appropriately.The thermal conductivity in the zone of simultaneous occurrence of IBP and GH was estimated at 2.7 W (m K) −1 .It is lower value than that for pure IBP (3.4 W (m K) −1 because of the presence of GH with a low conductivity 0.45 W (m K) −1 at the expense of water ice with conductivity 2.2 W (m K) Discussion Paper | Discussion Paper | Discussion Paper | to form at time 5.9575 Myr ago reached its maximum thickness 125 m (between 255-380 m) shortly after the next GST change from −5.5 • C to −4.5 • C, which occurred 5.5 Myr ago (Fig. 7).The same is true with the IBP, which penetrated into the GH only 4 m (255-259 m).As can be seen in Fig. 7, following that GST warming, the last GH disappeared in the depth range 310-340 m at time 5.472 Myr ago when it was separated from the IBP base at the depth of 240 m by a 70 m thick melted zone.The GH formed again shortly after the next decrease of the GST from −5.0 • C to −5.5 • C 4.5 Myr ago.Since this moment, the GH survived all the following GST variations.Its upper constraint at 250 m was always within the GH stability zone and its maximum extent to 990 m was reached during the most recent glacials.
Discussion Paper | Discussion Paper | Discussion Paper | 1. IBP thickness did not exceed 200 m before the time 6 Myr ago 2. GST temperatures higher than −5.5 • C prevailing before the time 6 Myr ago did not enable the GH to form 3. First GH formed in the depth range 290-300 m shortly after the time 6 Myr ago when the GST dropped from −4.5 • C to −5.5.• C and started to spread both downward and upward.After another 100 ka, the upper GH boundary touched the IBP base propagating downward at the depth of 255 m 4. First GH disappeared after subsequent warming and formed definitively shortly after GST drop from −5 • C to −5.5 • C at time 4.5 Myr ago.Since this moment, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tial equations, J. Soc.Ind. Appl.Math., 3, 28-41, 1955.Polsson, K.: Climate Change -Carbon dioxide, http://www.islandnet.com/∼ kpolsson/climate/ carbondioxide.htm,2007.Safanda, J., Szewczyk, J., and Majorowicz, J. A.: Geothermal evidence of very low glacial temperatures on the rim of the Fennoscandian ice sheet, Geophys.Res.Lett., 31, Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig.1.Depth to GH base of Type I stability in the Mackenzie Delta (modified from Majorowicz & Hannigan 2000 [4]).Modelling undertaken here is for northern Richards Island and offshore (near centre of figure).

Fig. 1 .Fig. 1 .
Fig. 1.Depth to GH base of Type I stability in the Mackenzie Delta (modified from Majorowicz and Hannigan, 2000).Modelling undertaken here is for northern Richards Island and offshore (near centre of figure).

Fig. 4 .Fig. 5 .
Fig. 4. Surface temperature and bases of the IBP and the GH in the last 3 Myr according to numerical simulation of the last 14 Myr subsurface temperature forcing -Case 1 model.

Fig. 6 .
Fig.6.Transient temperature-depth profiles as a response to the GST cooling from -4.5 °C to -5.5 °C 6 Myr ago (and to all previous GST history since 14 Myr ago) followed by warming back to -4.5 °C 5.5 Myr ago.

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

Fig. 7 .
Fig. 7. Transient temperature-depth profiles as a response to the sudden GST warming from -5.5 °C to -4.5 °C 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.Curve for time 28 000 years (5.472Myr ago) shows the moment when the last GH disappears in the depth interval 310 -340 m, separated from the IBP base at the depth of 240 m by a 70 m thick melted zone.

Fig. 7 .Fig. 8 .
Fig. 7. Transient temperature-depth profiles as a response to the sudden GST warming from −5.5 • C to −4.5 • C 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 yr after the warming.Curve for time 28 000 yr (5.472 Myr ago) shows the moment when the last GH disappears in the depth interval 310-340 m, separated from the IBP base at the depth of 240 m by a 70 m thick melted zone.

Fig. 8 .Fig. 9 .
Fig. 8. Depth variations of the IBP and GH upper and lower boundaries during the last 14 Myr for Case 2 compared with results of Case 1, where the GH formation was restricted to the depth below 900 m.

Fig. 9 .
Fig. 9. Detailed part of depth variations from Fig. 9 for the IBP and GH upper and lower boundaries for the last 3 Myr before present.