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 (Reference BellBell, 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 (Reference RothlisbergerRothlisberger, 1972;Reference ShreveShreve, 1972) and these studies have become fundamental components of glacier hydrological theory (Reference Fountain and WalderFountain and Walder, 1998; Reference BellBell, 2008;Reference Cuffey and PatersonCuffey 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 (Reference ShreveShreve, 1972;Reference Flowers and ClarkeFlowers and Clarke, 2002;Reference Flowers, Bjornsson and PalssonFlowers and others, 2003, Reference Flowers, Bjornsson, Palsson and Clarke2004; Reference HewittHewitt, 2011), draw down water pressure in adjacent distributed systems and ultimately decrease glacier sliding speeds (Reference Mair, Nienow, Sharp, Wohlleben and WillisMair and others, 2002;Reference BellBell, 2008;Reference SchoofSchoof, 2010). Contrasting these theoretical predictions are field observations that water pressure can increase faster and to higher magnitudes in conduits than in distributed systems (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995; Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997;Reference Gordon, Sharp, Hubbard, Smart, Ketterling and WillisGordon and others, 1998;Reference Mair, Willis, Fischer, Hubbard, Nienow and HubbardMair and others, 2003), causing the highest ice velocities to occur above conduits rather than above adjacent distributed systems (Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997;Reference Mair, Willis, Fischer, Hubbard, Nienow and HubbardMair 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 (Reference Palmer, Shepherd, Nienow and JoughinPalmer 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 (Reference AndersonAnderson and others, 2003; Reference Kessler and AndersonKessler and Anderson, 2004; Reference WalderWalder and others, 2006;Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2008, Reference Bartholomaus, Anderson and Anderson2011).
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 (Reference Munson, Young and OkiishiMunson 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 (Reference Freeze and CherryFreeze 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 Reference ShreveShreve (1972).
Reference ShreveShreve (1972) proposed that the hydraulic potential ($) of glacier drainage systems could be described at any point (z) in a glacier by the equation:
where Φ0 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 (Reference SharpSharp 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 (Reference SharpSharp and others, 1993;Reference RichardsRichards and others, 1996; Reference Palli, Moore, Jania, Kolondra and GlowackiPalli and others, 2003;Reference RippinRippin and others, 2003;Reference Willis, Lawson, Owens, Jacobel and AutridgeWillis and others, 2009), the GrIS (Reference AhlstrømAhlstrøm and others, 2002, Reference Ahlstrøm, Mohr, Reeh, Christensen and Hooke2005; Reference MottramMottram and others, 2009) and the Antarctic ice sheet (Reference Fricker and ScambosFricker 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 (Reference GulleyGulley 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 (Reference SharpSharp 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 (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995;Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor 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 (Reference Hantz and LliboutryHantz 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 (Reference GulleyGulley 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 (Reference Boulton, Caban and Van GijsselBoulton 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 (Reference KambKamb, 1987), in areas of high strain such as along the edges of ice streams (Reference Beem, Jezek and Van der VeenBeem and others, 2010) or in areas of high subglacial geothermal flux (Reference Vogel and TulaczykVogel and Tulaczyk, 2006). In contrast, surface meltwater production rates can be quite high, on the order of 1000-10 000 mm a-1 (Reference Boulton, Caban and Van GijsselBoulton and others, 1995). Supraglacial streams collect large volumes of this meltwater from supraglacial drainage catchments (Reference McGrath, Colgan, Steffen, Lauffenberger and BalogMcGrath and others, 2011) and rapidly deliver it to discrete points at glacier beds via moulins (Reference GulleyGulley and others, 2009a;Reference Catania and NeumannCatania 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 (Reference RothlisbergerRothlisberger, 1972;Reference EwersEwers, 1982;Reference Groves and HowardGroves and Howard, 1994;Reference Gabrovsek and DreybrodtGabrovsek and Drey-brodt, 2001;Reference Rajaram, Cheung and ChaudhuriRajaram 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 (Reference RothlisbergerRothlisberger, 1972;Reference EwersEwers, 1982;Reference Groves and HowardGroves and Howard, 1994;Reference DreybrodtDreybrodt, 1996;Reference Gabrovsek and DreybrodtGabrovsek and Dreybrodt, 2001;Reference Szymczak and LaddSzymczak 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; Reference ThrailkillThrailkill, 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 (Reference Ford and EwersFord and Ewers, 1978;Reference WhiteWhite, 1988;Reference PalmerPalmer, 1991, Reference Palmer2007;Reference Ford and WilliamsFord 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 (Reference GulleyGulley 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. Reference ShreveShreve, 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 (Reference White, Palmer, Palmer and SasowskyWhite, 1999). Diffusely recharged conduits are therefore entirely dendritic because water always flows from the aquifer to the conduit (Reference PalmerPalmer, 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 (Reference White, Palmer, Palmer and SasowskyWhite, 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 (Reference EwersEwers, 1982; Reference PalmerPalmer, 1991, Reference Palmer2001;Reference Moore, Martin, Screaton and NeuhoffMoore 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 (Reference PalmerPalmer, 1991, Reference Palmer2001). 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 (Reference EwersEwers, 1982;Reference PalmerPalmer, 1991, Reference Palmer2001;Reference Moore, Martin, Screaton and NeuhoffMoore 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 (Reference WhiteWhite, 1988;Reference PalmerPalmer, 1991, Reference Palmer2007;Reference Ford and WilliamsFord 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) (Reference KlimchoukKlimchouk, 2000). Such relationships allow the hydrology of relict conduits to be inferred from morphology (Reference PalmerPalmer, 1991;Reference KlimchoukKlimchouk, 2000;Reference Ford and WilliamsFord and Williams, 2007; Reference Filipponi, Jeannin and TacherFilipponi and others, 2009) and conduit morphologies have proven to be valuable resources for reconstructing the paleohydrology of landscapes (Reference Granger, Fabel and PalmerGranger and others, 2001; Reference Stock, Anderson and FinkelStock and others, 2004;Reference Frumkin and FischhendlerFrumkin and Fischhendler, 2005; Reference Haeuselmann, Granger, Jeannin and LauritzenHaeuselmann 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 (Reference Sharp, Gemmell and TisonSharp and others, 1989) or spatial variability in the hydraulic conductivity of subglacial till (Reference PunkariPunkari, 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 (Reference HockHock, 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 (Reference GulleyGulley and others, 2009a;Reference Catania and NeumannCatania 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.
Conduits were mapped in October 2008 and again in September 2009 using speleological techniques modified for glaciers (Reference GulleyGulley, 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 km2 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 (Reference MooreMoore and others, 1999;Reference Benjumea, Macheret, Navarro and TeixidoBenjumea 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. Reference Gusmeroli, Jansson, Pettersson and MurrayGusmeroli and others, 2012).
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 (Reference SharpSharp and others, 1993).
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. Reference GulleyGulley 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).
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.
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).
Our conduit maps deviate in several important respects from predictions based on steady-state theories of glacier hydrology (Reference RothlisbergerRothlisberger, 1972;Reference ShreveShreve, 1972). Our conduit began at a glacier confluence, whereas hydraulic potential theory does not predict confluences to have drainage (Reference Pattyn, Delcourt, Samyn, De Smedt and NolanPattyn 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 (Reference ShreveShreve, 1972), conduits should not form at glacier confluences because confluences lack upstream ice that would control hydraulic potential (Reference Pattyn, Delcourt, Samyn, De Smedt and NolanPattyn 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 (Reference Willis, Lawson, Owens, Jacobel and AutridgeWillis and others, 2009). In many cases, the majority of moulins lie far from inferred major drainage axes (Reference Palli, Moore, Jania, Kolondra and GlowackiPalli 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 (Reference Benn, Gulley, Luckman, Adamek and GlowackiBenn and others, 2009;Reference GulleyGulley 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.
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 (Reference ShreveShreve, 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 (Reference Walder and FowlerWalder 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 (Reference EwersEwers, 1982;Reference Catania and PaolaCatania and Paola, 2001). Both braiding intensity and variability in flow- path direction increase as the rate of recharge to the system increases (Reference Catania and PaolaCatania and Paola, 2001). This process is directly analogous to how anastomotic caves form in limestone aquifers that are recharged from allogenic catchments (Reference EwersEwers, 1982;Reference PalmerPalmer, 1991, Reference Palmer2001). 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 (Reference Benn, Gulley, Luckman, Adamek and GlowackiBenn and others, 2009;Reference GulleyGulley and others, 2009a;Reference Catania and NeumannCatania and Neumann, 2010) or zones of subglacial till with higher subglacial hydraulic conductivity (Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997; Reference PunkariPunkari, 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 (Reference Fountain and VaughnFountain 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 (Reference FountainFountain, 1993;Reference Hock and HookeHock and Hooke, 1993;Reference KohlerKohler, 1995; Reference Schuler, Fischer and GudmundssonSchuler and others, 2004;Reference Werder, Schuler and FunkWerder and others, 2010), whereas dye traces that are injected into boreholes have residence times of days to weeks (Reference Moeri and LeibundgutMoeri and Leibundgut, 1986;Reference Iken and TrufferIken and Truffer, 1997;Reference Hock, Iken and WanglerHock 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. Reference Nienow, Sharp and WillisNienow 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 (Reference Moeri and LeibundgutMoeri and Leibundgut, 1986; Reference Iken and TrufferIken 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 m3 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 (Reference Moeri and LeibundgutMoeri 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 (Reference Tweed and RussellTweed and Russell, 1999;Reference BjornssonBjornsson, 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 (Reference BjornssonBjornsson, 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. Reference RippinRippin 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 (Reference Spring and HutterSpring and Hutter, 1981;Reference Arnold, Richards, Willis and SharpArnold and others, 1998; Reference ClarkeClarke, 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.
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 (Reference ClarkeClarke, 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).
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 (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995;Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997;Reference Mair, Nienow, Sharp, Wohlleben and WillisMair 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 (Reference Hock, Iken and WanglerHock 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 (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995;Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997;Reference Gordon, Sharp, Hubbard, Smart, Ketterling and WillisGordon and others, 1998). Attenuation and lag of water pressure with distance from conduits (Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995; Reference AlleyAlley, 1996;Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997; Reference Gordon, Sharp, Hubbard, Smart, Ketterling and WillisGordon and others, 1998;Reference Martin and DeanMartin and Dean, 2001) could explain why pressures measured in boreholes do not always correlate with sliding (Reference Harper, Humphrey, Pfeffer, Fudge and O'NeelHarper 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 (Reference Vieli, Jania, Blatter and FunkVieli 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 (Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997) but not all (Reference Hantz and LliboutryHantz 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 (Reference Vidonish, Wooley, Cathles, Amundson, Darnell and MacAyealVidonish 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 (Reference SharpSharp 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 (Reference Mair, Willis, Fischer, Hubbard, Nienow and HubbardMair and others, 2003). Additionally, the highly channelized regions of the GrIS that experience summer increases in velocity are highly correlated with supraglacial drainage systems (Reference Palmer, Shepherd, Nienow and JoughinPalmer 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 (Reference GulleyGulley 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. Reference SchoofSchoof, 2010; Reference HewittHewitt, 2011) or the hydraulic potential (Reference ShreveShreve, 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 (Reference Ford and EwersFord and Ewers, 1978;Reference PalmerPalmer, 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 (Reference Pimentel and FlowersPimentel 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 (Reference BellBell, 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. Reference RothlisbergerRothlisberger, 1972;Reference ShreveShreve, 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. Reference ShreveShreve, 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. Reference Bartholomaus, Anderson and AndersonBartholomaus and others, 2008).
J. Gulley acknowledges funding from the US National Science Foundation in the form of a Graduate Research Fellowship, Nordic Research Opportunity, and an EAR Postdoctoral Fellowship (grant No. 0946767). Additional funding was provided to J. Gulley by the American Philosophical Society Lewis and Clarke Fund for Exploration and Field Research, the Evolving Earth Foundation, the University of Florida and the University Centre in Svalbard. M. Grabiec, P. Glowacki and J. Jania acknowledge funding from the ice2sea project of the European Union 7th Framework Programme (grant No. 226375). The Polish Polar Research Station, Hornsund, and the Polish Academy of Sciences Institute of Geophysics are thanked for logistical support. Alison Banwell, Doug Benn, Annelie Bergstrom and Artur Adamek are thanked for field assistance. Conversations with Doug Benn and Ian Willis helped improve the quality of the manuscript. We also thank Bob Anderson and an anonymous reviewer for reviews that improved the quality of the manuscript.