Abstract
Gas hydrates (GHs) are a prominent subsurface feature on the Canadian Arctic continental margin. They occur both onshore and offshore, although they formed generally terrestrially, during the last glacial sea level low-stand, both in a region that was persistently glaciated (Queen Elizabeth Islands Group, Canadian Arctic Archipelago (QEIG)), and in a region that was not persistently glaciated (Mackenzie Delta-Beaufort Sea (MD-BS)). Parts of both regions were transgressed in the Holocene. We study the dynamic permafrost and GH history in both regions using a numerical model to illustrate how changes in setting and environment, especially periodic glacial ice cover, affected GH stability. MD-BS models represent the Mallik wellsite and these models successfully match current permafrost and GH bases observed in the well-studied Mallik wells. The MD-BS models show clearly that GHs have persisted through interglacial episodes. Lower surface temperatures in the more northerly QEIG result in an earlier appearance of GH stability that persists through glacial-interglacial intervals, although the base of GH base stability varies up to 0.2 km during the 100 ka cycles. Because of the persistent glacial ice cover QEIG models illustrate pressure effects attributed to regional ice sheet loading on the bases of both permafrost and GHs since 0.9 MYBP. QEIG model permafrost and GH depths are 572 m and 1072 m, respectively, which is like that observed commonly on well logs in the QEIG. In order to match the observed GH bases in the QEIG it is necessary to introduce ice buildup and thaw gradually during the glacials and interglacials. QEIG sea level rose 100–120 m about 10 ka ago following the most recent glaciation. Shorelines have risen subsequently due to isostatic glacial unloading. Detailed recent history modeling in QEIG coastal regions, where surface temperatures have changed from near zero in the offshore to −20°C in the onshore setting results in a model GH stability base, that is, <0.5 km. These coastal model results are significantly shallower than the inferred average GH base about 1 km in wells, Smith and Judge (1993). QEIG interisland channels are generally shallow and much of the previous shoreline inundated by the Holocene transgression was above the glacial sea level low-stand during the last ice age, resulting in a QEIG setting somewhat analogous to the relict terrestrial GH now transgressed by the shallow Beaufort Sea. It is also possible that the marine conditions were present at emergent shorelines for a shorter time or that the pretransgression subsurface temperatures persisted or were influenced by coastal settings, especially where lateral effects may not be well represented by 1D models.
1. Introduction
Natural petroleum gas hydrates (GHs) are a significant feature of Canada’s continental margins and Arctic regions [1]. Two major types of gas hydrate occurrences are found on Canada’s continental margins at a large range of pressure-temperature stability conditions [1, 2]. The first are essentially terrestrial GHs that formed in Arctic regions on, or adjacent to the, continental shelf in response to low ground surface temperatures during glacial intervals in the Mackenzie Delta-Beaufort Sea (MD-BS), [3] and Queen Elizabeth Islands Group (QEIG) of the Arctic Archipelago [1], (Figure 1, areas A and B). Part of both areas were transgressed in response to Holocene sea level rise and, additionally, part of the QEIG coastline has emerged more recently as a result of isostatic uplift following ice sheet removal [4, 5]. The second type occurs on sub-Arctic continental margins, along the Pacific and Atlantic-Baffin continental shelves, where water column pressure and temperature contribute to gas hydrate stability [6–8]. The second type of GHs probably exist on the margins of the Arctic Ocean Basin, but a lack of data, other than a “BSR” from the Alpha Ridge [8] precludes current description of much of the Arctic continental margin.
Majorowicz and Osadetz [1] estimated a conservative volume of 1010–1012 m3 GHs with an associated methane gas potential in the range of 1012–1014 m3 and pointed out that large part of gas hydrate occurrences is related to deeper hydrocarbons found in the Mackenzie Delta-Beaufort Sea and Sverdrup basin (Arctic Islands). Lesser known are the biogenic deposits.
The volume of methane in hydrates in Canada is geographically distributed in the following regions (Figure 1):(A)the Mackenzie Delta—Beaufort Sea,(B)the Arctic Archipelago,(C)the Atlantic margin, and(D)the Pacific margin.
The total in situ amount of methane in Canadian hydrates is estimated to be 0.44– m3, as compared to a conventional Canadian in-situ hydrocarbon gas potential of approximately m3. Osadetz and Chen [9] reevaluated Mackenzie Delta-Beaufort Sea gas hydrate resources as a function of the spatial variation of pore space gas hydrate saturation (Sgh). Deterministically they mapped m3 raw initial natural gas in place (GIP), where, if the average Sgh is either >30% or >50%, the inferred GH resource volumes are m3 to m3 GIP, respectively. They also provided a probabilistic resource appraisal with an expected total resource 10.2 m3 GIP, which similarly was studied as a function of Sgh. If the average Sgh is either >30% or >50%, then the respective GH resource volumes inferred probabilistically are m3 and m3 GIP, expected, respectively, a result comparable to the mapped deterministic estimates, and not greatly different from the regional estimate of Majorowicz and Osadetz [1].
In this study we focus on GH occurrences that are at or near the Canadian margin of the Amerasian Basin of the Arctic Ocean, its marginal seaways in the Canadian Arctic Archipelago, including the Queen Elizabeth Islands Group (QEIG), and the contiguous onshore areas (Figure 1, areas A and B). GHs occur both onshore and offshore in both regions studied. The primary difference between the two regions is that the QEIG was persistently and thickly glaciated as indicated by a strong postglacial isostatic uplift record [4, 5]; while, the MD-BS is not inferred to have been either persistently or thickly glaciated and it lacks significant postglacial isostatic uplift [5]. The QEIG is completely underlain by the upper Paleozoic to Tertiary Sverdrup Basin succession [10, 11] while the MD-BS is underlain by the largely Cretaceous and Tertiary successions of the Beaufort-Mackenzie Basin [3, 12].
Using a 1D numerical model we model geothermal regime and infer how these might affect gas hydrate stability and ice-bonded permafrost regime and distribution in both regions. Gas hydrate and ice-bonded permafrost model studies require knowledge or inference of the surface temperature history.
Majorowicz, Osadetz, and Safanda [13] developed an effective thermal history model for the analysis of GH and permafrost stability in the onshore Beaufort-Mackenzie basin constrained by Mallik project results [3] (Figure 1 area A). They showed that the formation and history of thermogenic petroleum GHs in permafrost regions can be modeled numerically using a model of surface temperature forcing due to both glacial-interglacial history and future climate change. Their model employed both a detailed Holocene history and a future climate where atmospheric CO2 has doubled, due to the warming inferred to accompany current climate change trends. The initial model runs were constrained by deep heat flow estimates consistent with bottom-hole temperatures in deep wells and observed type I GH and permafrost thicknesses and characteristics from Mallik wells.
The models considered the pressure—depth dependence of ice and GH thawing points over the entire modeled extent of both GH and permafrost layers. Initial model results, in areas of thick permafrost, showed that expected warming in neither onshore nor offshore Beaufort Mackenzie Basin will destabilize GHs completely. In addition, where the permafrost is thick, a thinned GH stability zone is inferred to have persisted through previous interglacial periods. Previously Taylor and others [14, 15] modeled permafrost and GH stability in the present and near future environment, including both the recent glacial-postglacial history and projected future climate change for both the MD-BS (area A in Figure 1) and the offshore Pacific margin.
In this work we use the same numerical modeling approached applied by Majorowicz, Osadetz, and Safanda [13] to the MD-BS, but applying it to model the characteristics of permafrost and GH origin, growth, and persistence, as a function of temperature history and ice sheet loading (where applicable) during the last 3 million years in order to compare and contrast cases, considering specifically: (1) MD-BS (Figure 1()) in the vicinity of Richards Island and the Mallik wells, an area that is inferred not to have been persistently or thickly glaciated and which lacks significant post-glacial isostatic uplift [5], versus (2) QEIG (Figure 1()) in Sverdrup Basin, which was persistently and thickly glaciated and where there is a strong postglacial isostatic uplift record [4, 5] such that the glacial ice sheet thermal and pressure effects were significant.
2. Method
Solution of the transient heat conduction equation gives the temporally dependent subsurface temperature change in response to surface temperature forcing: where is the temperature, K is the thermal conductivity, is the volumetric heat capacity, A is the rate of radiogenic heat generation per unit volume, z is the depth, and 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. [16]. Within, model (1) is solved numerically by an implicit finite-difference method [17]. The upper boundary condition is the temporally varying ground surface temperature (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–10000, 10000–15000 m, respectively. Time steps vary between 0.5 to 50 yrs, depending on the amplitude of GST changes. The finite-difference scheme of (1) on the depth and time grids and , respectively, has a form: where the subscript and the superscript denote values at the kth depth step and the 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 , , (where subscript k and superscript n denote similarly the k-th depth and the th time steps) within a tridiagonal matrix, which was solved by the forward method of Peaceman and Rachford [18].
To estimate effective thermal conductivity values and volumetric heat capacity, we consider the respective geometric and arithmetic averages of the constituent values for the rock matrix, water, ice, and GH in proportion to their fractional volumes [17]. A consumption or release of latent heat, L, accompanying either the water/ice (334 kJkg−1) or GH (430 kJkg−1) transitions, accompanying thawing or freezing was included. The effects of interstitial ice and GH were accounted for using apparent heat capacity [19], when the volumetric heat capacity is increased in the model depth sections where the thawing and freezing occurs, that is, where the temperature is within the thawing range between the temperature of solidus, ,and liquidus, , during the actual simulation time step.
The liquidus and solidus temperatures of water/ice and GH are depth and hydrostatic pressure dependent [17] and solidus temperatures were 0.2°C lower than liquidus temperatures. A contribution to the heat capacity from the latent heat = ) is considered, where ρ is the density of either ice or GH and Φ is a volume fraction occupied each phase. In the ice-bonded permafrost (IBP) zone we infer the 30% rock matrix porosity to be occupied fully by water at temperatures above , and by ice at temperatures below . Within the GH stability zone the GH saturation in the rock matrix porosity was inferred to be 30 or 60% range for two models, where 60% is the maximum saturation found in part of the Mallik well section (Scott Dallimore, personal communication, 2011). The formation water salt concentration, 9 g/L, was considered constant with depth and the pressure-temperature phase curves were adjusted to this value. The salinity of 9 g/L 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 (3): This formula was taken from paper by Galushkin (1997), [17]. If there is a fresh water within Mallik sediments, the liquidus temperature would be by 0.58°C higher. If there is a sea water with 40 g/L, the liquidus would be by 1.98°C lower. For temperature gradient 20 K/km it would mean a shift of the permafrost base by 30 m downward and by 100 m upward, respectively.
Numerical code performance was tested by comparing results against the analytical solidification problem solution [19], where a molten half-space at liquidus temperature, 1300°C, is in contact with a solid half-space at zero temperature and which releases latent heat of 477 kJkg−1 in the temperature range 1100–1300°C. Differences between the numerical and analytical temperature profiles were found to be within ~20°C. Assuming that the magnitude of the numerical/analytical difference is proportional to the temperature range the error expected for the IBP and GH numerical simulations should be about 100 times smaller (i.e., two tenths of a °C). A similar error range was estimated by halving the time and/or depth steps.
Our model uses deep heat flow, thermal conductivity, present IBP, type I (CH4) GH thicknesses, and a surface melting temperature (−0.576°C) that also considers formation water salinities (9 g/L). It employs latent heat effects throughout the IBP and GH layers, which improves upon previous models [14, 15]. The models are constrained by deep heat flow [20, 21], current IBP thickness and observed type I GH thicknesses [22, 23], as described above. The models consider the pressure-depth dependence of ice assuming hydrostatic and GH thawing points over the entire expected extent of the IBP and GH layers [24]. Previously published models have considered only a thin layer using a constant dissociation temperature [15].
In both the MD-BS and QEIG regions, it is uncommon to find GHs within the IBP, however, there are no methods that define GH within the ice-bonded layer. Logging signatures of GH and ice are similar and shallow section is not routinely logged, in large part due to drilling and hole stability problems in the IBP zone. This means we have no effective way to infer GH concentrations within the ice-bonded permafrost [25, 26]. Commonly the well surface casing is cemented in the stable indurated rock strata that occur below IBP base, but above the first GH occurrences. Rarely, GH occur in the IBP [25, 26], but the general lack of GHs in the IBP layer is consistent with models where most GHs form due to the in situ transformation of conventionally entrapped natural gas accumulations rather than a dynamic GH or IBP trap [1, 13]. This would appear consistent with other inferences regarding the origin of GH accumulations [1, 10]. In simulating the temporal subsurface changes in models where IBP and GH occurrences are separate and distinct, a consumption/release of the latent heat during decay/formation of IBP or GH was considered separately. In fact, GH formation/degradation occurs where the P-T field works in tandem with the complex petroleum system geological factors as illustrated by the work of Chen and Osadetz [10], and Smith and Judge [2]. Both show that rich GH accumulations are rare and the rich accumulations are commonly associated with conventional petroleum accumulations [10] as part of the thermogenic petroleum systems.
The model division into IBP and GH stability zones was prescribed explicitly. Because the IBP generally lacks GHs and the GH occur generally below the IBP only the water ice latent heat and P-T phase curve is considered in the IBP zone and that only the GH latent heat and P-T phase curve were considered in the GH stability zone.
The above conditions may not be appropriate in models where GHs are intercalated within the IBP layer. In such cases [25] at each depth point, it is prescribed which fraction of the pore space can be occupied by water ice or GH. Then the code checks independently the P-T conditions for IBP and for GH at each depth point, in each time step, and appropriately adjusts the volumetric heat capacity. 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, whereas GH forms first below 250 m.
3. General Paleoclimatic History
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. The evidence for past GST comes mainly from isotopic data and their analysis, especially O studies. Our models of surface forcing of temperature are generally based on known models of surface change [27, 28], (Figure 2). The recent and the Quaternary surface temperature forcing uses a detailed Holocene glacial-interglacial history compiled from other sources [15]. While we project the model results into the future we consider the climate effects based on a doubling of atmospheric CO2, that results in a local mean surface temperature increase of 2°C/100 yrs [29, 30].
The Paleocene-Eocene Thermal Maximum (PETM) could have been caused by abrupt releases of methane hydrates under the bottom of the ocean [31, 32].
During that time, the Arctic Ocean may have reached levels more typically associated with modern temperate (i.e., midlatitude) oceans. Following the Paleocene to early Eocene peak warming the climate cooled towards the Pleistocene glacial environment. However, 3 to 6 Myr ago, global average temperatures were still higher than today with surface temperature at poles higher than currently. During this interval the Northern Hemisphere likely lacked continental glaciers [27, 28]. Climate changed dramatically during the following 3 Myr before present and significant amplification of the climate response to orbital forcing (i.e., Milankovitch cycles) began. This caused changes in surface temperature forcing resulting in drastic oscillations between ice ages and warmer periods during the past 1 Myr. These resulted in cycles of glacials and interglacials marked in the Northern Hemisphere by the growth and retreat of continental ice sheets at frequencies initially 41 kyr and later 100 ky within a progressively deepening ice age. The gradual intensification of this ice age during the last 3 Myr has been associated with declining concentrations of the CO2 which partly influences temperatures changes. El Nino was continual rather than intermittent until 3 million years ago. The temperature change model that is applied to our study areas is shown in Figure 2.
We consider the variation between 0°C to −6°C (i.e., “warm” versus “cold” models), as a probable temperature range for the time 3 Myr ago, prior to cooling. In the “cold” model the amplitude of variations increased to 4°C by the onset of the “small” ice ages of the 41 ka period at time 2.5 Myr ago, and their mean temperature was decreasing stepwise from −8°C to −9°C and then to −9.5°C at time of the transition to the “big” ice ages of the 100 ky cycles 900 ka ago. The temperature variations during the “big” glacial cycles were 13°C. For the era prior to 2.5 Myr ago 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. In the “warm” model the simulation starts at time 3 Myr ago with steady-state T-z profile corresponding to GST of 0°C. GST is decreasing linearly from 0°C at 3 Myr ago to −10°C at the beginning of the large 100 ka glacial cycles 0.9 Myr ago. The GST history during the 100 ka glacial cycles is the same in both models. Using these constraints we model the permafrost and GH thickness in the changing surface temperature environment for the MD-BS and QEIG.
It is essential to understand the origin, growth, and persistence of subpermafrost GH accumulations, as a function of temperature history, to characterize a geological natural gas resource model. Our analysis will use the characteristics of IBP and GH origin, growth, and persistence, as a function of temperature history, to 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 providing a tool that will assist the appraisal of risks related to surface imposed changes on the GH layer.
4. Modeling of Permafrost and Gas Hydrate Stability Histories in the Mackenzie Delta-Beaufort Sea Region
Depth to GH base of type I stability in the Mackenzie Delta is well described [31] and it serves as an important constraint for the acceptance of model results. Recent analysis of the Holocene-Pleistocene temperature history using a 1D model of IBP and GH layers in an onshore Mackenzie Delta setting during the last 600 k years, Majorowicz et al. [13] showed that due to the buffering effect of the IBP and the retardation of GH dissociation due to latent heat effects GH can persist through interglacial intervals despite large variations in surface temperatures.
To better understand 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 and future climate change. 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 (Figure 1) under continental and continental relic IBP (accumulated during the Pleistocene marine low-stand prior to the Holocene sea level rise). More detailed descriptions of the MD-BS geothermal environment GHs and the Mallik well experiments are found elsewhere [23, 33–42].
4.1. P-T Time Envelope and Formation of Permafrost and Gas Hydrate: Preliminary Considerations
Considering the IBP and GH formation and due to the large uncertainty of surface temperature models for the remote past we used the temperature history range (“cold” versus “warm” model) shown in Figure 3 (i.e., Mallik case). Pleistocene surface history of glacials and interglacials is after Muller and MacDonald [27]. The more recent surface climate history for the end of the Wisconsinan and Holocene is after Taylor et al., [15]. The range of temperature critical for GH formation is shown in preliminary test for the steady-state in Figure 4.
Majorowicz et al. [13] showed 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, which is different than the onset time for the IBP. The equilibrium IBP thickness is some 250 m for a GST of −5°C. 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 GH formation.
4.2. Model for the 3 Myr History of Surface Temperature Change: Results
We have simulated the downward propagation of the surface warming and cooling attending the cyclical glacial and interglacial models for the eastern Richards Island/Mallik location based on the above-mentioned 2 models as shown in Figure 3.
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 ice are 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 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 temperatures and observed 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 surface-subsurface coupling. GH formation can start only when the long-term GSTs drops below −5°C.
4.2.1. Case 1: Gas Hydrate Formation Controlled by the Geological Gas Entrapment Conditions
In Case 1 we allow GH occurrence only in the porous and permeable geological reservoir confined under a preexisting impermeable seal (i.e., GH formation was considered below 900 m only). At Mallik, the interbedded character of the GH bearing strata also indicates a stratigraphic and lithologic control on GH occurrence. GH layers occur in coarse-grained sandstone beds separated by thin nonhydrate-bearing, fine-grained siltstone and claystone beds. The gas and GH at Mallik appears to be entrapped in association with the anticline structure. The fault bounding the structure is the likely conduit for migration of the gas into the GH stability zone. This is the probable situation for rich MD-BS GH in the underlying Beaufort-Mackenzie Basin (BMB). There GHs occur in close spatial association with underlying conventional petroleum accumulations, such as at the Mallik boreholes.
Two general types of gas were observed in the Mallik wells. Microbial gas is characterized by dry gas compositions (i.e., C1/(C2 + C3) > 1000) and methane carbon isotopic ratios between −70 to −93‰ [3]. Thermogenic gas is wetter and has carbon isotopic ratios for methane of around −35 to −45‰. The carbon isotopic ratios for thermogenic ethane and propane are −31‰ and −26‰, respectively. Methane isotopic compositions of 12 GH samples averaged −42.7‰ and clearly indicate a thermogenic source according to Lorenson (personal comm., 2005) and [3]. Thermogenic gas likely migrated upward along listric-normal growth faults from larger depths where gas generation is possible due to availability of the source rock and higher temperatures. The deep upward migrating gas could have been trapped in the anticline and/or tilted fault blocks and converted, much later, into GH when regional cooling took place. This hypothesis is consistent with the coincidence of GH occurrences with known underlying conventional petroleum accumulations [1, 10].
We relate our model to the above findings. The division into IBP and GH stability zones was prescribed in the model explicitly such that we have vertically separated the occurrence of IBP from GHs. A consumption/release of latent heat during the decay/formation of either IBP or GH is considered separately. As mentioned above this means that only the latent heat and P-T phase curve of water ice were taken into account in IBP zone, and only GH latent heat and P-T phase curve were used in the GH stability zone. GH formation was considered below 900 m only and the maximum likely GH saturation of 60% has been used.
In the resulting models GH did not form during the coldest 41 ka cycles (from −11.5°C to −7.5°C), because the steady-state geotherm corresponding to the mean surface temperature of the cycle, −9.5°C, enters the GH stability region just above 900 m (see Figure 4) 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. The model simulates the current base of IBP and GH at 541 m, and 1060 m, respectively. These are only slightly above the observed positions of IBP and GH in the Mallik wells at 600 m and 1107 m, respectively. The surface temperatures during the glacial, −15°C, and the interglacial were the same as that used in previous including model 3 of Majorowicz et al. [13] as modified from Taylor et al. [15].
The 3 Myr “warm” history considered as an alternative to the more detailed “cold” model of the ground surface temperature used in simulations at the Mallik wells gives very similar results to that of the “cold” model for the time interval covering the second half of Pleistocene, which is characterized by 100 ka long glacial cycles (Figure 5). The simulation starts at time 3 Myr ago with a steady-state T-z profile corresponding to GST of 0°C. GST is decreasing linearly from 0°C at 3 Myr ago to −10°C at the beginning of the large 100 ka glacial cycles 0.9 Myr ago. The GST history during the more recent 100 ka glacial cycles is the same in both models. Despite the warmer temperatures of the 3 Myr history, which is up to 6°C warmer for most of the period between 3 Myr and 0.9 Myr ago, in comparison to the previous 14 Myr model, the simulated thicknesses of the IBP and the GH layer differ by at most a few tens of meters between the two models. It means that the present temperature-depth distribution is not very sensitive to the remote GST history and that the results obtained by the 14 Myr history for the end of Pleistocene and Holocene are effectively similar for the period after 3 Myr ago. Despite our use of quite general surface temperature histories, the agreement between the two model predictions and the observed bases of GH and IBP is reasonably good.
4.2.2. Case 2: Considering Simultaneous Occurrence of Permafrost and Gas Hydrate
In Case 2 we permit the simultaneous occurrence of IBP and GH. This has some geological grounds as the geological seals are discontinuous laterally [3] and at least one such occurrence has been observed [25]. The model considered 30% porosity, with all the pore water frozen. The thermal conductivity of the frozen and thawed rock was 3.4 W/(m.K) and 2.1 W/(m.K), respectively, and the latent heat of water freezing is 0.334 MJ/kg. The numerical simulation of the subsurface temperature response to changes of the surface temperature forcing was calculated over the last 3 Myr (Figure 6). The GST history is the same as that for the Case 1 model (Figure 3).
In this case of simultaneous IBP and GH occurrence, the 3 Myr model of linear GST decrease between 3 Myr and 0.9 Myr illustrates well the interplay of the downward propagating surface cooling and the latent heat release during permafrost and gas hydrate formation. The first gas hydrate forms at the depth of 290 m, when the GST drops to −5.6°C. The latent heat of hydrate formation starts to be released, which decelerates substantially the downward propagation of the permafrost base that is, at that moment, at 250 m.
The formation of permafrost accelerates again only when the upward migrating upper boundary of hydrate meets the sediment occupied by permafrost. At that moment, the hydrate formation and therefore also its latent heat release stop there. Subsequently, the latent heat is released only at the downward-migrating bases of the IBP and GH layers.
It is questionable, what could have happened with the methane contained in the freezing rock. If all pore water had frozen, the methane 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 thermokarst areas exist under water bodies of the present Mackenzie River Delta.
Figure 6 shows that depth variations of the base of both the IBP and GHs are 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 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 cause its larger depth variations compared to model—Case 1.
In order to illustrate how the simulations are sensitive to hydrate saturation, two alternative levels, 30% and 60%, were considered. The results for simultaneous IBP and GH occurrence in the Mallik geothermal model and the “warm” 3 Myr ground surface temperature history are shown in Figure 7. As expected, before hydrate formation begins, the downward propagation rate of the permafrost base does not depend on the saturation. However, since the first hydrate formation at 290 m, both the IBP and GH bases migrate downward faster when the saturation is lower due to the smaller latent heat released during less saturated hydrate formation. This deficit cannot be compensated by the latent heat of ice released in the layer of simultaneous IBP and GH occurrence from the volume occupied in the higher saturation model by hydrate, because the value of latent heat of ice is lower than that of gas hydrate. During the 100 ka glacial cycles, the IBP and GH bases are by about 30 m and 40–50 m, respectively, deeper in the model with lower saturation. It is caused mostly by higher thermal conductivity within the layer of simultaneous IBP and GH occurrence, where the more conductive ice substituted for less conductive gas hydrate. The higher amplitude of the hydrate base variations in the less saturated hydrate model is caused by a smaller damping effect of the latent heat.
The Case 2 model provides no mechanism for methane to escape to the surface from below the older, but overlying IBP. This model predicts a simultaneous occurrence of IBP and GH below 250 m, which, however, is not generally observed [3]. The simultaneous occurrence also implies lower conductivity in this zone, which means a higher temperature gradient, which shifts the IBP and GH bases upward, above their currently observed position. The indicated upper GH boundary above 250 m shown in Figure 7 is a model prediction rather than an observed GH boundary. It was calculated according to the P-T phase curve. In reality, no GH was considered there because all the pore space above 250 m is occupied by water ice and GH cannot displace a volume already occupied completely by ice.
5. Modeling of Permafrost: Gas Hydrate Stability History in the Queen Elisabeth Islands, Canadian Arctic
GHs are inferred to occur in only 24% of the QEIG wells. In many cases the inferred GH occur at deep depths [40] that are much deeper than the inferred depth to the base of the Type I GH stability zone determined using similar geothermal data. Permafrost thickness in the QEIG are well described [21]. There permafrost base varies with depth. In the terrestrial part it is from tens of meters to 726 m but it is only few meters under the sea offshore, based on precise temperature logs in the area [21, 34–39].
Gas hydrate thickness onshore is to a large degree controlled by the extent of the permafrost and depth to the base IBP (circa −1°C). GH occurrence in the QEIG is relatively well known although the number of wells is small for such a large region [11, 42, 44]. According to Brent [44], his study of seismic data on southern Ellef Ringnes Island indicated that the permafrost base of 470 m is directly underlain and in contact with a deep GH layer that extends down to 930 m. In some cases GH detection from well logs exceeded the projected depth interval for the methane hydrate stability [11, 42]. However, well-log-based GH detection can be uncertain. The average depth to the detected GH base is 1,020 m and it is 368 m for the top of GHs.
Heat flow in Sverdrup Basin underlying the QEIG is typical of a continental margin sedimentary basin setting and it varies mainly between 50–80 mW/m2 [45, 46]. In such an environment of relatively high heat flow the IBP and GH stability is controlled by paleoclimate history of very low surface temperature and pressure forcing, partially due to the pressure attributed to the glacial ice sheet load.
In the QEIG area, the deepest base of subsurface terrestrial permafrost is inferred at 726 m [21]. GHs occur as deep as 1.5 km mainly a result of variable low surface temperature and past ice load history. Present occurrence of GH and overlaying permafrost is not necessary in equilibrium with present P-T conditions for GH stability as derived from present knowledge of temperature depth (T-z) and formation pressure (P). The best example of this is the dependence of permafrost thickness, which is highly variable in throughout the Arctic, and which depends on the history of marine regression and shoreline emergence during the last 9 Myr [47, 48], both of which are likely due to terrestrial rebound after ice unloading.
The QEIG began to rebound following the Wisconsinan glaciation [5, 49, 50]. Change in temperature from offshore submarine to terrestrial settings was up to 20°C and this caused a gradation of permafrost onshore as thick ice permafrost is not found offshore [39]. According to Taylor, [14] thick ice-bearing permafrost is not observed currently beneath the deeper channels of the central Queen Elizabeth Islands. Analysis of a precision temperature log obtained from an offshore well near Ellef Ringnes Island [40] indicates that the thermal regime beneath the seabed is in equilibrium with today’s marine environment. If thick permafrost similar to that observed on land today had existed in the Pleistocene in areas that are presently offshore, then such permafrost must have started melting no later than 25,000 years ago [14]. This suggests that the inter-island channels must have been water-filled at least by that date [14].
Currently, there are only 2-3 small ice caps from the Wisconsin glaciations in the QEIG. There was previously a controversy, now resolved, regarding the thickness extent in space and duration in time of any glaciations of the QEIG. In spite of the substantial postglacial uplift [4, 5, 50–52], which suggested a thick and continuous ice sheet [5, 53, 54] some used geomorphic evidence to suggest that no such pervasive ice sheet existed [50, 52, 55]. However, the proponents of the discontinuous Franklinian Ice Complex have recently acknowledged [56] presence and extent of the Innuitian Ice Sheet Model as correct.
Recent model results [57] give several important results relevant to glacial history in the Arctic: (1) an important characteristic of the 100-kyr cycle is its asymmetric structure, with long (90 kyr) fluctuating ice-growth phases followed by rapid (10 kyr) terminations; (2) 60–80% of the Laurentide Ice Sheet was cold-based (frozen to the bed) at the LGM, and therefore unable to undergo large-scale basal flow. The fraction of warm-based ice increases significantly through the ensuing deglaciation, with only 10–20% of the Laurentide Ice Sheet being cold-based by 8 kyr BP.
The study of Ice Sheet history shows that once the ice sheet formed following peak interglacial periods, at least its core remained intact through the entire glacial cycle, with ice thickness of ~2 km inferred for the Foxe Basin region. The volume of ice required to maintain this core of the Laurentide Ice Sheet represents at least 10 m of sea level. The Innuition ice sheet was cold-based (comparable to or colder than Greenland Ice Sheet which is today less than −13°C according to Tarasov and Peltier [56]).
Basal temperatures in the glacial ice models are controlled by air temperatures and complex ice sheet dynamics (e.g., accumulation phase, internal sliding, etc.). Recently, Rolandone et al. [58] report data that suggests that the Laurentide Ice Sheet was not frozen to its bed along its southern margin, implying basal sliding as suggested first by Payne [59]. In such a case the glacial ice is in a completely steady state. The ice sheets do not insulate the ground from severe air temperatures, because ice has a thermal conductivity very similar to the average thermal conductivity of earth materials. This effect is often referred to as “insulating” and we will use it here. It needs to be understood though, that it refers to steady state of ice load and that more realistic is a buffering by the ice in a dynamic process of its history. The basal ice sheet temperature is the upper boundary condition for GH model temperatures.
We will use present observation of permafrost and GH base depths to constrain and the surface forcing model for the glacial-interglacial history. These models, 21 constrained by present observations, will provide insight into the past environments that led to the initiation, growth and persistence of the linked occurrences of both permafrost and GH accumulations in the subsurface of the QEIG (Canadian Arctic Archipelago) models of the surface air temperature and ice cover dynamics (base temperature and pressure).
5.1. Recent Postglacial Rebound: Gas Hydrate-Permafrost Model
According to the recent model of Taylor and Wang [60, and references therein], the cooling of the terrestrial surface of the QEIG occurred soon after coastal emergence, a result of isostatic rebound within 200 years [60]. Different portions of the islands emerged at different times before the present mainly during the interval of hundreds of years ago to several thousands of years ago. This process is well described by the relationship between elevation and time of emergence [14, 60]. This process of temperature change from subsea, −1°C, to harsh terrestrial arctic temperature conditions, −20°C, controls, to a large extent, the subsequent IBP and GH formation. It is known from temperature-depth logs in QEIG wells (e.g., Sabine Peninsula, Melville Island, and Ellef Ringnes Island) that wells near the coasts are influenced by the recent history of emergence and some even show negative thermal gradient in the shallow section, while wells further inshore appear to have reached thermal equilibrium. Some [14] consider thick permafrost in the island interiors to be in equilibrium. In that sense, recent history would be the only factor controlling thickness of IBP and GH. We check that hypothesis with our modeling.
The QEIG model of a pure postglacial rebound considers only the effect of the islands emersion from the marine environment with a steady-state surface temperature of −1°C to the harsh terrestrial conditions of −20°C. It is assumed that similarly to [60], the cooling occurred within 200 years of emergence.
The geothermal model is:(1)rock porosity 30%;(2)gas-hydrate saturation 60% of the pore space (i.e., 18% of the total volume);(3)the remaining 40% of the pore space (not hydrate) is filled with either water or ice (i.e., 12% of the total volume);(4)conductivity of pure permafrost 3.4 W/(mK);(5)conductivity of the mixture of permafrost and GH 2.7 W/(mK);(6)conductivity of the GH 2.1 W/(mK);(7)conductivity of melted rock 2.1 W/(mK).
It is assumed that simultaneous occurrence of permafrost and GH is possible. The models considered the pressure-depth dependence of ice and GH thawing points over the entire expected extent of the GH and overlying permafrost. The first model assumes heat flow of 60 mW/m2 [45]. This model of surface forcing and the resulting variations in the permafrost and GH below is shown in Figure 8(a).
(a)
(b)
In Figure 8(a), the yellow part of the plot marks the depth range down to 255 m, where permafrost forms before the P-T conditions for GHs are attained. If we assume that all interstitial water is frozen in permafrost, then there is no space for the hydrate formation within this interval. The P-T phase curve of hydrate is first touched by the sinking temperature curve in the depth interval 74–91 m at 1240 years, but as mentioned above, the permafrost already occupies the available pore space. The depth interval, where the actual temperature is lower than the hydrate phase curve spreads from this time both upward and downward, but no hydrate can form until 2700 years, when the model lower boundary of GH stability crosses the permafrost base at the depth of 255 m. Subsequently, the GH forms prior to permafrost below 255 m.
Due to the assumed 60% GH saturation of the pore space the remaining pore space interstitial water, 40%, can freeze and a mixture of hydrate and permafrost can form as the cooling proceeds.
We also constructed models for heat flows values of 70 and 80 mW/m2, which could be possible values in the QEIG [45] for the first 10 k years of the model runs. Comparison of the results of the emergence model shows that the results for three versions of the heat flow of 60, 70, 80 mW/m2 differ only slightly due to differences in temperature gradient. The steady-state initial thickness of permafrost decreases from 23.1 m to 17.5 m with increasing heat flow. The time when the temperature curve touches the P-T curve of the hydrate increases with the heat flow (60, 70, 80 mW/m2) 1240 years, 1280 years, 1340 years, respectively. The time when the temperature curve drops to the crossing point of the gas/hydrate and water/ice phase curves increases with increasing heat flow from 2700 years via 2940 years to 3200 years. These are the times when the GH starts to form below the depth of 255 m prior to the water freezing (if the free gas is available, of course). The depth of the permafrost and the GH bases 5000 years after the land immersion decreases with increasing heat flow from 303 m and 335 m via 294 m and 319 m to 286 m and 304, respectively. Figure 8(b) shows T-z profiles versus stability curves for the time interval from 0 to up to 10,000 years. Time 0 equals a steady state at sea temperature −1°C; the time 1240 years equals the time when the temperature curve touched for the first time the phase curve of the hydrate; time 2700 years is the time when the temperature curve passes through the crossing of the ice/water and hydrate/water P-T curves.
According to these simulations, GH should not occur at sites which emerged less than 2700 years ago where the permafrost thickness should be up to 250 m. In areas, which emerged earlier than 2700 years ago, GH could occur just below the permafrost at depths greater than 250 m.
In some cases our models indicate current permafrost to a depth of 0.7 km and the hydrate layer is over 1 km [42] or close to 0.9 km deep, which is like that illustrated on Ellef Ringnes Island [44]. These depths are higher than ones possible for the models as shown above. It is possible, that marine conditions lasted only a short time and that the deep premarine subsurface temperatures could persist, at least partially. The sea level rose by some 100–120 m some 10 ka ago shortly after the end of recent glacial period. If the inundated area was above the glacial sea level during the last ice age, then the terrestrial conditions prevailed and controlled the permafrost hydrate thickness there for most of the glacial cycle. The inundation may have lasted only several thousand years. Alternatively, there may be lateral thermal effects between the exposed terrestrial setting and the inundated marine setting that are not well represented by our 1D models of coastline emergence.
5.2. Simulation of the Previous 3 Million Years
5.2.1. No Ice Insulation Case
In the first simulations of the last 3 Myr for the QEIG two surface temperature scenarios are considered, only without considering the phase curve changes due to the ice sheet pressure. The more complicate models that consider ice sheet pressure effects are discussed below. First scenario is simply a surface model constrained by the current −20°C in terrestrial conditions in the QEIG. Model 1 is shown in Figure 9.
It is apparent that assuming no ice sheet insulation, the typical heat flow results in a very thick GH stability zone base throughout the last 6 mln years and that the GH layer reaches 1.5 km in the most recent 100 ka year glacial maximum. While such base GH stability zone depths have been inferred for the Present [2, 40]; these depths are not typical and most GH in the QEIG occurs above depths of 1 km.
5.2.2. Ice Sheet Cover Case
The second surface forcing model (2) assumes that during the last 0.9 Myr, at a time 20 ka after the beginning of each of the 100 ka glacials episodes, the area was covered by an ice sheet and that the temperature increased due to the insulation from −29°C to −15°C, which is a cold-based glacier scenario like the current Greenland ice sheet. In model 2 the the interglacial temperatures are the same as those in model 1. The model 2 results for the IBP and GH layers are shown in Figure 10, where it is compared to model 1.
It is apparent that the introduction of ice sheet insulation effects, in model 2, reduces the thickness of both the IBP and Type I GH layers significantly. The IBP and the GH stability zone bases are 0.7 km and 1.1 km, respectively, which is closer to the currently observed values then the Model 1 result predicted. On the other hand, the introduction of ice sheet insulation increases the variability of permafrost and GH thickness between interglacial and glacial intervals.
5.2.3. Ice Buildup and Thaw Case
The subsequent simulations consider the ice sheet basal temperature and thickness during the last 0.9 Myr, when the large 100 ka glacial cycles occurred. We are specifically interested to produce models that illustrate how the base of the GH stability zone responses to both ground temperature and pressure changes for each of the nine 100 ka glacial-interglacial cycles assuming progressive ice sheet formation and melting. Below we show model results that illustrate ice sheet pressure effects on GH stability zone depth and history.
Figure 11 shows that the downward migration of the GH zone base decelerates with increasing ice thickness because of the steeper position of the phase curve at higher pressures. The ice sheet thickness increase from zero to 2 km by a step of 500 m means the base shifts downward from 1344 m to 1604 m by successive increases of 97 m, 67 m, 52 m, and 44 m. The maximum effect of the glaciation on GH thickness increase is therefore less than 250 m.
The modeling process that considers glaciation and ice sheet loading is illustrated in Figure 12. The sudden pressure increase due to the ice load causes the hydrate equilibrium temperature to increase at a given depth. The old and new positions of the base of the GH stability zone are given by intersections of the T-z profile with the old and new phase curves. In this case the base of the GH stability zone drops by 97 m and amounts to 1344 m and 1441 m, respectively (Figure 12).
If we assume that there is enough gas in the water-filled pores of this 97 m thick sedimentary rock layer, the process of hydrate formation within it starts immediately. Hydrate crystallization is an exothermic phase change. For the assumed porosity of 30% and the hydrate saturation of 60%, the amount of crystallization heat, that would be released by crystallization of all the possible hydrate, is by an order of magnitude larger that the heat necessary to warm the surrounding rock from the actual temperature to the equilibrium phase temperature. It means that the temperature of the rock, water, gas, and hydrate mixture in this zone reaches the equilibrium temperature long before all the possible hydrate can be formed. The fraction of the hydrate formed at a given depth by the time when the equilibrium temperature is reached is proportional to the difference between the equilibrium temperature and the original temperature. The software used for the simulations is based on the concept of apparent specific heat, where it is assumed that the crystallization heat is released or consumed during hydrate formation/decay in a temperature range constrained by the solidus () and liquidus () temperatures. In the new simulations for the ice thaw and refreezing, the position of the hydrate phase curve is time dependent and it is prescribed in each time step. If the curve position changes between the two consecutive time steps, the actual temperature within the “disequilibrium” layer (in the Figure 11 between 1344 m and 1441 m) is increased to the (, ) temperature range. The different fraction of the formed hydrate in the individual depths of the zone is taken into account by setting the temperature between (, ) in such a distance below , which is proportion to the hydrate fraction at the given depth. The (, ) range used in the calculations was 0.2°C. An analogical process/treatment, but in a reverse direction, was considered and applied for the opposite case of ice sheet melting and dissipation.
5.2.4. “Realistic Ice” Model
Variations of ground surface temperature and the bases of IBP and GH is referred here as the “realistic ice” model including pressure effects. The model differs from the previous two surface temperature models, the “Mallik –14°C” and the “simple ice” without pressure effect models (Figures 7 and 8) only in the last 0.9 Myr in the period of the large, 100 ka long, glacial cycles. The model is the same for each of the 9 glacial cycles, as is shown and explained in detail in Figure 13.
The Figure 13 model uses the repetition of the cycle shown in Figure 14 (last 0.9 Myr). The first 75 ka of the 100 ka cycle is the glacial part followed by 25 ka long interglacial. The upper time scale shown in that figure is the time of one 100 ka cycle and the lower time scale is that of our era. It assumes that at the present (time 0) we are 13.5 ka after the end of the last ice age and remaining 11.5 ka of the interglacial is still ahead. The model assumes the interglacial temperature course similar to the Mallik case, but shifted by −14 K. The glacial temperature at the ice sheet base depends on the thickness of the ice and was considered as:(i)−29°C for the first 20 ka of the glacial, when no ice was assumed,(ii)−24°C between 20–25 ka and ice thickness 0.5 km,(iii)−21°C between 25–30 ka and ice 1 km,(iv)−18°C between 30–35 ka and ice 1.5 km, and(v)−15°C between 35–70 ka and ice 2 km.
At the end of the glacial, the ice thickness drops from 2 km to 1 km at the time of 70 ka and to zero at 75 ka. The pressure effect of the ice overburden on the hydrate equilibrium temperature was considered in the first approximation as a hydrostatic one, ignoring the density difference between water and ice. The depth (pressure) dependence of the equilibrium temperature T(z) at depth z was given by (4): Both, the depth and the ice thickness are given in meters. The pressure dependence of the water ice (permafrost) equilibrium temperature was not considered. Because of its very small (negative) vertical gradient, considering this dependence would influence results of the simulations negligibly.
In the model that considers ice pressure we assume that the pressure increase is hydrostatic and that the effect of a 1 km thick ice sheet corresponds to a 1 km thick water column. As can be seen, the base of the GH stability zone migration decelerates with the increasing ice thickness because of the steeper position of the phase curve at higher pressures. The ice sheet thickness increases from zero to 2 km by a step of 500 m meaning that the base moves downward from 1344 m to 1604 m by increases of 97 m, 67 m, 52 m, and 44 m. The maximum effect of the glaciation on hydrate thickness increase is therefore, less than 250 m, which is rather small. A variety of glacial ice models, including: no ice cover; “simple ice” cover without pressure effects and a “realistic ice” cover with pressure effects illustrated for comparison in Figure 15.
In order to bring the model GH base upward to more accurately predict the currently observed GH base of about 1 km, an alternative surface temperature history was considered. It assumed an 8°C higher temperatures during the ice free periods of the glacial cycle than those in the model “Mallik –14°C”, that is, in the periods 0–20 ka and 75–100 ka. Ground surface temperature below the ice sheet was considered −15°C, independently of the ice sheet thickness. Ground surface temperature used in the simulations and the calculated depth variations of the permafrost and GH bases during the last 100 ka are shown in Figure 16.
Currently the simulated depth of IBP and GHs occurs at 572 m and 1072 m, respectively. It is slightly more than the depths indicated by geophysical research. The ground surface temperatures can be hardly higher in the ice-free periods than those considered in the +8 K version. The simulated values could be brought closer to those observed on the seismic profile (0.9 km) by considering higher basal temperatures under the ice sheet, (−10°C). Another explanation might be that there is an overestimation of model rock thermal conductivity values, although these are only speculations.
6. 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 the last 14 Myr of GST history. In the MD-BS (i.e., Mallik) case our model shows that the onset of GH formation always begins after the IBP layer is initiated. IBP provides a potentially impermeable seal that could entrap any upwardly migrating gas. The IBP seal is inferred to have come into existence prior to when P-T conditions permit the first GH to have formed. Therefore, if it were common for GH to form due to the IBP providing a vertical migration barrier it would be common for GH to be found intercalated with water ice within the IBP, since the IBP grows subsequently downward through the GHs. While intercalated GHs have been observed it is rare [25] and it appears more likely that most GH accumulations result from the in situ transformation of conventionally entrapped natural gas below a lithologic or stratigraphic seal.
The results of the Case 1 MD-BS/Mallik models permit the formation of GH from natural gas previously entrapped under a deep geological seal and the Case 2: models permit simultaneous GH and IBP formation due to a persistent gas flux into the free pore space result in two entirely different times when GHs first begin to form. For Case 1, GHs could have begun to be formed at 0.9 km only about 0.9 Myr ago. For Case 2, the first GHs formed earlier
Where the IBP layer is thick, it is likely that sub-IBP some GHs persisted through the 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. This is also likely scenario for the offshore MD-BS case where hydrates can and have probably existed under relic IBP despite marine transgression [13]. Likewise, the impact of future climate change attending a doubling of atmospheric CO2 does not completely destabilize the GH layer. Majorowicz et al., [13] predicted that future GH layer thinning is very small, and within the range of previous natural cycle variations, in spite of the accelerated surface warming accompanying climate change. This model results appear consistent with the implications of recent observations of methane isotopic compositions from ice cores for the “clathrate gun” hypothesis of Sowers et al., [61]. However, the GHs can destabilize rapidly in response to environmental change in the marine non-IBP GHs [62, 63].
In the QEIG Sverdrup Basin case, our models show that thermal and pressure effects of the repeated formation and disappearance of the Innuitian Ice Sheet during glacial and interglacial periods is needed to successfully match the observed base of the GH stability zone at about 1020 m. Without the Innuitian Ice Sheet the exposed surface temperature would have been much lower resulting a hydrate stability zone thickness much thicker than the observed values of up 1.5 km. Detailed modeling of the recent history of coastal emergence due to isostatic uplift and the attending marine to terrestrial (−20°C) temperature changes results in a base of the GH stability zone depth that is less than 0.5 km. These are lower than the observed average depth of base of the GH stability zone of just over 1 km [42]. It is possible that marine conditions have lasted for only a very short time in the shallow interisland channels and the deep preinundation subsurface temperatures persisted, to some degree, while sea level rose by some 100–120 m some 10 ka ago, shortly after the end of the most recent glacial period. If the subsequently inundated area was above the glacial sea level during the last ice age, then those terrestrial conditions could have prevailed and controlled the permafrost and gas hydrate thickness for most of the glacial cycle. The Holocene coastal inundation lasted probably only several thousand years. Alternatively, there may be lateral thermal effects between the emergent islands and the inundated shore face that are not well represented in our 1D models.
One of the uncertainties of the model scenarios is the subice temperature. Temperatures below the ice are assumed based on the observation that current temperatures below the Greenland ice sheet may be as low as −15°C. The thickness of Greenland’s ice sheet is 2-3 km. We employed representatively lower temperatures assuming a thinner ice coverage over most of the QEIG.
Acknowledgments
This research was funded by the Geological Survey of Canada. The authors would like to make special thanks to Alan Taylor for his helpful suggestions and informed advice, which helped them to build the models and develop their interpretation. Scott Dallimore is to be thanked for his in-depth review of the Mallik case scenario, which led to significant changes of our previous models as to surface forcing scenarios and saturation range. Other errors or omissions are the authors own.