Controls on the basal water pressure in subglacial channels near the margin of the Greenland ice sheet

Abstract

Assuming a channelized drainage system in steady state, we investigate the influence of enhanced surface melting on the water pressure in subglacial channels, compared to that of changes in conduit geometry, ice rheology and catchment variations. The analysis is carried out for a specific part of the western Greenland ice-sheet margin between 66° N and 66°30′N using new high-resolution digital elevation models of the subglacial topography and the ice-sheet surface, based on an airborne ice-penetrating radar survey in 2003 and satellite repeat-track interferometric synthetic aperture radar analysis of European Remote-sensing Satellite 1 and 2 (ERS-1/-2) imagery, respectively. The water pressure is calculated up-glacier along a likely subglacial channel at distances of 1, 5 and 9 km from the outlet at the ice margin, using a modified version of Röthlisberger’s equation. Our results show that for the margin of the western Greenland ice sheet, the water pressure in subglacial channels is not sensitive to realistic variations in catchment size and mean surface water input compared to small changes in conduit geometry and ice rheology.

Introduction

The study of basal drainage beneath glaciers and ice sheets has been pursued from a theoretical and experimental point of view for many years, with selected sites ranging from well-studied mountain glaciers (e.g. Hock and Hooke, 1993; Iken and Truffer, 1997; Arnold and others, 1998) over intermediate-sized ice caps such as Vatnajökull, Iceland, (e.g. Björnsson, 1982; Flowers and others, 2003) to Antarctic ice streams (e.g. Alley and Bindschadler, 2001). In contrast, relatively few studies of basal drainage have focused on the Greenland ice sheet (Thomsen and others, 1988; Iken and others, 1993; Ahlstrøm and others, 2002), although interest has increased with the recent investigations of Krabill and others (2000) and Thomas and others (2003) showing that the Greenland ice-sheet margin is, in many places, thinning at a rate which cannot be explained solely by an increase in the surface melt. Changes in ice dynamics must account for at least part of the thinning and it is likely that an increase in the surface melt might influence the dynamics by changing the basal conditions. The type of basal drainage system that might be active underneath the Greenland ice sheet is largely unknown, as the methods usually applied in similar studies of mountain glaciers have been either unsuccessful or not feasible. Dye-tracing studies have been attempted as a way of delineating the hydrological basins of the Greenland ice sheet, with moderate success for smaller well-defined marginal basins, but experiments aiming to delineate the larger basins have so far failed (personal communication from H.H. Thomsen, 2003).

Our aim in this study is to reveal how sensitive the water pressure in subglacial channels near the ice margin is to changes in surface water input compared to changes in channel geometry, ice rheology, channel sinuosity and catchment area, respectively, when applied to real ice-sheet geometry. To accomplish this, we utilize new surface and bed elevation data from the western margin of the Greenland ice sheet (see Fig. 1) to locate the likely positions of subglacial drainage channels and their contributing basins using Björnsson’s (1982) simplification of Shreve’s (1972) theory. We assume that the derived drainage system near the ice margin is channelized in the form of an arborescent network of pressurized channels, varying in cross-sectional geometry from semicircular to broad and low (termed H channels by Hubbard and Nienow, 1997), and that the ice sheet rests on hard bedrock as opposed to a layer of till. The theoretical considerations of Röthlisberger (1972) showed that for a subglacial channel in steady state, higher discharge would be equivalent to a lower water pressure, thus favouring larger channels over smaller ones. Ice viscosity and channel roughness was shown to be several orders of magnitude more important than channel discharge in determining the water pressure. However, Röthlisberger (1972) also pointed out that the discharge of a channel may regain importance from its extreme variability. We have taken up this point and we examine the importance of discharge variations on water pressure due to climatically induced changes in surface water input and catchment variations. The water pressure is calculated in a specific subglacial channel using a slightly extended version of the solution of Hooke and others (1990) taking into account a range of possible channel geometries and the resistance of basal drag to channel closure.

Fig. 1. Overview of the region studied. The grey lines on the ice sheet show the extent of the airborne survey carried out in 2003, and the black lines signify which part was utilized for this work. The point marked A shows the position of a weather station on the ice sheet and also marks the outlet of the hypothetical drainage channel used for calculations of basal water pressure.

Elevation Models

One of the major obstacles to the study of basal drainage has been the lack of the necessary detailed topographical information on the terrain above as well as beneath the ice sheet. Thorning and Hansen (1987) used an ice-penetrating radar carried on a helicopter for a study near the ice margin, using visible features of the landscape to determine the position. However, large-scale regional studies of sufficient precision as envisaged in this work were not feasible until accurate positioning was made possible with the global positioning system (GPS). The investigated region was surveyed by aircraft in 2003, as shown by grey and black lines in Figure 1. Data for the two short lines oriented from southwest to northeast shown at the bottom of Figure 1 were supplied by W. Krabill. Surface altitude and ice thickness were measured by combining laser altimetry and ice-penetrating radar with measurements from two differential GPS instruments and an inertial navigation system (INS) (Christensen and others, 2000; Ahlstrøm and others, 2002). In the present study, only the bed elevation derived from the ice-penetrating radar was used. This dataset was merged with bed elevation data derived from combining the surface elevation model described below with ice-thickness data (Gogineni and others, 2001) collected in the area between 1993 and 2001 (the two short black lines of Figure 1 oriented from southwest to northeast). Only data from the black flight-lines shown in Figure 1 were used. The spacing of the flight-lines was chosen to be roughly 2.5 km to enable resolution of the typical undulations of the bed as derived from a twodimensional Fourier analysis of the neighbouring ice-free terrain. In addition, the position of the ice margin was imposed on the bed dataset to force the bed elevation to match the known elevation of the adjacent ice-free terrain. Subsequently, the merged bed elevation dataset was averaged over 250 m and gridded using ordinary block kriging, a linear Gaussian least-squares interpolation method (Carr, 1995; Flowers and Clarke, 1999). To facilitate the kriging, directional semivariograms were calculated in eight directions to make use of the a priori knowledge contained in the dominant topographical features of the bed. Semivariograms describe the covariance of a dataset with itself as a function of the lag distance between pairs of data points. The resolution of the interpolated bed elevation grid was chosen to be 250 m, a compromise between the very high along-track resolution and the 2.5 km across-track resolution.

The surface elevation model was based on repeat-track interferometric synthetic aperture radar (InSAR) analysis of two image pairs from the European Remote-sensing Satellite 1 and 2 (ERS-1/-2) tandem mission. The ERS-1 images were acquired on 20 October and 29 December 1995. The derivation of the elevation model is described in Ahlstrøm and others (2002) and will not be discussed further here. The precision was found to be within 10 m rms of the elevations derived from the airborne laser altimetry survey in 2000. The original resolution of 100 m was reduced to match the resolution of the bed elevation model. Contour maps of the bed and surface elevation models are shown in Figure 2.

Fig. 2. Contour maps of the bed and surface elevation models. See text for details.

Conduit Localization

Our aim is to localize a representative subglacial channel draining a significant ice-sheet surface area. To find the most likely positions of subglacial drainage channels within the study area, Shreve’s (1972) formulation of the englacial water routing was used with the simplifying assumption of Björnsson (1982) that all water at the surface reaches the bed through moulins and crevasses, and drains along the base of the ice sheet which is assumed to be impermeable. The assumption that the surface water reaches the bed at the place where it is formed is supported by a map of individual moulin catchments drawn from aerial photographs over the Greenland ice-sheet margin around 69°30′ N, 50°00′W by Thomsen and others (1988) showing that the mean surface flow distance is on the order of a few kilometres.

The above-mentioned assumptions allow the direction of the water flow at the base of the ice sheet to be determined by a potential Φ_b defined as

(1)

where ρ _w and ρ _i are the densities of water and ice, respectively, Z _b is the bed elevation, Z _s the elevation of the ice-sheet surface and g the gravitational acceleration. The term k ρ _i g (Z _s – Z _b) in Equation (1) is the subglacial water pressure which is proportional to the overburden ice pressure, with the factor k = 0 corresponding to atmospheric pressure in subglacial channels, and k = 1 corresponding to subglacial water pressure equalling overburden ice pressure. The water at the base of the glacier will flow in the direction of the maximum gradient of the potential Φ _b. This knowledge makes it possible to derive the most likely positions of subglacial channels from a surface and bed elevation model and an estimate of the k factor (Hagen and others, 2000). Just as it is possible to derive likely positions of individual channels from the potential Φ _b, it is also possible to derive the total hydrological drainage basin fed by these channels. The basin extent, which naturally has a strong impact on the discharge Q through an individual channel, is therefore dependent on k.

Initially, the probable channel network was derived for values of k ranging from 0 to 1.05. This exercise revealed that one particular channel outlet, situated underneath the station marked A in Figure 1, drained the largest area for all values of k. Using this outlet, the contributing basin was calculated for each value of k in steps of 0.05. These basins can subsequently be used to estimate the water discharge through the channel outlet into the Tasersiaq lake system as a function of the k factor once the surface melt has been determined over the region.

Obviously, the quality of the elevation models is important for the basin derivation. In studies in the same region, Ahlstrøm and others (2002) and Ahlstrøm (2003) found a strong variation of the basin size was mainly due to artificially introduced topographical features in the bed elevation model. The erroneous features arose from merging two bed elevation models of different quality and led to an extension of the planned survey of 2003 with the two southernmost northwest–southeast-oriented flight-lines seen in Figure 1 and the inclusion of data from Gogineni and others (2001) to cover the critical region. Ahlstrøm and others (2002) also showed that completely different basins were derived using existing elevation models of lower resolution, based on a sparser dataset, compared to the highresolution models derived from the dense datasets of regional airborne surveys.

The simplified theory of Equation (1) implies that the ratio of conduit closure to conduit radius is dependent on ice thickness to the third power (assuming Glen’s flow law with an exponent n = 3). However, the change in the derived river network with k is minor for k > 0.5, and the position of the major outlet is nearly constant for all k, suggesting that little would be gained from introducing a more elaborate expression for the dependence of the closure/radius ratio on ice thickness.

The surface and bed elevation profiles along the main channel were extracted and used in the subsequent calculation of channel water pressure for all values of k. Only the part of the channel within 10 km of the ice-sheet margin was used, partly to limit the analysis to the main channel and partly because the assumption of channelized drainage is most likely to be fulfilled near the margin.

Water Pressure in Conduits

Hooke and others (1990) presented a modified version of Röthlisberger’s (1972) steady-state relationship between water pressure in a conduit P _w, overburden pressure of ice P _i, slope of the conduit β, water discharge Q, Glen’s flow-law parameter B, conduit closure by ice flow and wall roughness, expressed along a horizontal x axis pointing against the flow of water in the conduit. In this work we use the modified relation of Hooke and others (1990) with an additional modification to include the effect of basal drag τ _b on conduit closure at the bed. Basal drag is included by formulating the effective stress as This means that setting the basal drag to zero reduces the relation to the solution of Hooke and others (1990). The basal drag τ _b is set to 50 kPa in the present work (Paterson, 1994). Finally, we arrive at the following modified version of Röthlisberger’s (1972) relationship:

(2)

where

(3)

and

(4)

where L is the latent heat of fusion. This relation is based on a balance between wall melting and closure by ice flow and relies on the Manning equation to describe the relation between Q and tunnel roughness. The latter is expressed through the Manning roughness coefficient k. The modification of Hooke and others (1990) compared to Röthlisberger (1972) lies in Ω which is a conduit-shape factor. The conduit shape is approximated by the space between the chord of a circle and the arc subtended by that chord (Hooke and others, 1990). Modifying the relation for conduit closure used by Röthlisberger (1972), Hooke and others (1990) arrived at a relation between the conduit-shape factor Ω and the angle Ə subtended by the chord of the circle spanning the conduit. A semicircular conduit thus corresponds to Ω = 180°, whereas the conduit approaches a thin film for Ω near zero. The height-to-width ratio of the tunnel thus influences the water pressure in the conduit directly.

The analysis above only allows for steady-state solutions of a channelized drainage system, meaning that the temporal evolution of the drainage system is not captured. The present work is thus restricted to a comparison of the influence of various parameters on a theoretical steady-state situation. To best accommodate the assumption of a fully developed channelized drainage system near equilibrium, the investigated period should preferably be during, or just after, the peak of the ablation season. To allow the use of in situ observations for the calculation of the meltwater discharge Q, July 1999 was chosen. The steady-state situation is here assumed to exist during July 1999, ignoring sub-diurnal peaks and even fluctuations over several days in the water flux. The present work thus assumes a temporal scale of 1 month in an attempt to model a mean equilibrium state during the warmest part of the ablation season, even if such an equilibrium state may never be reached in nature.

Meltwater Discharge

The meltwater discharge Q in Equation (2) has been calculated as the surface meltwater from all the upstream grid elements contributing to the conduit chosen for the analysis. The production of surface meltwater has been modelled using the close relationship between positive degree-days and surface melting (Ohmura, 2001). This approach requires a gridded temperature field which can be adjusted to allow for the influence of climate change. To realistically model the effect of a rise in regional air temperature on meltwater discharge, the dampening effect of the zero-degree surface of the ice sheet on the screen-level air temperature over the ice should be taken into account. Thus, the assumed air-temperature rise was imposed on measurements from an automatic weather station (AWS) located within the Tasersiaq lake basin (station 105 in Fig. 1) but away from the ice sheet. Subsequently, the empirical transfer function of Ahlstrøm (2003) was used to convert the air temperature from station 105, off the ice sheet, to station A on the ice sheet. Finally, an ice-sheet July temperature lapse rate of 0.0051 Km⁻¹ (personal communication from J. Box, 2001) was applied to obtain a gridded temperature field. Converting the temperatures into surface melt requires a suitable choice of degree-day factor. Data from an AWS located on the ice-sheet surface in 1999 (25 June–5 August) at 886 m a.s.l. directly above the channel (marked A in Fig. 1) indicate a degree-day factor for ice (DDF_ice) of 0.014 m per positive degree-day. This value is high compared to most studies (Hock, 2003) but is comparable in magnitude with other observations near the same altitude on the western margin of the Greenland ice sheet during the melt season, by Van de Wal (1992) who obtained a DDF_ice of 0.020 at 1028 m a.s.l. and Ambach and others (1988) (cited by Hock, 2003) who obtained a DDF_ice of 0.0186 at 1013 m a.s.l. Thus, the observed degree-day factor of 0.014 m per positive degree-day was used to calculate the surface melt over the region from the modelled temperature field.

The meltwater discharge Q through the chosen channel can now be derived by combining the calculated surface melt with the channel catchments previously calculated for various values of k. The meltwater discharge then becomes a function of k and the assumed temperature forcing. The discharge is shown as a function of k in Figure 3. The impact of the temperature forcing intended to mimic climate change is shown in Figure 4 for k = 0.7.

Fig. 3. The channel meltwater discharge Q as a function of the k factor in Equation (1). The k factor causes changes in the drainage basin extent, and thus indirectly in the meltwater discharge Q.

Fig. 4. The channel meltwater discharge for the basin corresponding to a k factor of 0.7 as a function of the 3 m level air-temperature forcing off the ice sheet.

Results

As previously mentioned, only the part of the channel within 10 km of the ice margin was included in the model runs. Specifically, the 10 km is measured from a point at a fixed short distance from the ice margin, where the pressure in the channel is assumed to be atmospheric, and up-glacier along the subglacial channel.

For illustration of the effect of the distance from the terminus, points at a distance of 1, 5 and 9 km up-glacier were selected. Many parameters might influence the water pressure in the subglacial channels, but to maintain an overview, each parameter is varied individually relative to a reference run with parameters as given in Table 1. The change in water pressure resulting from varying a number of parameters is subsequently shown in Figure 5 as the fraction of the subglacial water pressure needed to cause flotation of the ice sheet.

Table 1. Model parameters used in the reference run

Fig. 5. (a–e) The influence of various parameters on the water pressure at three distances from the ice margin in a selected subglacial channel. The legend in (a) is valid for all plots. See text for a detailed explanation of the plots. (f) Calculated mean water speed in the subglacial channel as a function of channel geometry.

Figure 5a shows how the pressure varies as B ranges from 36 to 200 kPaa^1/3. There are many factors influencing the flow-law parameter, B, such as ice temperature, water content, impurities and crystal size. The ice temperature is unknown in the region studied, but Thomsen and Olesen (1990) found a basal ice temperature of – 0.1 °C at an ice depth of 350 m, 9.5 km from the ice margin in the Pâkitsoq area of the Greenland ice sheet (69°30′N, 50°00′W). The area studied by Thomsen and Olesen (1990) had an active basal drainage system as indicated by fluctuations in the water pressure of several boreholes. The area studied here most likely derives its basal ice from warmer areas at lower altitude compared to the area studied by Thomsen and Olesen (1990). It is thus reasonable to expect that the basal ice temperature is above –1°C and that it is possibly temperate. An ice temperature of –1°C corresponds roughly to B = 200 kPa a^1/3, whereas temperate ice would correspond to B = 167 kPa a^1/3 (Paterson, 1994). A further lowering of B could be due to higher water content, which is known to have a strong impact on viscosity. Even though the temperature of the basal ice at Dye 3 upstream of the study area was found to be around –13°C by Dahl-Jensen and Johnsen (1986), the effect might still be important closer to the margin, as suggested by modelling (e.g. Huybrechts and others, 1991). An important influence on the rate factor B of the basal ice could arise from the impurity content and crystal size of Wisconsinan/Weichselian ice compared to Holocene ice as suggested by Dahl-Jensen and Gundestrup (1987). These authors found that the rate factor B decreased to <70% of its expected value at Dye 3, corresponding to a drop in B from 167 kPa a^1/3 to 116 kPa a^1/3 due to this effect alone. The basal ice in the region studied could also be Wisconsinan/Weichselian ice unless extensive bottom melting takes place.

Figure 5b shows how pressure varies with k from Equation (1). As previously discussed, the discharge Q depends on the extent of the channel catchment, which in turn depends on k (see Fig. 3). The variation in pressure in Figure 5b is therefore due to the variation in the discharge Q. The pressure is seen to be rather insensitive to variations in the discharge Q, implying that even large variations in the basin extent due to a varying k would have a minor impact on the subglacial water pressure for a steady-state situation. Variations in k could arise from variations in surface water input or changes in, or development of, the subglacial drainage system.

Figure 5c shows the pressure variation as a function of the air-temperature forcing imposed off the ice sheet. The air-temperature forcing causes a change in the channel meltwater discharge Q through an increase or decrease in the surface melt as shown in Figure 4. Even though the change in Q resulting from the temperature forcing is significant, with a deviation ranging from 84% to 137% of the reference run, there is hardly any change in the subglacial water pressure. Thus, a central result from Röthlisberger’s (1972) analysis, namely that higher discharge Q means lower water pressure, is evident in Figure 5c, but of little consequence.

Figure 5d illustrates the dependence of water pressure on the degree of sinuosity of the channel path along the bed, which should be small in the present case since the channel path has been derived on a 250 m grid and the assumed channel is probably rather wide. The dependence is rather significant, but it is difficult to estimate channel sinuosity. Possibly, eskers left exposed by the retreated ice sheet could provide some evidence of channel sinuosity.

Figure 5 e shows the influence of channel geometry, expressed as the ratio of channel height to channel halfwidth, on the water pressure using the approximation of Hooke and others (1990) for linking Ω to channel geometry. As discussed in Hooke and others (1990), this approximation was off by an average of 22% in estimated closure rate when compared to a numerical solution, but is considered adequate for the current purpose. The basal water pressure approaches the flotation pressure as the assumed channel geometry goes from semicircular to near-film. Figure 5f accompanies Figure 5e and shows the calculated mean water speed in the subglacial channel for a meltwater discharge of 153 m³s⁻¹ as a function of channel geometry. Generally, the mean water speed for a near-film channel geometry is reduced to approximately one-third of the semicircular case irrespective of the distance from the ice margin.

Discussion

An important assumption in the present approach is that we are only dealing with a steady-state situation. Thus, the discharge is calculated as a mean July value for a realistic range of air temperatures and channel catchments. In reality, both seasonal and diurnal discharge variations are much larger than the changes examined here. Cutler (1998) applied a transient model of subglacial tunnel evolution driven by surface water input variations calculated for Storglaciären, Sweden. He concluded that typical diurnal fluctuations played a minor role in the channel evolution and that channels responded on a time-scale of days to fluctuations in inflow. He also concluded that under transient conditions of diurnally varying discharge, the water pressure varied in phase with the discharge, but that over longer time-scales there was broad agreement with Röthlisberger’s (1972) scenario for phase differences between water pressure and discharge. This scenario states that the water pressure should be high in the spring and early summer when the drainage system expands, and then decline after the peak of the ablation season as the channel closure lags behind the decrease in discharge. Accepting Röthlisberger’s scenario and Cutler’s conclusion that channel evolution is insensitive to diurnal variations suggests that the steady-state assumption made in this work is reasonable at the peak of the ablation season. The water pressure in the channel would then fluctuate around the derived solution with diurnal variations in surface water input, as the channel would be unable to accommodate the rapid changes.

Another important assumption is that the drainage system is channelized. The evolution of the drainage system for the ice-sheet margin studied here is unknown, and difficult to infer from the limited amount of available information. The subglacial channel hypothesized in this work drains below the largest ice lobe of the region into an extensive lake system, precluding any attempt to validate the existence of this particular outlet. Satellite images show that great supraglacial meltwater lakes with diameters of several kilometres are abundant near the equilibrium-line altitude around 1400 m a.s.l. At lower altitudes, an extensive supraglacial drainage system efficiently drains the surface meltwater until it terminates in crevasses before reaching the ice-sheet margin. The crevassed region varies in extent depending on subglacial topography, but is sufficiently widespread over the margin to facilitate the transport of water from the surface to the bed through development of moulins and propagation of water-filled crevasses (Boon and Sharp, 2003). Unfortunately there is no solid evidence for moulin drainage, as the resolution of the imagery does not allow for the localization of individual moulins as done by Thomsen and others (1988) in the Pâkitsoq area further north. However, observations from aircraft confirm that most of the surficial meltwater enters the englacial drainage system before reaching the ice-sheet margin. Measurements of the change in surface velocity over the year, water-level fluctuations in boreholes connecting to the subglacial drainage system or hydrochemical analysis of water emerging from the channels would each be a valuable tool in estimating the type of drainage system and possible transitions between channelized and distributed drainage.

In the present work, we have assumed a hard bed. This assumption is not supported by any observations, and there may be a layer of till at the bed. For the case of a till-layered bed, Walder and Fowler (1994) showed that for the small surface slopes typical of ice sheets, drainage would occur in a distributed non-arborescent network of canals (N channels) incised into the till, whereas for larger slopes both canals and regular channelized drainage (R channels) could occur. This conclusion points towards the ice-sheet margin as the most likely area for development of a channelized drainage system due to the larger slope. The drainage system beneath the Greenland ice-sheet margin could also be a distributed system of linked cavities as envisaged by Kamb (1987). As discussed by Kamb (1987), the linked-cavity system can remain distributed with a high water pressure, whereas the channelized system is inherently unstable and must evolve towards fewer and larger channels with a progressively lower steady-state water pressure. In the present work, we have focused on the near-marginal region at the peak of the ablation season. Furthermore, we have examined only the largest channel, assumed to drain a large part of the region. Since a high discharge over an extended period favours a transition to a channelized drainage system (Iken and Truffer, 1997), we have selected the time and place where the channelized system is most likely to have developed. However, if the ice sheet moves over deformable till, the channelized drainage system is likely to be restricted to near the ice-sheet margin. This would imply that each channel outlet would drain a much smaller catchment and carry substantially less water than assumed here.

Conclusion

We have addressed a question posed by Röthlisberger (1972) regarding the influence of realistic discharge variations on the water pressure in subglacial channels for the case of the Greenland ice-sheet margin. To this end, we have expanded on the method of Hooke and others (1990) to evaluate the controls on the basal water pressure in subglacial channels and applied it to a realistic channel configuration beneath the Greenland ice-sheet margin assuming steady-state channelized drainage over a hard bed. The variation in the water pressure has been calculated as a function of various physical input parameters at three distances up-channel from a hypothesized major outlet. Changes in the water pressure were calculated for variations in the flow-law parameter B, cross-sectional channel geometry, channel sinuosity, discharge variations due to air-temperature changes and discharge variations due to changes in channel catchment.

The results show that the influence on subglacial water pressure of realistic changes in the mean discharge during the peak of the ablation season is insignificant compared with that of changes in the flow-law parameter, crosssectional channel geometry or sinuosity of the meandering channel.