Corps de l’article

Introduction

A large part of the boreal zone of the western Canadian Arctic is underlain by discontinuous permafrost, much of which is at temperatures close to 0 °C (Burgess and Smith, 2000). It is an area subject to considerable climatic fluctuations (Kistler et al., 2001) and is prone to climatic warming (Zhang et al., 2000), both of which have major implications for landscape and hydrological processes and on the construction and maintenance of northern infrastructure including oil and gas pipelines. Thawing of ice-rich permafrost in response to climate warming or surface disturbance can lead to settlement of the ground surface due to a large loss of volume when the ice is melted. This often has significant impacts on the integrity of northern infrastructure (Smith et al., 2001; Couture et al., 2003). Characterization of present and future permafrost conditions is therefore critical to facilitate the appropriate design of northern infrastructure and also to evaluate the environmental impacts associated with northern development.

An existing buried oil pipeline traverses the discontinuous permafrost zone of the Mackenzie Valley and transports oil from Norman Wells, Northwest Territories to Zama, Alberta (Fig. 1). In addition, there is a proposal to construct pipelines to carry natural gas and natural gas liquids from the Mackenzie Delta south through the Mackenzie Valley to northern Alberta. Both the existing and proposed pipeline routes traverse the entire width of the boreal, discontinuous permafrost zone, an area that is sparsely populated. Since the mid 1980s, the Geological Survey of Canada (GSC) has maintained a permafrost thermal monitoring network in the Mackenzie Valley that includes 20 monitoring sites along the Norman Wells pipeline corridor. The data generated from the network provide information on present permafrost conditions and recent trends. The network of observed permafrost conditions however is sparse. A simple yet physically-based model is desired to simulate thawing of the active layer. This study proposes to use Stefan’s algorithm to calculate the thawing of ice-rich soils and to investigate the effects of summer temperature and near-surface moisture on the seasonal thawing of materials commonly found along the Mackenzie Valley pipeline corridor. Data from the GSC thermal monitoring network are used to validate the model performance.

Physical Setting Of The Mackenzie Valley

The present climate of the boreal region within the Mackenzie Valley is characterized by long winters with mean January air temperatures (based on Environment Canada’s 1971-2000 climate normals; available online at www.climate.weatheroffice.ec.gc.ca/climate_normals) ranging from ‑26.5 °C in the central Mackenzie valley (Norman Wells; 126° 48’ W, 65° 16’ N) to ‑21.6 °C in northern Alberta (High Level; 117° 09’ W, 58° 37’ N), and mean July temperatures reaching 17 °C in the central valley and 16.2 °C in northern Alberta. Temperature generally increases southward although local sites can be anomalously warmer or cooler due to topographic and vegetation influences. Normal annual total precipitation ranges from 290 mm in the central valley to 394 mm in northern Alberta, about half of which comes as snowfall that often stays on the ground from October until April. Snow depth exhibits considerable spatial variability and is controlled by local conditions. Snow cover is effective in insulating the ground against extreme winter coldness, as is shown by a comparison of the air and near-surface (0.05 m depth) temperatures at a site near Wrigley (Fig. 2).

The present physical landscape of the Mackenzie valley is primarily a result of the last continental glaciation which covered most of the region about thirty thousand years ago and most areas are underlain by unconsolidated glacial and post glacial deposits. A series of large glacial lakes were impounded by the Laurentide Ice Sheet as it retreated rapidly into the present Mackenzie valley about 13 ka BP (Duk-Rodkin and Lemmen, 2000). Thick (up to 30 m thick in some places) extensive deposits of glaciolacustrine and lacustrine silt and clay are associated with these lakes (Aylsworth et al., 2000). The post-glacial landscape comprises morainic and fluvial landforms of the northern Interior Plains. Most parts of the undulating valley are covered by boreal forest dominated by spruce, shrub undergrowth and a moss-lichen floor. Many areas have impeded drainage which favours the formation of wetlands, ponds and lakes. Peat development is extensive particularly in the southern portion of the valley (i.e., south of Fort Simpson) (Zoltai et al., 1988; Aylsworth and Kettles, 2000), covering the mineral soils with an organic cover that has thicknesses varying from centimetres to metres. Permafrost is found frequently under thick peat which prevents summer heat penetration and preserves ground frost. Poor drainage ensures frequent saturation of the peat, silt and clay, which when frozen, have a large ice content. The ice-rich lacustrine sediments and peat are highly sensitive to disturbance as thawing of the ground is commonly accompanied by subsidence (Burgess and Smith, 2003) or slope failure (Dyke, 2000).

Figure 1

Boreal region of western Canada bounded by the Arctic tree-line in the north and the temperate forest in the south. Also shown are permafrost zones (from Heginbottom et al., 1995), Mackenzie Valley pipeline routes (both the existing oil pipeline south of Norman Wells and the approximate route for the proposed gas pipeline north of Norman Wells) and sites mentioned in the text.

La région boréale du Canada occidental est délimitée par la limite des arbres au nord et la forêt tempérée au sud. Les zones de pergélisol (d’après Heginbottom et al., 1995), les tracés des pipelines de la vallée du Mackenzie (tant le tracé actuel du pipeline de pétrole au sud de Norman Wells que le tracé approximatif du pipeline de gaz proposé au nord de Norman Wells) et les sites mentionnés dans le texte sont également indiqués.

-> Voir la liste des figures

Figure 2.

Comparison of daily air and near-surface (0.05 m below ground) temperatures for a site (kp182) near Wrigley, NWT, for 2001.

Comparaison des températures journalières de l’air et de la surface du sol (0,05 m sous le sol) d’un site (kp182) près de Wrigley, Territoires-du-Nord-Ouest, en 2001.

-> Voir la liste des figures

Thaw Process and Stefan’s Equation

Heat flux in the soil occurs through convection and conduction. Heat convected with snowmelt and rain water is usually small in magnitude in an Arctic environment, though the heat released by refreezing of snowmelt water that infiltrates the frozen soil in the spring can rapidly raise the temperature of large parts of the active layer to around 0 °C (Woo, 1986). Conduction is the principal mechanism of ground heat transfer but the amount of heat flow is strongly dependent on the thermal conductivity of the soil material. Furthermore, for ice-rich soils, a large portion of this heat is consumed by the melting of ground ice (Woo and Xia, 1996). For these soils, the near-zero temperature can remain for a protracted time (days) at a particular depth where ice is abundant, producing the zero-curtain effect (in Figure 3, the temperature at 1 m depth hovered around ‑0.1 °C for eleven days after the thawing front arrived on May 25th). Normally, a small quantity of the ground heat is spent to raise the soil temperature but if the frozen substrate is much below 0 °C, a large amount of the heat may be conducted away from the active layer to warm the permafrost below.

Figure 3

Ground temperature measured at several depths at site kp182 during 2001 (top), and comparison of simulated thaw front with ±0.5 °C isotherms interpolated from these temperature measurements (bottom).

Température du sol mesurée à plusieurs profondeurs au site kp182 en 2001 (haut) et comparaison du front de dégel simulé avec des isothermes de ±0.5 °C interpolées à partir de ces mesures de température (bas).

-> Voir la liste des figures

Stefan’s equation, originally formulated for lake ice melt, has been adapted for ground freeze-thaw calculations (Jumikis, 1977):

where dzt /dt is the rate of thaw front descent (m/s), kt is the thermal conductivity of the thawed soil (J/m s K), Dt is the thawing degree days, λ is the volumetric latent heat of fusion of water (J/m3), and θz is the volumetric fraction of soil moisture content at depth z.

Thawing is initiated at the top of the soil column and heat flux to the freeze-thaw front is conveyed through conduction. Heat convection is ignored and all available heat is used for phase change, not for warming or cooling of the soil. Ground near-surface (0.05 m) temperature is the forcing that drives the heat to the freeze-thaw front through conduction, which depends on thermal properties of the soil.

Freeze-thaw calculated using Stefan’s algorithm has been tested on the soils at a number of sites that range from Arctic, through subarctic, to temperate locations (Woo et al., 2004). The following procedures are used to compute the position of the thaw front.

(1) The soil column is divided into n slabs (e.g., n = 10) of ∆Z thickness, each with its specified physical properties including bulk density (ρb), porosity (φ), mineral and organic contents (fm and fo), and minimum unfrozen moisture content (θ min) as a lower limit that does not freeze in the environment under consideration.

(2) The soil moisture content of each slab is estimated from field data (or for saturated soil, it is set equal to porosity) to enable the calculation of kf(thermal conductivity of the soil above the front, in J/m s K) (Farouki, 1981):

where

with k(j) and f(j) being the thermal conductivity and volumetric fraction of the j th type of substance making up the soil (including mineral and organic materials, ice, water and air).

(3) The thermal resistance of various slabs is calculated to represent the resistance to the transmission of the thawing temperature down the soil column.

where zf is the thickness of the soil above the front and R has the unit of m2 s K/J.

(4) When the front reaches slab m, the warmth needed to thaw this slab is:

where 0≤∆x<∆Z is the thickness of the frozen zone within slab m.

(5) The thawing degree-days that descends from the top layer is (T-To )∆t, with T being the mean temperature of the top soil layer for the period ∆t, and To being the freezing point of the soil moisture. This is expressed as the mean ‘degree-day’ condition, D = (T-To )∆t (in s K), which is compared with the warmth needed to thaw the slab. If D>N, slab m will thaw, and the front will proceed to slab m+1, while the remaining ‘degree-days’ available to thaw the lower slab is (D-N). Through repeating these steps, the thaw front continues to descend until the remaining ‘degree-days’ are exhausted.

(6) In the case of D<N, there is only partial thawing of the slab and the thaw depth within layer m is:

This value is added to ∑∆Z for the m‑1 slabs to get the depth of the new thaw front. In the following day, if the new temperature exceeds the freezing point of the soil moisture, further descent of the thaw front will be calculated by repeating steps 4 to 6. Calculation terminates on the day when T<To.

Thawing of frozen ground implies the loss of ice in the soil and this does not necessarily correspond to the soil temperature reaching 0 °C and above. Complications are due to depression of the freezing point for soil-water with impurities, and the frequent presence of unfrozen water in frozen soils (Burt and Williams, 1976). For algorithms such as Stefan’s equation, ground thaw is defined on the basis of total elimination of ground ice. There is therefore a discrepancy between the thaw front thus obtained and the 0 °C isotherm. Additionally, there are uncertainties due to measurement error and interpolation of isotherms so that for practical considerations, the position of the thaw front is considered to reside between the ‑0.5 and +0.5 °C limits.

Data

Two types of data were obtained: (1) ground temperature for simulation and verification of results, and (2) soil properties for estimating the required parameters for model runs. Ground temperatures to depths of up to 20 m have been measured since the mid 1980s by the GSC at several sites representing a variety of soil, vegetation, ground ice and thermal conditions. Temperature cables with YSI44033 thermistors have been placed in PVC casing installed in boreholes and filled with silicone oil. The accuracy of the temperature sensors is ±0.1 °C while the measurement system allows for a resolution of ±0.01 °C. Temperature measurements are made both manually or with dataloggers. Near-surface temperatures are measured every 4 hours with single channel loggers placed just below the ground surface. The near-surface temperature is recorded to an accuracy of ±0.5 °C with a resolution of ±0.3 °C.

Soil cores were obtained during drilling of the boreholes in 1984 and 1985 at a number of sites. These cores were used to determine bulk density and moisture contents. Porosity estimates were based on bulk density and assumed particle density for peat and mineral soil. In addition, time domain reflectometry (TDR) was used to determine unfrozen water contents and its variation with temperature. Further information on laboratory analysis of soil cores can be found in Patterson and Riseborough (1988) and Patterson et al. (1988 and 1991).

Model Testing

The model for ground thaw was tested on a permafrost site (kp182; 124° 28’ W, 64° 16’ N) between Tulita and Wrigley, Northwest Territories. This site provides daily temperature measurements at depths of 0.05 (near-surface), 0.75, 1.00, 1.125, 1.25, 1.50, 1.70, and 1.90 m. The near-surface data were used to drive the model and all temperature data were used to locate the depth of the thaw front through interpolation of the ±0.5 °C isotherms. The site is in a spruce woodland with shrub undergrowth (Fig. 4A). The soil is ice-rich silty-clay. Table I summarizes the measured and estimated soil parameters for various depths used in the simulation. Since measured near-surface temperature obtained below the vegetation mat is used as the forcing, the insulating effects of snow and of the thin living plant layer are already accounted for in the input data.

Ground thaw was computed for the thaw periods of 2001, 2002 and 2003. Only the result for 2001 is presented because the performance of the model is similar for the remaining years. Figure 3 (bottom) shows that the simulated thaw front lies within the bounds of the ±0.5 °C isotherms for most of the summer except during the initial thaw period. The simulated maximum thaw penetration of 1.8 m is in agreement with observation. Retardation of the advance of the thaw front after July reflects the presence of abundant ground ice (and the associated large latent heat flux), indicated by the zero-curtain featured in the temperatures measured at these depths. The discord in the initial descent of the thaw front is attributed to heat convection by snow meltwater entry into the frozen soil, a phenomenon that cannot be represented easily without coupling to the complex land surface and hydrological processes. Given that the primary purpose of the simulation is to obtain the maximum thaw depth using a limited amount of data, the model performance is considered to be satisfactory.

Sensitivity Simulations

To further explore the effects of climate and soil types on active layer depth, the Stefan’s model was applied to sites in the Mackenzie valley where ground surface temperatures are available for the summers of 2001 and 2002. These sites (Fig. 4) include the Norman Wells area (126° 48’ W, 65° 17’ N) at the northern end of the oil pipeline, the Fort Simpson area (120° 42’ W, 61° 11’ N) at the middle of the pipeline route, and the Zama area (119° 37’ W, 59° 32’ N) at the southern end (Fig. 1). We tested the sensitivity of active layer thaw in three typical soil conditions: thick peat, mineral soil, and thin (0.2 m) peat on mineral soil. The peat assumes a low bulk density of 60 kg/m3 and a porosity of 0.95. The mineral soil is typically silty fine sand with a bulk density of 1 800 kg/m3 and a porosity of 0.32. The unfrozen water content in all frozen materials is taken as 5%.

Figure 5 shows the results of the simulations. The observed near surface temperatures in 2001 (Table II) for the three sites were applied to the soil columns, first considering the entire columns to be saturated and then re-simulating under drier condition where the top 0.2 m has a moisture content of only 10%. Note that the most southerly site (near Zama) has a lower near surface temperature than the most northerly site near Norman Wells. This is partly due to the effects of the higher elevation associated with the Alberta Plateau and also to differences in exposure, daylight hours, precipitation and drainage conditions at each site. Under full saturation, a situation commonly found in areas of impeded drainage, one obvious effect of soil materials on thaw penetration is the large contrast between peat and mineral soils, the former thawing to a maximum depth of about 0.5 m but the latter thawing to over 1.5 m (Table III). Such a difference is due to the large amount of ice in the peat that has to be melted and the associated latent heat requirements. The presence of a thin peat layer on top of mineral soil also slows down the progress of the thaw front until the front passes below 0.2 m. Then, the descent quickens in the mineral substrate. Under the condition of a dry top layer, the effect of insulation by the porous peat layer is evident, due to the considerable presence of air in its pores (thermal conductivity of air is 0.025 J/m s K, compared with a conductivity of 0.570 J/m s K for water and 2.200 J/m s K for ice). Thus, compared with the situation of complete profile saturation, maximum thaw is only slightly reduced in the mineral soil (about 0.1 m), but is halved in the peat. In the case of peat-over-mineral soil, the thaw front does not pass much beyond the peat layer so that the mineral soil remains largely frozen throughout the season. A peat layer thickness of 0.5 m or more would hinder thaw penetration into the mineral soil layer below, as was noted by Burgess and Smith (2000) at a number of sites near Fort Simpson. This is illustrative of permafrost preservation under dry peat, a condition common in the discontinuous permafrost region. In fact, permafrost is largely limited to organic terrain at the southern fringe of the discontinuous permafrost zone and this permafrost is likely a relic of the Little Ice Age (Vitt et al., 1994; Halsey et al., 1995).

The summer of 2002 had a mean near-surface temperature very close to that of 2001, but the duration of the thaw season was about 20 days shorter (Table II). Both the temperature and the number of thaw-days affect ground heat input, which is represented by the thawing degree-day factor in Stefan’s equation. In response to the difference in heat input, the simulated maximum thaw depths for 2002 are shallower than their 2001 counterparts, though the general patterns of thaw front descent are similar. The reduction in maximum thaw depth is especially noticeable in the mineral soil. This shows that in areas with large climatic variability such as the Mackenzie valley (Kistler et al., 2001), much fluctuation can be expected of the maximum frost table depth in mineral soils.

Figure 4

A

B

C

D

Photographs of monitoring sites referred to in the text: (A) kp182 (photograph: D. Riseborough), (B) Zama area (photograph: S. Smith), (C) Norman Wells area (photograph: D. Riseborough), (D) Fort Simpson area (photograph: S. Smith).

Photographies des sites de contrôle mentionnés dans le texte : (A) kp182 (photographie : D. Riseborough), (B) région de Zama (photographie : S. Smith), (C) région de Norman Wells (photographie : D. Riseborough), (D) région de Fort Simpson (photographie : S. Smith).

-> Voir la liste des figures

In general, simulations using the Stefan’s model demonstrate the sensitivity of active layer thaw to (1) soil materials due to their differential thermal properties, (2) moisture content which largely controls the latent heat requirement for phase change, and (3) inter-annual variations in near-surface temperature and duration of thaw season, which are governed by the climate.

Discussion and Conclusions

The vast boreal region has sparse thaw depth information, yet its discontinuous permafrost is sensitive to changes in the climate and the environment. The GSC operates an active layer monitoring network in the Mackenzie valley (Nixon, 2000) but only maximum annual thaw penetration is determined. Observations of the progression of thaw throughout the summer are recorded at only a few sites. The Stefan’s algorithm offers a robust method (a model which requires few assumptions) to calculate thaw penetration in frozen soils, requiring only near-surface temperature as the forcing variable, a description of the thermal properties and the moisture content of the soil columns. It has a strong physical basis, grounded on the assumptions that conduction is the primary mechanism of heat flow and all the available heat is used for melting of ground ice. Conduction is the only heat transfer mechanism included in the model and convective heat flow which may be associated with infiltration of snow meltwater into the soil is not accounted for. Model performance however is considered to be satisfactory as the simulated and observed thaw penetrations are in agreement.

We have used Stefan’s algorithm to assess the sensitivity of ground thaw to various factors including soil type (mineral vs. organic), moisture content and heat input (ground surface temperature), all of which are expected to vary considerably in the boreal region. In addition to natural variations in soil properties and hydrological setting, human activities, such as vegetation clearing, removal of the organic layer and other surface disturbance, alter the soil structure and drainage. For this area, large short-term variations and long-term changes are expected in the regional climate, as well as local modification of the microclimate due to wildfires, deforestation or changes in snow distribution. The Mackenzie Valley, like many other parts of the western boreal region, has strong potential for environmental changes induced by either natural forces or economic development. All these circumstances influence the factors that control active layer thaw, which can be easily assessed using the Stefan’s model. The model can facilitate the assessment of the impacts of northern development, climate change and other disturbance on active layer conditions and provide information required to mitigate and adapt to these changes.

Table I

Parameters used in the calculation of ground thaw at site kp182 near Wrigley, NWT

Parameters used in the calculation of ground thaw at site kp182 near Wrigley, NWT

A: entire soil column saturated.

B: top 0.2 m with 10% moisture content, bottom layer saturated.

-> Voir la liste des tableaux

Table II

Mean near-surface temperature, duration, and thawing degree days of the thaw period in 2001 and 2002, for three sites

Mean near-surface temperature, duration, and thawing degree days of the thaw period in 2001 and 2002, for three sites

A: entire soil column saturated.

B: top 0.2 m with 10% moisture content, bottom layer saturated.

-> Voir la liste des tableaux

Table III

Simulated maximum depth of thaw (m) for three soil profiles, using near-surface temperatures from three sites obtained in 2001 and 2002

Simulated maximum depth of thaw (m) for three soil profiles, using near-surface temperatures from three sites obtained in 2001 and 2002

A: entire soil column saturated.

B: top 0.2 m with 10% moisture content, bottom layer saturated.

-> Voir la liste des tableaux

Figure 5.

Descent of thaw front in three typical soil profiles, simulated using near-surface (0.05 m) daily mean temperatures measured near Norman Wells, Fort Simpson and Zama, for 2001 and 2002. Soil conditions are assumed to be entirely saturated (full lines), or with a dry 0.2 m top layer (dashed lines).

Abaissement du front de dégel de trois profils de sol typiques, simulée à partir des températures moyennes journalières proches de la surface du sol (0,05 m) près de Norman Wells, Fort Simpson et Zama, pour 2001 et 2002. Les sols sont présumés être entièrement saturés (lignes pleines) ou présenter une couche sèche de 0,2 m (lignes pointillées).

-> Voir la liste des figures