The effect of discrete recharge by moulins and heterogeneity in flow-path efficiency at glacier beds on subglacial hydrology

Abstract

Subglacial conduit systems are thought to consist of dendritic networks that exist at lower pressure than distributed systems and have locations that are determined by theoretical hydraulic potential. On glaciers with moulins, however, meltwater is delivered to glacier beds at discrete points, violating assumptions of uniform recharge needed to calculate potential. To understand how moulins affect subglacial conduit hydrology, we used speleological techniques to map 0.4 km of subglacial conduit at the base of a moulin in Hansbreen, Svalbard, and compared our observations with theoretical predictions. The conduit began in an area predicted to lack drainage, crossed equipotential contours at oblique rather than right angles and was locally anastomotic rather than dendritic. We propose moulin locations, which are determined by the locations of supraglacial streams and crevasses, control locations of subglacial recharge. Because conduits have no direct causal relationship with gradients in effective pressure, this recharge can form conduits in areas of glacier beds that may not be predicted by hydraulic potential theory to have conduits. Recharge by moulins allows hydraulic head to increase in conduits faster and to higher values than in adjacent distributed systems, resulting in an increase rather than a decrease in glacier sliding speeds above subglacial conduits.

Introduction

Where subglacial conduits form and how they interact with distributed systems exert important but poorly understood controls on basal water pressure and ice sliding speeds (Bell, 2008). Conceptual models of subglacial conduit hydrology have been deeply influenced by early assumptions of homogeneous englacial and subglacial water flow as well as steady-state conditions (Rothlisberger, 1972;Shreve, 1972) and these studies have become fundamental components of glacier hydrological theory (Fountain and Walder, 1998; Bell, 2008;Cuffey and Paterson, 2010). According to these theories, conduits should form by the convergence of flow paths within higher-pressure distributed systems, migrate up-glacier along flow paths determined by gradients in hydraulic potential (Shreve, 1972;Flowers and Clarke, 2002;Flowers and others, 2003, 2004; Hewitt, 2011), draw down water pressure in adjacent distributed systems and ultimately decrease glacier sliding speeds (Mair and others, 2002;Bell, 2008;Schoof, 2010). Contrasting these theoretical predictions are field observations that water pressure can increase faster and to higher magnitudes in conduits than in distributed systems (Hubbard and others, 1995; Harbor and others, 1997;Gordon and others, 1998;Mair and others, 2003), causing the highest ice velocities to occur above conduits rather than above adjacent distributed systems (Harbor and others, 1997;Mair and others, 2003). Additionally, areas of the Greenland ice sheet (GrIS) that experience summer increases in velocity appear to be channelized rather than broadly distributed (Palmer and others, 2011).

Discrepancies between observations and theory suggest conceptual models using steady-state and/or homogeneous conditions may not accurately represent all of the physical processes responsible for ice acceleration. A growing body of literature suggests that glacier hydrological systems never achieve steady conditions. Importantly, these studies indicate that non-steady delivery of meltwater to glacier beds by moulins or lakes can lead to reversed gradients in the subglacial system that may promote sliding events (Anderson and others, 2003; Kessler and Anderson, 2004; Walder and others, 2006;Bartholomaus and others, 2008, 2011).

In this paper we re-examine many of the assumptions that have guided the development of glacier hydrological theory, draw on similarities between glaciers and limestone karst aquifers, recast the results of other glacier hydrology studies and present direct observations of a subglacial conduit to suggest that heterogeneity and non-steadiness in glacier hydrological systems may exert important but little-recognized controls on formation, evolution and flow in subglacial hydrological systems. We propose that discrepancies between widely cited conceptual models and observations may occur because these conceptual models have not explicitly considered the effects of heterogeneity in glacier hydrological systems.

Hydraulic potential and water flow

Water in all hydraulic systems will flow from areas of higher potential to areas of lower potential (Munson and others, 2005). In groundwater studies, hydraulic potential is usually measured directly using networks of wells. These measurements are combined with assessments of heterogeneity in aquifer hydraulic conductivity and spatio-temporal variability in aquifer recharge to make estimates of groundwater flow (Freeze and Cherry, 1979). In studies of subglacial hydrology, however, measuring water pressure in dense networks of boreholes and determining the distribution of hydraulic conductivity at glacier beds is difficult, expensive and in the case of the latter, probably impossible. As a result, glacier-wide measurements of subglacial hydraulic potential have rarely been measured directly, but instead have been estimated using the theoretical model put forth by Shreve (1972).

Shreve (1972) proposed that the hydraulic potential ($) of glacier drainage systems could be described at any point (z) in a glacier by the equation:

(1)

where Φ₀ is a reference potential, ρ_i and ρ_w are the respective densities of ice and water, g is the acceleration due to gravity and H is the elevation of the glacier surface. The second term on the right-hand side is water pressure (assumed to be equal to ice pressure) and the third term is the elevation potential. The positions of subglacial conduits in glaciers around the world have been estimated by contouring theoretical hydraulic potential from glacier surface and bed digital elevation models (DEMs) and then determining the drainage axes with the lowest continuous hydraulic potential (Sharp and others, 1993). The process is directly analogous to how surficial watersheds and stream courses are defined from topographic data. This technique has been used to reconstruct conduit networks beneath glaciers (Sharp and others, 1993;Richards and others, 1996; Palli and others, 2003;Rippin and others, 2003;Willis and others, 2009), the GrIS (Ahlstrøm and others, 2002, 2005; Mottram and others, 2009) and the Antarctic ice sheet (Fricker and Scambos, 2009).

Calculation of equipotential surfaces using Eqn (1) requires several assumptions: (1) that the glacier has some intrinsic permeability;(2) that the distribution of permeability is homogenous and isotropic;(3) that the rate of conduit enlargement and creep are exactly the same at every point within the glacier;(4) that water pressure is equal to the ice overburden pressure; and (5) that recharge to the glacier is spatially and temporally uniform (Gulley and others, 2009a). While these assumptions are probably never met on actual glaciers, the hydraulic potential method has had some success in predicting conduit locations. Hydraulic potential was used to predict the general location of a subglacial conduit on Haut Glacier d'Arolla, Switzerland (Sharp and others, 1993), although it should be noted that the hydraulic potential measured in boreholes did not match hydraulic potential calculated from Eqn (1). Water pressure in the conduit beneath Haut Glacier d'Arolla ranged from atmospheric pressure to pressures that exceeded ice overburden on both seasonal and diurnal timescales (Hubbard and others, 1995;Harbor and others, 1997). However, other attempts to use hydraulic potential to locate subglacial conduits have not been as successful. The subglacial conduit on Glacier d'Argentiere, French Alps, famously changed locations after an expensive hydroelectric project had been constructed (Hantz and Lliboutry, 1983).

In most circumstances, neither conduit network reconstructions nor estimates of subglacial hydraulic potential have been validated with glacier-wide drilling campaigns. As a result, the accuracy and precision with which Eqn (1) simultaneously predicts both the subglacial hydraulic potential as well as the locations of subglacial conduits remains largely untested. Both the change in configuration of the subglacial conduit at Glacier d'Argentiere and discrepancies between subglacial hydraulic potential predicted by Eqn (1) and subglacial hydraulic potential measured in boreholes on Haut Glacier d'Arolla suggest that Eqn (1) might sometimes predict the locations of subglacial conduits but perhaps not because gradients in ice thickness and elevation are controlling subglacial water pressure (Gulley and others, 2009a).

Two key factors governing water flow at glacier beds have received only limited attention. These factors could exert important effects on subglacial hydraulic potential that would not be identified in estimates based on Eqn (1). First, there is the issue of spatio-temporal heterogeneity in recharge, used here to refer to water entering the subglacial drainage system. Beneath glaciers, there are two primary sources of recharge: (1) the glacier sole, where basal melting can occur in response to geothermal heat or frictional heat associated with ice motion;and (2) the glacier surface, where meltwater and rain are routed to the bed via moulins. These two sources of recharge have very different distributions, magnitudes and rates and these differences result in contrasting implications for subglacial water pressure. Basal melting rates are generally slow and steady, on the order of 1-100mma^-1 (Boulton and others, 1995). Basal melt tends to be distributed over wide areas of the bed, such as in patches on the upstream sides of bumps (Kamb, 1987), in areas of high strain such as along the edges of ice streams (Beem and others, 2010) or in areas of high subglacial geothermal flux (Vogel and Tulaczyk, 2006). In contrast, surface meltwater production rates can be quite high, on the order of 1000-10 000 mm a^-1 (Boulton and others, 1995). Supraglacial streams collect large volumes of this meltwater from supraglacial drainage catchments (McGrath and others, 2011) and rapidly deliver it to discrete points at glacier beds via moulins (Gulley and others, 2009a;Catania and Neumann, 2010).

The second key factor that has received limited attention is that spatio-temporal heterogeneity in the hydraulic conductivity of the glacier bed (taken here to include both the ice-bed interface and any underlying aquifer or subglacial till) will exert controls on patterns of water flow that are not captured by Eqn (1). For an equivalent hydraulic gradient, water fluxes increase as hydraulic conductivity increases. Conduits in soluble rocks, including limestone, salt, quartzite and ice, grow most rapidly along flow paths with the greatest solvent flux (Rothlisberger, 1972;Ewers, 1982;Groves and Howard, 1994;Gabrovsek and Drey-brodt, 2001;Rajaram and others, 2009), and the most conductive regions of a hydrological system that are connected to the greatest volume of recharge will grow the fastest (Rothlisberger, 1972;Ewers, 1982;Groves and Howard, 1994;Dreybrodt, 1996;Gabrovsek and Dreybrodt, 2001;Szymczak and Ladd, 2011).

Prediction of zones of the greatest subglacial water flux, and hence conduit locations, using theoretical hydraulic potential (Eqn (1)) implicitly assumes that subglacial recharge is spatially uniform and that subglacial hydraulic conductivities are homogeneous and isotropic. Interestingly, similar assumptions were made in models of conduit formation in limestone karst aquifers that were popular at the time that Shreve's theory was published (cf. Rhoades and Sinacori, 1941; Thrailkill, 1968). Later work, however, has shown that these simple models do not represent flow conditions in limestone aquifers. Heterogeneity and anisotropy in both aquifer recharge and hydraulic conductivity were found to result in flow paths that were structurally controlled and differed greatly from flow paths that were inferred from purely topographic analysis (Ford and Ewers, 1978;White, 1988;Palmer, 1991, 2007;Ford and Williams, 2007). Conduit networks were found to grow from a preexisting structurally defined architecture that concentrates water flow into discrete regions that have no direct causal relationship with surface topography. As a result, conduit density is primarily determined by the number of locations on the surface where water sinks into the aquifer along high- permeability pathways, rather than by the total surface area of the drainage basin.

Because glaciers and karst aquifers have similar hydrogeological settings (Gulley and others, 2009a), we use karst aquifer hydrology as an analogue system to explore how heterogeneities in meltwater flow at, and recharge to, glacier beds might result in conduit systems with different orientations and hydraulic gradients than predicted by hydraulic potential theory (cf. Shreve, 1972).

Karst aquifers and glaciers: similar hydrogeological settings

Conduits in limestone have two potential sources of water: (1) discrete recharge by surface runoff and (2) diffuse recharge from the aquifer matrix (Shuster and White,1971. In conduits that are recharged only by diffuse infiltration, hydraulic head is always lower in conduits than in the aquifer matrix. This lower head results because both the conduit and the aquifer are recharged at the same rate but the greater hydraulic conductivity of the conduit allows the conduit to discharge water more efficiently than the adjacent aquifer matrix. Because of continuity, the conduit must exist at lower pressure than the aquifer matrix (White, 1999). Diffusely recharged conduits are therefore entirely dendritic because water always flows from the aquifer to the conduit (Palmer, 1991).

While conduits were once thought to be solely low- pressure drains of a higher-pressure aquifer matrix, it was later discovered that conduits connected to surface recharge points such as swallets (analogous to moulins in glaciers) could sometimes function as higher-pressure sources of water to aquifers (White, 1999). When conduits are recharged more quickly than the adjacent aquifer matrix during storm events, conduit heads increase faster than groundwater heads, driving water from the conduit into the most highly conductive regions of the matrix, such as fissures and bedding planes, dissolving additional bypass routes to accommodate increased drainage (Ewers, 1982; Palmer, 1991, 2001;Moore and others, 2010). These bypass routes result in conduits with anastomosing plan forms, although this anastomosing pattern is typically superimposed upon an overall dendritic drainage pattern (Palmer, 1991, 2001). As recharge decreases following the storm event, water level drops more rapidly in the conduit than the surrounding aquifer matrix so that hydraulic gradients reverse and water drains back into the conduit (Ewers, 1982;Palmer, 1991, 2001;Moore and others, 2010).

The mode of conduit recharge during periods of conduit enlargement is therefore one of the strongest controls on conduit planform morphologies, and links between conduit morphology and hydrology are well established (White, 1988;Palmer, 1991, 2007;Ford and Williams, 2007). Because morphology depends largely on hydraulic rather than dissolution processes, similar relationships between morphology and hydrology have been found in conduits of widely varying lithologies (e.g. limestone, dolomite, salt, quartzite) (Klimchouk, 2000). Such relationships allow the hydrology of relict conduits to be inferred from morphology (Palmer, 1991;Klimchouk, 2000;Ford and Williams, 2007; Filipponi and others, 2009) and conduit morphologies have proven to be valuable resources for reconstructing the paleohydrology of landscapes (Granger and others, 2001; Stock and others, 2004;Frumkin and Fischhendler, 2005; Haeuselmann and others, 2007).

Unlike in recent developments in karst hydrology, however, the question of how meltwater delivery to glacier beds affects subglacial hydrology has received little attention. Spatial variability in the size and interconnectedness of linked cavities (Sharp and others, 1989) or spatial variability in the hydraulic conductivity of subglacial till (Punkari, 1997 introduces heterogeneity in the hydraulic conductivity of the subglacial hydrological system that violates assumptions of homogeneous flow conditions. Spatio- temporal variability in surface melt rates (Hock, 2003) violates assumptions of uniform recharge rate. Finally, moulins collect water from large surface catchments and deliver it to discrete points at glacier beds, introducing spatial variability to subglacial recharge (Gulley and others, 2009a;Catania and Neumann, 2010). How point recharge by moulins might interact with subglacial drainage systems that have a heterogeneous distribution of hydraulic conductivity, however, has not been evaluated. The similar hydrogeological settings between glaciers and limestone karst aquifers suggest that similar hydraulic processes might occur in both systems.

In this study, we use detailed geomorphic maps of a subglacial conduit beneath Hansbreen, a polythermal glacier in Svalbard, Norway, to infer hydraulic relationships between conduits and distributed systems. The subglacial conduit in our study extended from the base of a moulin that formed in the wall of an ice-marginal lake basin at the confluence of Fuglebreen and Hansbreen (Fig. 1). We compare the conduit map with glacier surface and bed topography, hydraulic potential calculated from Eqn (1), and theoretically predicted conduit morphologies (Shreve, 1971). These data suggest an alternative conceptual model for understanding the distribution and hydrology of subglacial conduits.

Fig. 1. (a) Hansbreen, a polythermal glacier near the Polish Polar Station in Svalbard. Our study area is indicated by a black box located at the glacier confluence indicated by the arrow. Elevation contours are in 50m intervals. (b) The supraglacial lake basin containing a relict hydrofracture; a person has been circled to indicate scale. The black arrow connecting the photo of the ice-marginal lake to its location on the map is also approximately the same as the view direction of the photograph. Three vertically oriented arrows point to a thermoerosion notch formed by a supraglacial lake in 2008. White stars denote sink-points for two ice-marginal streams.

Methods

Conduits were mapped in October 2008 and again in September 2009 using speleological techniques modified for glaciers (Gulley, 2009). Ice surface and bed topography data were collected from 18 radar profiles, totaling more than 5.5 km of survey lines, on 18 September 2009. We used a Mala Geosciences ground-penetrating radar (GPR) consisting of rough terrain antenna (RTA), control unit and data acquisition platform. The RTA consisted of an in-line flexible unit that was 9.25 m long and contained both transmitter and receiver. The transmitter center frequency was 50 MHz and the transmitter-receiver distance was fixed at 4 m.

Surveys were conducted in a 0.177 km² area of the glacier surface directly above or adjacent to the conduit (Figs 1 and 2). GPR signals were collected every 0.5 s and were simultaneously positioned by a dual frequency GPS receiver working in differential kinematic mode (DGPS). The processing of GPR data included adjustment of time-zero, correction of signal saturation and correction of amplitude and migration. The conversion to depth was obtained from estimated radio-wave velocity (RWV). RWV was calculated based on travel time to diffracting objects within ice, such as water- or air-filled englacial conduits or the glacier bed (Moore and others, 1999;Benjumea and others, 2003). The shape of hyperboles from radio-wave diffraction was used for RWV calculation. The mean RWV estimated based on all hyperbolas was 16.5 cm ns^-1; however, the average RWV within the cold ice layer in the vertical range 0-300 ns was 16.8 cm ns^-1. The differences in RWV are related to higher water content in the temperate ice layer than in overlaying cold dry ice. As the major fraction of ice in the area of the subglacial conduit was recognized as cold ice we adopt RWV 16.8 cm ns^-1 for further analysis (Fig. 2). Vertical accuracy of GPR sounding is estimated to be ±0.8 m, corresponding to a quarter of the wavelength, and horizontal errors are approximately ±0.1 m. The relative thickness of ice that was below the pressure-melting point versus the deeper temperate ice was determined from the depth that radargrams transitioned from transparent (cold ice) to scatter-rich (temperate ice) (cf. Gusmeroli and others, 2012).

Fig. 2. Migrated GPR survey line showing the glacier bed, conduit and polythermal structure (survey line and direction indicated by a white star and arrow in Fig. 6b). Additional survey lines depicted in Figure 6b were used to construct glacier surface and bed DEMs.

Data from DGPS/GPR tracks (Fig. 2) were used to create surface and bedrock DEMs with a 20 m grid. A mask was manually selected with an approximately 20 m margin around outer profiles. Surface and bed DEMs were used to calculate and contour subglacial hydraulic potential using Eqn (1). Conduits that are controlled by ice pressure should be normal to contours of equivalent potential, whereas conduits that form under atmospheric pressure should follow the bed slope (Sharp and others, 1993).

Results

Conduit description

In both 2008 and 2009, we reached the bed at a depth of 43 m below the ice surface by rappelling into a fracture that was located in the wall of an ice-marginal lake (Fig. 1). The fracture was vertical, oriented perpendicular to the ice flow direction and terminated at the glacier bed in ice that was below the pressure-melting point. Since new fracturing events were not observed during the study period, this conduit was interpreted to be a relict feature, most likely a hydrofracture. Two small supraglacial streams sank on either side of the hydrofracture entrance (Fig. 1b), forming cut-and- closure conduits that migrated upstream from the hydrofracture as nickpoints (cf. Gulley and others, 2009b). All portions of the conduit downstream from the hydrofracture displayed intact ceilings and appeared to have formed entirely subglacially.

In both 2008 and 2009 we mapped ~0.4 km of conduit to a maximum depth of ~95 m below the ice surface (Fig. 3). The conduit system was anastomotic and looked similar in both years (Fig. 3), although slightly more conduit was accessible in 2009 and conduit cross sections were slightly larger (Fig. 3). In both years, all conduits were followed until they either became too constricted by creep to continue or they became plugged with till from debris slides that commonly occurred in the steeply dipping conduit floors. Beginning at the base of the hydrofracture (Fig. 4a), the floors of the conduits consisted primarily of large boulders and were largely depleted of fine sediments (Fig. 4a-d and f). Where till was visible at the ice-bed interface, it was unsorted and had a dense matrix of fine sediments (Fig. 4b). Conduit cross sections varied from half-tubes incised upwards in the ice (Fig. 4b, c and f) to low and wide tunnels (Fig. 4d).

Fig. 3. Subglacial conduit maps from (a) October 2008 and (b) September–October 2009. Conduit cross sections are drawn 2_ plan view scale. The locations of photographs depicted on Figure 4 are marked with a ‘P’, followed by a lower-case letter indicating the corresponding panel. Unless otherwise indicated, the direction of view in the photographs is in the same direction as the cross section closest to the photograph. (c) Profile view of conduit mapped in October 2008.

Fig. 4. Photographs of subglacial conduit morphologies. Person shown for scale in all photographs. Locations of photographs refer to locations indicated in Figure 3. (a) Looking up the hydrofracture entrance. (b) Large well-sorted boulders on the floor of the conduit contrast with less sorted material in the till wall, at left. (c) R-channel morphology with highly scalloped conduit walls. (d) Conduit segment with low wide cross section and a boulder-covered floor. (e) Small ‘room’ formed by a landslide. Till in this area of the conduit failed, slipping into the conduit in the lower left-hand portion of the photograph and blocked additional exploration. (f) Classic R-channel morphology. While conduit cross sections may be similar to R-channels, the hydraulics of this system are different from that envisioned for R-channels as conduits are anastomotic not dendritic.

The main subglacial conduit extended from the base of the hydrofracture (Figs 1b and 4a) and plunged down the bed at an angle of ~35° (Figs 3 and 5), eventually splitting into several downstream continuations. We did not observe any upstream subglacial conduit extensions that were not associated with the hydrofracture. Similar to other conduits on Hansbreen, we crossed the cold-temperate boundary at an ice thickness of ~60m (Fig. 2). Cold-temperate transitions are readily identifiable in caves because the walls of conduits in temperate ice are wet to the touch and water is slowly dripping from the ceiling, whereas wet clothing freezes to the walls of conduits in cold ice. Despite the presence of temperate ice, no additional englacial conduits of any size were observed to feed into the subglacial conduit network.

Fig. 5. A 3-D view of the subglacial conduit showing spatial relationships between the subglacial conduit, ice and bed topography. The direction of view is shown in the small inset panel at the lower left-hand side. Survey data from September–October 2009. Universal Transverse Mercator (UTM) grid in meters. Conduit segment identifiers (i.e. (a), (b), (c)) correspond to the same conduit segment identifiers in Figure 6b.

Conduit relationship with ice thickness

In plan view, conduits extending northwest to southeast generally follow the trend of the decreasing ice surface elevation (Fig. 6a) and conduits trending to the northeast are generally oriented perpendicular to the ice slope and followed a 'gulley' in the bedrock topography (Fig. 6b). Some conduit segments cross ice surface contours at nearly right angles (Fig. 6a), continuing parallel to bed topography (Fig. 6b) and contours of ice thickness (Fig. 6c). These conduit segments followed the overall trend of hydraulic potential, but crossed contours at oblique rather than right angles (Fig. 6d). Other conduit segments were nearly parallel to contours of hydraulic potential (Fig. 6d) and generally followed the bed slope (Fig. 6b). East-trending conduit segments were generally oriented perpendicular to ice surface topography (Fig. 6a), closely following the strike of the cold-temperate ice boundary (Fig. 6e).

Fig. 6. Subglacial conduit map overlaid on: (a) ice topography; (b) bed topography (equivalent to hydraulic potential contours assuming subglacial water pressure was at atmospheric pressure); (c) ice thickness; (d) theoretical hydraulic potential (assuming water pressure and ice pressure were equal); and (e) percent cold ice thickness. In (b) dotted lines indicate GPR lines in the vicinity of the conduit. Some GPR lines extended outside the area shown in this figure and are not depicted. The figure area shown in each panel is equivalent to the area covered by the black box in Figure 1a. UTM data shown in meters. The star in (b) indicates the GPR line depicted in Figure 2 and the arrow indicates the direction of survey.

Discussion

Our conduit maps deviate in several important respects from predictions based on steady-state theories of glacier hydrology (Rothlisberger, 1972;Shreve, 1972). Our conduit began at a glacier confluence, whereas hydraulic potential theory does not predict confluences to have drainage (Pattyn and others, 2009). Mapped conduits were locally anastomotic (Fig. 3), not dendritic (although the overall architecture of the subglacial drainage system is very likely to be dendritic as there are many more moulins than there are proglacial discharge points). Finally, conduits extended parallel to, or crossed, contours of equivalent hydraulic potential at oblique angles (Fig. 6b and d). We explore the implications of each of these findings below within the context of published studies of glacier hydrology.

Conduit entrance location

If conduit network geometries and locations are controlled entirely by gradients in hydraulic potential (Shreve, 1972), conduits should not form at glacier confluences because confluences lack upstream ice that would control hydraulic potential (Pattyn and others, 2009). The subglacial conduit mapped in this study, however, began at the base of a moulin at a confluence (Fig. 1). It is important to note that the conduit began at the moulin as opposed to the moulin having intersected a conduit with upstream and downstream continuations.

Finding a major subglacial conduit in a glacier confluence suggests that the locations of recharge features, such as moulins, may be more important than theoretical hydraulic potential in controlling the upstream locations of subglacial conduits. Approximations of subglacial conduit networks using Eqn (1) do not consider the locations of moulins, and the locations of subglacial conduits are predicted solely on the basis of ice-thickness and elevation data. Moulins that fall close to inferred drainage axes are considered to be connected to the subglacial drainage system, whereas moulins that lie at some distance from these inferred drainage axes are considered to be unconnected (Willis and others, 2009). In many cases, the majority of moulins lie far from inferred major drainage axes (Palli and others, 2003) and it is not clear how these recharge features should interact with theoretically constructed subglacial drainage systems.

If, as our cave map suggests, subglacial conduits begin at moulins, calculation of hydraulic potential from ice topography and thickness alone does not take advantage of one of the most critical and easily obtainable pieces of information for the reconstruction of subglacial drainage systems: the locations of subglacial recharge. While hydraulic potential predicts confluences to be free of major subglacial drainage systems, confluences commonly have moulins that connect to glacier beds because they combine the glacier stress conditions and large meltwater supplies necessary to support hydrofracturing (Benn and others, 2009;Gulley and others, 2009a). The moulin in our study on Hansbreen has thus formed a recharge site to the glacier bed, as well as a subglacial conduit, in an area of the bed not predicted by hydraulic potential theory to have conduits.

Conduit relationships with hydraulic potential

Our maps of subglacial conduits only cover a small area of the total glacier bed (Fig. 1) and therefore reflect local rather than regional controls on subglacial water flow. Despite this local bias, several of our observations of relationships between subglacial conduits, glaciological data and reconstructions of hydraulic potential require discussion.

If gradients in theoretical hydraulic potential solely controlled water flow in our study, conduits should have followed ice surface topography and crossed contours of hydraulic potential at right angles (Fig. 6a and d). If bed topography solely controlled water flow in our study, then conduits should have crossed bed topography contours at right angles (Fig. 6b).

Comparisons of conduit maps with data derived from surface and bed surveys show that some conduit segments crossed ice surface contours at nearly right angles (Fig. 6a) and paralleled bed topography (Fig. 6b) and ice thickness (Fig. 6c). Some conduit segments extended in the general direction of hydraulic potential, but crossed hydraulic potential contours at oblique angles (Fig. 6d). Some of these segments followed the general direction of the bed slope (Fig. 6b). Conduits did not appear to be guided strictly by hydraulic potential (Fig. 6d) because conduit segments followed both the bed slope (Fig. 6b) and ice surface (Fig. 6a).

Individual segments of the conduit network roughly coincide with theoretical hydraulic potential theory. For instance, conduit segments that trended to the southeast generally followed gradients in theoretical hydraulic potential (Fig. 6d) and therefore could have formed when water pressure equaled the ice overburden pressure. The position of a thermal erosion notch in the ice-marginal lake basin is consistent with a conduit formation pressure near ice overburden pressure (Fig. 1b). Conduits that generally follow bedrock topography could have formed later when water pressure decreased to atmospheric pressure. Several observations suggest, however, that heterogeneity in flow paths at the bed may be important in controlling which areas of the bed enlarge to form major conduits. First, southeast-trending conduits crossed equipotential surfaces at oblique rather than right angles (Fig. 6d). While this oblique intersection could result from the fine scale at which we are examining spatial relationships, it is worth noting, particularly given that the subglacial conduit crossed contours of ice topography at angles that were much closer to 90°. Second, the northeast-trending conduit, segment (a), could not have formed at atmospheric pressure (Figs 5 and 6b) and therefore been controlled solely by flow down the glacier bed slope. The only source of recharge for segment (a) would have been when segment (b) was transmitting water, which would have required segment (c) to be completely water-filled (Fig. 5). Because segments (a) and (c) plunged down the bed at a ~35° angle, water pressure must have increased down the conduit during formation, resulting in a pressurized hydraulic system.

Conduit planform

Subglacial conduits are widely cited to be dendritic and to exist at lower pressure than adjacent distributed systems because they have been conceptualized to be recharged by a homogeneous englacial drainage system (Shreve, 1972). In contrast to these predictions, the conduit in our study was, at least locally, anastomotic. Anastomosing conduits are thought to form under ice streams with low surface gradients where uniformly distributed basal melt generates high water pressures at deformable glacier beds (Walder and Fowler, 1994). These conditions do not occur beneath Hansbreen, however, and thus cannot account for anastomosing in our conduit system.

Studies of conduit formation in bench-top models have shown that anastomosing conduits can form in planar discontinuities, such as bedding planes or ice-bed interfaces, when water is injected into the discontinuity more quickly than it can be evacuated (Ewers, 1982;Catania and Paola, 2001). Both braiding intensity and variability in flow- path direction increase as the rate of recharge to the system increases (Catania and Paola, 2001). This process is directly analogous to how anastomotic caves form in limestone aquifers that are recharged from allogenic catchments (Ewers, 1982;Palmer, 1991, 2001). We hypothesize these same hydraulic processes are responsible for the anastomosing conduit planform beneath Hansbreen.

We propose the following scenario to explain the formation of the anastomosing subglacial conduit in our study. At the beginning of the melt season, the supraglacial drainage basin delivered meltwater to the moulin at a rate that exceeded the hydraulic capacity of the subglacial conduit. As a result, meltwater backed up in the subglacial and englacial drainage systems, ulttimately spilling over into the ice-marginal lake basin and allowing it to fill with water. Water levels in the lake were stable for a sufficiently long period of time to melt a thermoerosional notch (Fig. 1b). Water that could not be accommodated by the primary drainage pathway was forced into the most highly conductive regions of the ice-bed interface, which were preferentially enlarged to form the anastomosing conduits that we mapped.

Reuse of the same preferential subglacial flow paths between years could have allowed the conduit in our study to form in the same location both years. Potential preferential flow paths include constricted subglacial conduits that are kept open by water backing up in englacial conduits (Benn and others, 2009;Gulley and others, 2009a;Catania and Neumann, 2010) or zones of subglacial till with higher subglacial hydraulic conductivity (Harbor and others, 1997; Punkari, 1997). As a result, feedbacks between subglacial conduit flow paths, winnowing of fine-grained material and increased hydraulic conductivity of subglacial till could allow conduits to follow similar flow paths each year, even if they were closed entirely by creep during winter.

Additional support for the existence and importance of preferential subglacial flow paths comes from several published dye-tracing studies. First, some dye traces injected into crevasses on Cascade Glacier, Washington, USA, were found to cross a subglacial drainage divide that was predicted by theoretical hydraulic potential (Fountain and Vaughn, 1995). Crossing drainage divides defined by theoretical potential is only possible if heterogeneity in the hydraulic conductivity of glacier beds is more important in determining the direction of water flow than gradients in ice thickness. Second, dye traces that are injected into moulins are widely reported to have residence times of hours (Fountain, 1993;Hock and Hooke, 1993;Kohler, 1995; Schuler and others, 2004;Werder and others, 2010), whereas dye traces that are injected into boreholes have residence times of days to weeks (Moeri and Leibundgut, 1986;Iken and Truffer, 1997;Hock and others, 1999). These studies indicate that there are strong heterogeneities in flow paths at glacier beds and that moulins may be connected to regions of higher subglacial hydraulic conductivity than other areas of the bed. If preferential flow paths were not important, dye traces conducted in moulins at the beginning of melt seasons should have breakthrough curves that bear some similarity to breakthrough curves of dye traces that are conducted from boreholes. However, dye traces that are conducted even in the earliest parts of melt seasons have transit times that are measured in hours instead of the days or weeks required for borehole traces (cf. Nienow and others, 1998).

Finally, an interesting dye trace study was conducted on Findelengletscher, Switzerland, in an area of the glacier where ice was ~165 m thick (Moeri and Leibundgut, 1986; Iken and Truffer, 1997). Two different dyes were injected into two separate boreholes during the winter when surface melting was absent. One dye was injected into a borehole that was drilled to the bed through intact glacier ice. The second dye was injected into a borehole that was drilled to the bed through a moulin. Both traces were flushed with 1.5 m³ of water. Dye injected into the borehole drilled in intact ice began emerging from the glacier 86 hours after injection, and two separate concentration maxima were recorded 137 hours and 4 weeks after injection. In contrast, dye injected into the moulin began emerging at the glacier front 3 hours after injection, and peak concentration was reached 7 hours after injection (Moeri and Leibundgut, 1986). We argue that this last study provides strong support that moulins are connected to preferential flow paths at glacier beds and that these preferential pathways can persist through winter.

In addition to providing a possible explanation for the formation of an anastomosing conduit that was located in an area of the bed that had been predicted by hydraulic potential theory to lack conduits, preferential flow paths could also explain why many ice-dammed lakes drain before flotation values are reached (Tweed and Russell, 1999;Bjornsson, 2010). According to hydraulic potential theory, hydraulic gradients should direct water flow towards lakes where water pressure is lower than the ice overburden pressure, thereby prohibiting lake drainage. Preferential flow paths, however, could allow conduits to form as soon as head in the lake basin increased sufficiently to drive water along preferential flow paths to initiate conduit formation. While effective pressures would still be positive if water pressure was below ice overburden pressure, conduit enlargement rates could still outpace rates of creep closure due to the large meltwater supply in the lake and the reduced rates of creep closure associated with lower effective pressures. Such a formation mechanism would also allow conduits to form and to drain lakes before the water level increased sufficiently to float the ice dam to create accommodation space for conduit formation (Bjornsson, 2010).

We hypothesize that conduits begin at lakes, moulins or other similar recharge points and then form along preferential flow paths as soon as rates of enlargement exceed rates of creep closure. In this scenario, the effects of subglacial water pressure are secondary, acting to constrict conduit diameters along their chosen flow paths rather than forcing a reconfiguration of conduit networks based on the new effective pressure (cf. Rippin and others, 2003). The ability of conduit systems with fixed flow paths to persist in spite of seasonally and diurnally fluctuating effective pressures has been widely investigated with pipe-flow models (Spring and Hutter, 1981;Arnold and others, 1998; Clarke, 2003). Such stability highlights the capacity of conduit network configurations to persist in fixed configurations under positive effective pressures because conduit enlargement rates can offset creep closure.

How recharge-discharge relationships may affect subglacial water pressure

Hydraulic relationships between preferential flow paths (including conduits) and adjacent areas of glacier beds could be determined by the mechanism of recharge, similar to karst aquifers. Where glaciers are recharged entirely by diffuse basal melt and moulins are not present, such as beneath Antarctic ice streams, conduits would be expected to exist at lower water pressure than distributed systems (Fig. 7a and b). Lower conduit pressure would occur because both the conduit and the distributed system are recharged at nearly the same rate. Because the conduit is more hydraulically conductive than the distributed system, however, water discharges faster and, because of continuity, must have a lower pressure.

Fig. 7. Different recharge-discharge relationships should affect the distribution of subglacial hydraulic head differently. Top panels are plan views of a subglacial hydrological system and black dots represent boreholes. Bottom panels show the distribution of hydraulic head (in cross-sectional view) of each type of subglacial drainage represented in the upper panels. The spatial variability of hydraulic head is shown for low rates of meltwater delivery (low recharge rate) as well high rates of meltwater delivery (high recharge rate). (a) Spatially uniform recharge to a distributed subglacial drainage system with a homogeneous hydraulic conductivity. Hydraulic head is uniformly high at high recharge rates and uniformly low at low recharge rates. (b) Spatially uniform recharge to distributed system that is drained by a conduit. Hydraulic head is higher in the distributed system and lower in the conduit when recharge rates are high as well as when recharge rates are low. (c) Discrete recharge to a subglacial drainage system by a moulin. The highest heads occur in the conduit when recharge to the moulin is increasing, such as during the day. As recharge begins decreasing at night, hydraulic head drops faster in the conduit than in the distributed system.

Where conduits or other preferential flow paths (i.e. the most hydraulically conductive regions of linked cavity systems) are recharged by moulins that drain large supraglacial catchments, surface runoff should increase hydraulic head in the conduit faster than basal melting can increase head in the distributed system (Fig. 7c). As a result, recharge to moulins will create hydraulic gradients that drive water from the conduit into the distributed system (Figs 7c and 8a). Such reversals in gradient should occur even at the beginning of melt seasons when subglacial water pressure might be at the ice overburden pressure because of basal melting in a low-capacity drainage system. Exchange can occur because meltwater that backs up to the level of the ice surface in moulins would elevate hydraulic head in preferential flow paths. Because of the difference in the density of ice and water, maximum water pressures at the base of the moulin and in connected high hydraulic conductivity flow paths would be ~10% greater than the maximum possible water pressure in the distributed system. Hydraulic head within high conductivity flow paths or conduits could be higher than the ice thickness down the flow path if head losses were small and the surface topography was steeply dipping (Clarke, 2003). During the recession limb of hydrographs, hydraulic head decreases most rapidly in the most conductive regions first, causing hydraulic gradients to reverse and allowing surface water that was injected into the distributed system during the day to drain back into the conduit (Fig. 8b).

Fig. 8. (a) In conduits recharged by moulins (i.e. Fig. 7c), rapid delivery of large volumes of surface water to the upstream end of conduits causes head to increase faster in the conduit than in the distributed system during the day. This increase drives water from the conduit into the distributed system, ultimately forming anastomosing conduits. (b) As recharge rates to moulins decrease at night, conduit heads decrease more quickly than heads in the distributed system because the conduit has a greater hydraulic capacity. Head gradients reverse and water flows back into conduits.

While we argue for diurnally varying head gradients between conduits and adjacent areas of the bed on the basis of anastomosing conduits (Fig. 6), diurnal gradient reversals have been observed in studies of subglacial hydrology. A series of investigations that occurred over several years on Haut Glacier d'Arolla found water pressure in boreholes installed in a subglacial conduit increased faster and to higher values than water pressure in boreholes installed in adjacent areas of the bed during the day (Hubbard and others, 1995;Harbor and others, 1997;Mair and others, 2002). Conduit water pressure decreased faster and to lower values than water pressure in adjacent boreholes at night.

Additional evidence that diurnally varying head gradients between conduits and distributed systems are being driven by discrete recharge to preferential flow paths comes from comparisons of studies in which dye was injected into boreholes (Hock and others, 1999) with studies in which dye was injected directly into moulins (Nienow and others,1997. All dye injected into moulins discharged from the glacier within hours of injection, even at the beginning of the melt season, indicating that the dye flowed rapidly through interconnected zones of high hydraulic conductivity. Conversely, dye injected into boreholes discharged over a period of weeks and diurnal dye fluxes were inversely proportional to diurnal discharge. Little or no dye discharged from the glacier during the rising limb of proglacial hydrographs or at peak discharge, suggesting that the dye was trapped in the distributed system by head gradients that were oriented away from the conduits (Fig. 8a). Dye flux increased on the falling limb of proglacial hydrographs, reaching a maximum during or near proglacial hydrograph minima, reflecting a reversal of head gradients (Fig. 8b). Dye was therefore trapped in the distributed system during the day when moulins rapidly delivered meltwater to conduits and increased conduit heads relative to the distributed system (Fig. 8a). As melt rates and conduit heads decreased at night, hydraulic gradients were redirected towards conduits and small amounts of dye would be able to flow into the conduit (Fig. 8b). These diurnal changes in head are a plausible explanation for the observed pulses in dye discharge.

Exchange of water between conduits and distributed systems could also explain why ice sliding speeds do not always show strong correlations with subglacial water pressure that is measured in boreholes. With flow from conduits to distributed systems (Fig. 8a), head increases in the distributed system would be attenuated as a positive function of distance from the conduit and an inverse function of hydraulic conductivity (Hubbard and others, 1995;Harbor and others, 1997;Gordon and others, 1998). Attenuation and lag of water pressure with distance from conduits (Hubbard and others, 1995; Alley, 1996;Harbor and others, 1997; Gordon and others, 1998;Martin and Dean, 2001) could explain why pressures measured in boreholes do not always correlate with sliding (Harper and others, 2005) since boreholes may be located far from conduits (Fig. 8a). This same connection between moulins, conduits and distributed systems would explain correlations between moulin water level and sliding speed (Vieli and others, 2004). This link between conduits and distributed systems indicates that comparisons of water pressure in moulins and boreholes would be required to fully understand the spatial and temporal distribution of subglacial water pressures.

Theoretical hydraulic potential revisited

While we propose that hydraulic potential theory overlooks the importance of recharge location and heterogeneity in determining conduit network locations, hydraulic potential has predicted the locations of conduits in some (Harbor and others, 1997) but not all (Hantz and Lliboutry, 1983) cases.

Subglacial hydraulic potential theory might sometimes predict the locations of subglacial drainage systems, but perhaps it does so for the wrong reasons. We consider the most likely reason that hydraulic potential sometimes locates conduits is because glacier surface topography controls the locations of supraglacial drainage that ultimately recharges conduits at moulins. Troughs in surface topography, which are often controlled by subglacial topography, control the positions of supraglacial streams and lakes (Vidonish and others, 2010), which ultimately recharge moulins. Since moulins form where crevasses intersect supraglacial streams, both moulins and crevasses control the positions of the upstream ends of subglacial conduits. These features are essentially in a fixed reference frame relative to a moving glacier and would allow recharge to occur in the same places year after year. Drainage reconstructions using theoretical hydraulic potential may therefore actually reflect the influence of supraglacial recharge on subglacial hydrology. In support of this hypothesis, subglacial conduit networks reconstructed using hydraulic potential have been shown to be closely similar to supraglacial drainage systems (Sharp and others, 1993), and the subglacial conduit that was predicted using hydraulic potential and studied intensively on Haut Glacier d'Arolla was located just down-glacier of a field of moulins (Mair and others, 2003). Additionally, the highly channelized regions of the GrIS that experience summer increases in velocity are highly correlated with supraglacial drainage systems (Palmer and others, 2011).

Implications for modeling

Our findings have important implications for how the hydrological systems of glaciers with moulins are modeled. Because moulins are the primary source of surface melt- water recharge to glacier beds (Gulley and others, 2009a), the density and spatial distribution of subglacial conduits should be controlled by the number and locations of moulins rather than effective pressure (cf. Schoof, 2010; Hewitt, 2011) or the hydraulic potential (Shreve, 1972). As a result, conduits may not be concentrated near glacier center lines, as is frequently predicted from theoretical potential, but will begin at moulins and follow preferential flow paths at glacier beds. Although there may be no surface expression of these preferential flow paths, they should dominate subglacial flow similar to groundwater flow karst aquifers (Ford and Ewers, 1978;Palmer, 1991).

The subglacial conduit in our study did follow the overall trend indicated by theoretical hydraulic potential once the moulin had routed water to the bed (although it was located where hydraulic potential had predicted a lack of conduits). This finding suggests that models using hydraulic potential theory might be improved by treating each moulin as a recharge point and connecting moulins to known discharge points along lines of theoretical hydraulic potential, even if calculations of theoretical hydraulic potential may lack any physical meaning.

Realistic models of hydraulic interactions between conduits and distributed systems require inclusion of supraglacial recharge to glacier beds at discrete points and inclusion of subglacial flow paths that have a higher hydraulic conductivity than adjacent areas of the bed. Spatial heterogeneity in recharge and subglacial water flow results in hydraulic interactions that are not captured in models where water is uniformly and diffusely recharged to glacier beds that have initial homogeneous hydraulic conductivities.

Because moulins connect the glacier surface to the bed, the maximum attainable head will be determined by the elevation of the moulin, and excess recharge will overflow onto the glacier surface. Such overflow is not possible in models that inject water at prescribed rates to the base of ice sheets (Pimentel and Flowers, 2011), which is why these models generate unrealistically high subglacial water pressures.

While it is widely cited that conduits decrease, and distributed systems increase, subglacial water pressure and rates of sliding (Bell, 2008), moulins recharging conduits provide a mechanism whereby high basal water pressure and increases in glacier ice velocity can occur in conduit drainage systems.

Summary and Conclusions

The characteristics of subglacial drainage systems as determined from our direct exploration and the results of previously published studies are inconsistent with subglacial hydrological models based on hydraulic potential theory (cf. Rothlisberger, 1972;Shreve, 1972). Fundamental assumptions used to calculate subglacial hydraulic potential (i.e. that recharge to the bed is spatially and temporally homogeneous and that the hydraulic conductivity of the ice-bed interface is homogeneous) are inconsistent with observations. Discrete recharge by moulins and heterogeneous flow paths at glacier beds result in flow paths that may be different from flow paths predicted by hydraulic potential theory and interactions between conduits and distributed systems that are much more complex.

Subglacial conduits begin at moulins and form along zones of the highest hydraulic conductivity (preferential flow paths) that link moulins to discharge points at ice margins. The distribution of subglacial conduits is thus controlled to first order by the locations of moulins and the locations of preferential flow paths at glacier beds. These preferential flow paths may occur as zones of highly sorted till or in efficient pathways through systems of linked cavities. Because of these controls, preferential flow paths may not coincide with drainage axes defined by theoretical hydraulic potential (cf. Shreve, 1972). Although ice pressure does not directly control the locations of subglacial conduits, ice pressure does exert an important control on subglacial conduit hydrology. The effects of ice pressure are secondary, however, acting to constrict conduit diameters along their chosen flow paths during periods of increased effective pressure.

Discrete recharge to preferential flow paths (including conduits) by moulins establishes hydraulic conditions that can explain many glacier hydrologic phenomena. Variations in recharge to moulins during diurnal heating and freezing cycles will drive observed diurnal cycles of exchange between preferential flow paths and distributed systems. On the rising limb of moulin recharge hydrographs, conduit heads increase faster than heads in the distributed system, driving water from the conduit into the distributed system, and will form anastomosing conduits. Recharge to the distributed system therefore originates from conduits. As a result, the highest water pressure and sliding speeds should be above conduits instead of distributed systems. Lags in the response of the distributed system to recharge pulses originating in conduits may explain why records of subglacial water pressure measured in boreholes sometimes correlate poorly with records of ice velocity. On the falling limb of the moulin recharge hydrograph, water that was recharged to the distributed system flows back into the conduit. Increased subglacial water pressure and sliding can occur whenever the recharge rate exceeds the hydraulic capacity of the system and the configuration of the drainage system is likely to be largely irrelevant (cf. Bartholomaus and others, 2008).