Introduction
Recently initiated coal mining in Svea Nord, Svalbard, has encountered problems due to a subglacial talik (an unfrozen region within an otherwise frozen region (permafrost)). Mining in permafrost regions, where it may encounter taliks, is a known challenge (Reference Tolstikhin and TolstikhinTolstikhin and Tolstikhin, 1976), since, as in this case, a permeable talik enables the leakage of water into the mine. The leakage from the usually continuous hydraulic system of the subpermafrost aquifer can be difficult, or even impossible, to stop. Currently the mining company, Store Norske Spitsbergen Kullkompani (SNSK), pumps the intruding water out of the mine, a procedure that consumes much energy. Hence, the issue arises as to whether the water can be drained gravitationally along the glacier bed.
Høganesbreen (77°57’ N, 16°42’ E) is a small High Arctic valley glacier, located on the northern side of Van Meijenfjorden, south Spitsbergen, Svalbard (Fig. 1). It is 13.4 km 2 in area, 1–1.5 km wide and 6.6 km long, possessing four small tributary glaciers (Reference Hagen, Liestøl, Roland and JørgensenHagen and others, 1993). Høganesbreen is part of the drainage from Gruvefonna (900ma.s.l.) and Sjaktbreen (700ma.s.l.). Its tongue terminates at approximately 140 m a.s.l. in a front moraine complex. Like most glaciers in Svalbard, Høganesbreen is presumably polythermal (Reference LiestølLiestøl, 1977). The large ice-cored moraines in front of Høganesbreen and the partially debris-covered glacier tongue suggest a cold base. Icings are found regularly in front of Høganesbreen. Reference LiestølLiestøl (1977) associated such icings with subglacial water drainage underneath the cold-based tongue. Previous studies have shown that subglacial drainage systems can exist underneath cold-based glaciers (Dowdes-well and Drewry, 1989; Reference Pfirman and SolheimPfirman and Solheim, 1989; Reference Spring and HutterStone and Clarke, 1993; Reference Wadham, Hodson, Tranter and DowdeswellWadham and others, 1998).
Fig. 1. Topographic map of Høganesbreen (grid: UTM zone 33, unit: metre), with GPRprofile (black lines), start- and end-points of drainage tunnel in bedrock (black stars), transport road up to the mine entrance (dashed shaded line) and mine (hatched area).
Fig. 2. Temperature profile of a through-glacier borehole (140 m) in the vicinity of the assumed outlet of the drainage tunnel in bedrock. The thermistor chain ends at 120 m depth.
Most of the Svea Nord coal mine is situated in the bedrock beneath Gruvefonna (a small ice cap of 8 km2 ranging up to 900ma.s.l.). The ground-water enters the mine as coal production extends beyond the permafrost areas into a talik underneath the ice cap. SNSK observes an inflow of approximately 0.03 m3 s–1 in late winter, which increases to1.38 m3 s–1 during summer. This indicates a possible connection between surface meltwater from the accumulation area and the ground-water entering the mine.
The local geology is dominated by gently dipping (north-east–southwest) layers of Tertiary sediment rocks (sandstone, siltstone, shales and coal). The thickness of the coal layer ranges from 3 to 5.5 m. SNSK has carried out extensive surveys in Svea Nord. Several core drillings into the bedrock have shown that Høganesbreen has eroded through the coal layers, where the bedrock is lower then 300–340ma.s.l.
To reduce the cost of pumping water out of the mine, SNSK has suggested drilling a drainage tunnel through the bedrock to the glacier bed of Høganesbreen (Fig.1). The formation of a subglacial conduit connected to the end of the tunnel is then initiated by pumping heated water to the glacier bed. Such a conduit might permit gravitational drainage of water. In this study, we examine the requirements for maintaining such a drainage channel in open-channel conditions throughout the year. This requires an understanding of the subglacial plumbing system, which is a complicated affair (Reference PatersonPaterson, 1994) and demands a huge amount of field data and theorizing. In this paper, we present field data, which are required in order to define the boundary conditions of a numerical model. This model is then used to analyze the hydraulic conditions of a possible drainage route.
The objectives of this paper can be summarized as follows:
(1) Determine the topography of the glacier bed by mapping the ice thickness of Høganesbreen.
(2) Highlight possible subglacial drainage pathways based on a hydraulic-potential approach, derived from digital elevation models (DEMs) of the glacier bed and the surface.
(3) Assess the pressure variations in a subglacial channel at a given discharge by model calculations.
Field Data
The thickness of Høganesbreen was already known,to some extent, as a result of through-glacier geologic altest drillings conducted by SNSK and a ground-penetrating radar (GPR) profile driven along the medial moraine between Høganesbreen and Stollbreen (personal communication from T. Larsen, 1998). However, neither the geological drillings nor the GPR profileprovidedenough data to assess the possibility of draining the water from the mine subglacially. For this purpose, we performed additional measurements of ice thickness and temperature in the vicinity of the planned junction of the drainage tunnel and the subglacial channel.
Ice thickness and bedrock elevation attained with GPR
The radar system used was a Ramac/GPR made at Malå GeoScience, Sweden, equipped with a 25 and 50 MHz antenna. Receiver and transmitter arms were oriented in line (E-plane). Noise reduction was achieved by averaging the received signal eight times before recording. The radar antenna was mounted 2 m behind a sledge carrying global positioning system (GPS) and GPR units. The GPR data were sampled with a fixed time interval of 0.5 s. The measurement interval along a transect is dependent on the snowmobile speed (~4 m s–1) and averages ~2 m.
Borehole measurements
A temperature profile through the glacier in the vicinity of the planned drainage tunnel was measured by a thermistor chain lowered down a borehole. The temperature measurements were carried out using individual calibrated thermistors with electrical resistance of about 32 kΩ at a temperature of 0°C. The accuracy of the temperature measurements is estimated to be ±0.2°C. The total depth of the borehole is 152 m and ends, most probably, at the glacier bed. However, the GPR profile indicates a glacier thickness of 135–140 m. Hence, the borehole must be inclined. Further, the length of the thermistor cable is only 130 m and consists of 13 thermistors, which leaves a gap of 22 m at the bottom of the borehole, or 20 m above the glacier bed.
Model Theory
Hydraulic potential
To investigate the possible structure of the subglacial drainage system, we followed the widely applied hydraulic-potential approach (e.g. Reference SharpSharp and others, 1993; Reference Hagen, Etzelmuller and NuttallHagen and others, 2000). The hydraulic potential H, expressed as the height of a water column, is given by
where z is the topographical potential of the glacier bed and h is the hydraulic head. In the long-term scale, the value of h varies between 0 and the ice-flotation level h* such that Equation (1) can be written
The factor k indicates the flow conditions, which are allowed to vary between open-channel (k = 0) and flotation (k = 1).
Water flow through a subglacial conduit
To assess the possibility of draining the water gravitationally through a subglacial conduit connecting the mine to the glacier terminus, we applied a model that was formulated by Reference ClarkeClarke (1996) and slightly modified by Reference SchulerSchuler (2002). With this model, we simulated the distribution of water pressure along a subglacial conduit, and its transient evolution in response to a time-varying discharge input.
The lack of detailed information on the real course and geometry of a subglacial conduit requires that we make the simplest possible assumptions. Hence, we assume that the conduit follows a straight line and has a circular cross-section. Additionally, we consider a uniform distribution of pressure within the conduit cross-section and treat the problem in one dimension, the direction of flow. Furthermore, the formulation implies that the energy released by the water flow is consumed locally and the contribution to discharge from melting the conduit walls is neglected.
Using these assumptions, we follow Reference Spring and HutterSpring and Hutter (1982) and use the Darcy–Weisbach equation to describe the loss of hydraulic potential ΔH. For turbulent flow through a pressurized conduit element having length Δx and cross-sectional area A, it is stated that
where v is the mean flow velocity of the water, g is acceleration due to gravity and f is a dimensionless friction factor. Using Equation (1) and v = Q/A yields an expression for the loss of pressure head Δh:
Where Q exceeds the maximum amount that can be drained gravitationally through a conduit of given size, flow conditions are pressurized. Below that limit, open-channel flow takes place and hence h = 0. The threshold Qmax for the transition from pressurized to open-channel flow is approximated using Equation (4) for Δh = 0:
To complete the model, we need a relation that describes the evolution of the conduit size A. Reference RöthlisbergerRöthlisberger (1972) suggested that a subglacial conduit can adjust its size to the prevailing hydraulic conditions. He described a steady state by balancing two opposing mechanisms: melt enlargement due to dissipation of energy by the water flow, and creep closure of the viscous ice. Reference NyeNye (1976) and Reference Spring and HutterSpring and Hutter (1982) developed a time-dependent theory, on which our simplified model is based. The evolution of A with time is described by
where the first term on the righthand side expresses the melt enlargement and the second term describes creep closure of a cylinder surrounded by ice (Reference NyeNye, 1953). Here, p w and P i denote the water and ice density, respectively. Reference RöthlisbergerRöthlisberger (1972) introduced the dimensionless constant γ = ctcwPw, where ct expresses the change in melting temperature per unit pressure and cw is the specific heat of water. B and n are parameters in Glen’s flow law (Reference NyeNye, 1953). If the conduit is not filled by water, the melt enlargement is reduced by a factor C, which represents the ratio to which the conduit perimeter is wetted, though C = 1 for pressurized flow (Reference CutlerCutler, 1998). For further details, we refer the reader to Reference SchulerSchuler (2002).
Application of the model
The calculation of the hydraulic potential and associated subglacial flow pattern is based on the surface and bedrock DEMs. We used the spatial analysis capabilities of a geographic information system (GIS) to reveal the basic hydrological characteristics present beneath Høganesbreen. The procedure closely follows that described by Reference Etzelmuller and BjornssonEtzelmuller and Björnsson (2000) and Reference Hagen, Etzelmuller and NuttallHagen and others (2000). A short description is given here. First, the hydraulic potential was calculated applying Equation (2). We considered two situations: k = 1 and k = 0, characteristic limits of possible pressure variations. In the first case, k = 1, it is assumed that the water pressure is everywhere equal to the ice overburden pressure. In this case, the direction of water flow will mainly be controlled by the surface gradient. In the second case, k = 0, water flows in an unpressurized, open channel and the water flow follows the topographic gradient of the bed.
By mapping the hydraulic-potential surface and assuming that the water flow is always directed down the maximum potential gradient (Reference ShreveShreve, 1972), the possible location of subglacial flow paths can be derived. For our purpose, we used the contribution area approach in which, for each individual gridcell, the total is built of all gridcells in the upstream direction that drain to the considered gridcell. Based on the resulting potential map, the local flow direction is determined by adaptive filtering techniques, namely the D8 algorithm (Reference O’Callaghan and MarkO’Callaghan and Mark, 1984). This algorithm determines the direction of the flow and quantifies the area that contributes to the considered gridcell. Preferential drainage pathways can be defined as those areas for which the contributing area exceeds a defined threshold. Furthermore, this approach describes a static situation and does not take into account changes of mass balance or water pressure in channels. For a more detailed presentation of the method, we refer to Reference Etzelmuller and BjornssonEtzelmuller and Björnsson (2000) and Reference Hagen, Etzelmuller and NuttallHagen and others (2000).
Equations (4) and (6) have been discretized on a numerical grid with Δx = 100 m for which the bedrock elevation z and ice thickness h* are prescribed along the 4100 m long conduit trajectory shown in Figure 3b (k = 1). For our calculations, we assume that the conduit is straight and we neglect any influence of curvature on the water flow. With the boundary condition h = 0 at the glacier terminus (x = 0), the model equations can be solved numerically using a finite-difference scheme (Reference ClarkeClarke, 1996). Equation (5) was used to discriminate between open-channel and pressurized conditions. Where flow was pressurized, the set of Equations (4) and (6) was solved simultaneously, otherwise h was set equal to 0, and only Equation (6) was solved. The results were calculated with a temporal resolution of 1 hour.
Fig. 3. Predicted flow pathways beneath Høganesbreen for (a) unpressurized (k = 0) and (b) pressurized (k=1) conditions (Equation (2)). Stars indicate the start-and endpoints of the drainage tunnel in bedrock. In (b), a grey zone is added around the 4100 m long conduit trajectory used in the numerical modeling (Equations (4) and (6) ).
The discharge input Q was described in close relation to observations made in the mine. Q varies from 0.03 m3 s–1 during winter to 1.38 m3 s–1 during summer. The transition from winter to summer conditions is described by a linear increase over 10 days, whereas the termination of the 40 day ablation period is abrupt.
In our computations, we had to assign values for several parameters. The parameter B in Glen’s flow law was varied from 2.4 x 10–24 to 6.8 x 10–24 s–1 Pa–3 to represent ice of different stiffness (Reference PatersonPaterson, 1994). Further, we used n = 3, g = 9.8 ms–2, ρ w= 1000 kgm–3, P i = 900 kgm–3, γ = 0.313 (Reference ClarkeClarke, 1996) and f =05. The latter value was used by Reference Spring and HutterSpring and Hutter (1982) and lies within the narrow range of roughness parameters determined for subglacial conduits by other authors (Reference NyeNye, 1976; Reference ClarkeClarke, 1982; Reference Spring and HutterSpring and Hutter, 1982; Reference NgNg, 1998). Nevertheless, this value is very large, and, intuitively, one would expect an ice-walled conduit to be much smoother. However, we have to bear in mind that in our model formulation we have neglected any explicit description of curvature of the conduit axis and irregularities of the cross-sectional shape and size. These effects would cause local head losses in addition to frictional head losses. Here, we use a single parameter f to describe frictional, as well as local, head losses, so the value of f can be significantly larger than if it were to describe frictional losses only.
Results
Surface and bed topography
An area of 3 km2 was surveyed with GPR/GPS during the fieldwork in September 2001. GPR traces were recorded at 20000 locations (Fig. 1). Positions of individual radar measurements were determined by differential GPS to within ±1m. Ice thickness was determined at all trace locations using the measured two-wave travel time and assuming a constant velocity through the ice of 169 m μs– 1 . The bedrock topography was determined by subtracting the ice thickness from the surface elevation, which was derived from the GPS survey. A DEM of the glacier bed with 10 m grid spacing was produced by interpolating the dataset combining GPR data, surrounding bedrock elevations along the glacier margins and bedrock elevation derived from the SNSK core drillings.
A DEM of the surface topography was derived from a digital map of scale 1:10 000 with 10 m contour lines (© SNSK), made from air photographs taken in August 1990. Gridding to 10 m resolution was carried out following a Reference HutchinsonHutchinson (1989) procedure.
Field observations
No water-level changes were observed in the borehole after the assumed glacier bottom was reached. After 3 days the water was slowly squeezed out, which normally is interpreted as the freezing of the water column surrounded by subfreezing ice. Since no obvious variations were discovered, the borehole most probably did not connect to any existing drainage system.
The temperature profile (Fig. 2) measured at Høganesbreen is similar to other ablation-zone temperature measurements carried out on Svalbard glaciers (Reference Ødegård, Hagen and HamranØdegård and others, 1997). The internal temperature distribution increases with depth, but the temperature never reaches the pressure-melting point. Further, our temperature profile shows that the glacier has a temperature of–1.3°C at 120m depth, 20m above the bottom. The observed temperature gradient indicates that pressure melting requires that the ice must be at least 180 m thick. It is therefore unlikely that the glacier bed is temperate in the area where the mine outlet conduit should connect to the glacier bed. Two earlier temperature records higher up on Høganesbreen (500 m a.s.l.) show a cold surface layer of–2.5°C at 80 m depth (personal communication from J. O. Hagen, 1988). Our temperature record shows –2.0°C at 80 m. Hagen’s records only cover two-thirds of the glacier thickness, and a bottom temperature is not available. Bjornsson and others (1996) showed that GPR internal echoes occur in the transition zone between warm and cold ice. Such echoes were found in the thickest parts of Høganesbreen, but were not detected in the area of interest.
Hydraulic potential mapping
The k value (k € [0,1]) can be tuned to mimic the evolution of different drainage situations. Figure 3a and b show the results of two different simulations. For k = 0 the water channels follow the bottom topography. k = 0 indicates one main stream 150 m west of the junction between the drainage tunnel and the glacier bed, whereas k = 1 indicates that the water from up-glacier (relative to the drainage tunnel) drains into two rivers: one east, and the other west, of the intersection of the glacier bed and the drainage tunnel. For the pressurized case, k = 1, the model does not take into account the fact that channels of different pressures will converge to the one with the lowest pressure. Both simulations indicate that a main drainage pathway is located 100–150 m away from the outlet point and not in the direct vicinity of the mine outlet.
One should bear in mind that, for both simulations, only the major channels are shown, because the area contributing to each individual cell has to exceed a defined threshold to be plotted in Figure 3. Furthermore, both maps represent steady-state situations and do not take in account any temporal or spatial pressure variation in the glacier hydrological system.
Computation of subglacial water pressure
We started the model calculations by assuming steady state at Q = 0.03 m3 s–1 (Fig. 4a). The resulting pressure at the place where the mine outlet connects to the conduit is 5 bar for the stiffer ice and 8 bar for the softer ice (Fig. 4b and c). At the transition to summer discharge, the calculated pressures in both scenarios jump to unrealistically high values before they approach a new steady-state value. Being 4– 6 bar, the new equilibrium pressures are lower than the winter values. With the sudden termination of summer discharge, the flow regime switches to open-channel conditions, and the cross-sectional area of the conduit experiences increased closure rates. With the progressive decrease of the conduit size, the flow regime reverts to pressurized conditions, and water pressure increases until a new steady-state size is reached. In the scenario calculated for stiff ice, this adjustment requires about 100 days, whereas 50 days are required in the soft-ice scenario (Fig.4b and c).
Fig. 4. (a) Discharge input; (b, c) time series of water pressure (solid line) and conduit cross-sectional area (dashed line) as calculated using two different values of the flow-law parameter B: B= 6.8 x 10–24s–1Pa–3 (b) and B = 2.4 x10–24s–1Pa–3(c).
Concluding Discussion
We have used known field techniques and found Høganesbreen to be below the pressure-melting point through the entire thickness in the research area. Using numerical models, we investigated the possibilities of draining water from a mine along the base of a glacier. The results of the hydraulic-potential approach show that the existence of a channelized drainage system in the direct vicinity of the outlet point from the mine is unlikely. This finding suggests that water at the site flows, rather, through a spatially distributed and hydraulically inefficient system. However, 100–150 m away from the outlet point, a preferential flow path might exist. A connection of the mine outlet to this channel might be accomplished by melting an artificial channel using a circulation pump connected to a heater. In the following, we discuss the results of modeling the behaviour of such a conduit.
Our goal was to assess the possibility of maintaining an artificial subglacial conduit through the period of low discharge. Therefore, we simulated a seasonally varying discharge through a subglacial conduit. The discharge input function (Fig. 4a) was selected to mimic the seasonality of melt in Svalbard, with ablation occurring over 24 h per day. Using a simplified geometry of Høganesbreen, we calculated the evolution of water pressure at the junction between the mine and the glacier bed in response to the discharge input (Fig. 4).
It is not possible for the model to take into account storage of water during periods when the water pressure reaches the ice overburden pressure. Therefore, the resulting high pressure values in this situation are unrealistic. Further, it can be argued that the sudden drop in summer discharge is unrealistic. However, it represents the situation in which the duration of open-channel conditions is maximal. If the decline of discharge is slower, the closure of the conduit required to obtain pressurized conditions will be smaller and the conduit will be pressurized earlier or even throughout the period.
The mine has a higher topographic potential of 10 m relative to the place where the outlet conduit should connect to the glacier bed. This corresponds to a back-pressure of 1bar. To drain the water out of the mine, the pressure in the subglacial conduit must be 5 1 bar, otherwise water would back up in the mine. Modeling results show that both maximum and minimum values of Q reveal pressures substantially higher than this limit. The results indicate further that the closure rate of even stiff ice would be too rapid to permit open-channel conditions for a period longer than a few weeks. On the other hand, the water pressure at the glacier bed is likely to be much higher than the pressure in the mine, and subglacial water would flow into the mine. Drainage systems different from a conduit system usually operate at substantially higher pressures (Reference Röthlisberger, Lang, Gurnell and ClarkRöthlisberger and Lang, 1987; Reference Fountain and WalderFountain and Walder, 1998), which would add to this effect.
We have advised SNSK against proceeding with plans for draining the mine subglacially. The results of this study show that the planned drainage tunnel reaches the glacier bed where a subglacial channel is unlikely to exist. Furthermore, if an artificial subglacial drainage system could be established, it would be pressurized most of the time. The water-pressure gradient from the mine to the glacier bed would be negative most of the time. Therefore, rather than drain water away from the mine, an artificial conduit would flood it, thus compounding the mining company’s problem.
Acknowledgements
This project was initiated and supported by SNSK, who also provided the logistics during the fieldwork. Thanks to M. Vaksdal and T. Chareyron who helped with the fieldwork, and T. Abrahamsen and J. Stenvold at SNSK for all the support in Sveagruva. The authors wish to thank B. Etzelmuller for providing computer code for the rapid calculation of hydraulic potential. We gratefully acknowledge the comments made by F. Ng and an anonymous referee, which helped to clarify an earlier draft of the manuscript.
