A numerical study of hydrologically driven glacier dynamics and subglacial flooding

Sam Pimentel, Gwenn E. Flowers


A hydrologically coupled flowband model of ‘higher order’ ice dynamics is used to explore perturbations in response to supraglacial water drainage and subglacial flooding. The subglacial drainage system includes interacting ‘fast’ and ‘slow’ drainage elements. The fast drainage system is assumed to be composed of ice-walled conduits and the slow system of a macroporous water sheet. Under high subglacial water pressures, flexure of the overlying ice is modelled using elastic beam theory. A regularized Coulomb friction law describes basal boundary conditions that enable hydrologically driven acceleration. We demonstrate the modelled interactions between hydrology and ice dynamics by means of three observationally inspired examples: (i) simulations of meltwater drainage at an Alpine-type glacier produce seasonal and diurnal variability, and exhibit drainage evolution characteristic of the so-called ‘spring transition’; (ii) horizontal and vertical diurnal accelerations are modelled in response to summer meltwater input at a Greenland-type outlet glacier; and (iii) short-lived perturbations to basal water pressure and ice-flow speed are modelled in response to the prescribed drainage of a supraglacial lake. Our model supports the suggestion that a channelized drainage system can form beneath the margins of the Greenland ice sheet, and may contribute to reducing the dynamic impact of floods derived from supraglacial lakes.

1. Introduction

In recent years, several authors have stressed the importance of including interactions between subglacial hydrology and ice dynamics in prognostic models (e.g. Bartholomaus et al. 2008; Bell 2008; Stearns et al. 2008). This is relevant to reducing uncertainties surrounding glacier and ice-sheet contributions to eustatic variations in sea level (Solomon et al. 2007). Glaciers and ice sheets respond to both forcing from meteorological inputs as well as dynamical changes at the glacier bed or margin, including those driven by the subglacial drainage system (e.g. Bartholomaus et al. 2008). Though many numerical ice-flow models have neglected aspects of hydrological transience on glacier motion, promising early approaches include those of Arnold & Sharp (2002) and Marshall et al. (2005), who calculate subglacial water pressures and use a water-pressure-dependent sliding law to calculate ice-sheet velocities. More mature modelling approaches are now beginning to emerge; in particular, the framework of Bueler & Brown (2009) holds much promise for large-scale ice-sheet deployments (Bueler et al. submitted). Bueler & Brown (2009) and Bueler et al. (submitted) take a hybrid approach in which the shallow-shelf approximation is used as a sliding law within a model that uses the shallow-ice approximation, and pore water pressure is dynamically related to basal melt rates. Model simulations show a good deal of correspondence with surface velocity observations over the Greenland ice sheet (Bueler et al. submitted).

For valley glacier settings, we have developed a hydrologically coupled higher order flowband (Blatter/Pattyn type) ice-dynamics model (Pimentel et al. submitted) in which water pressures are computed assuming a diffusive subglacial drainage system, and related to sliding rates via a Coulomb friction law. The present paper describes a major extension to this work, whereby the evolution of the subglacial drainage network is modelled as a morphological transition between a slow or distributed drainage system and a fast or channelized drainage system (e.g. Iken & Bindschadler 1986; Nienow et al. 2005). We demonstrate the enhanced capability of a model that includes both drainage system types using hypothetical drainage scenarios inspired by recent observations. These simulations comprise the seasonal drainage evolution of an Alpine-type glacier (§4a), summer drainage at a Greenland-type glacier (§4b) and a supraglacial lake drainage event (§4c). The significance of this work is in the application of a hydrologically coupled model to problems in which the evolution of the subglacial drainage system plays an important role. To date, hydrologically coupled models of ice dynamics have not included explicit treatment of interacting elements representing channelized (fast) and distributed (slow) drainage.

2. Background

Recent research has highlighted the key role of hydraulic processes in generating rapid changes to ice-sheet behaviour (e.g. Das et al. 2008; Stearns et al. 2008). Observations along the western margin of the Greenland ice sheet indicate ice-flow acceleration in summer occurring as a result of surface meltwater-enhanced lubrication of the glacier bed (e.g. Zwally et al. 2002; and if considering only regions dominated by ice-sheet flow; Joughin et al. 2008). Indeed, in a particular case, diurnal flow variations were associated with periodic changes in surface melting, whereas flow changes over a longer time period were linked to the episodic drainage of supraglacial lakes (Shepherd et al. 2009). Theoretical predictions (Alley et al. 2005; van der Veen 2007; Krawczynski et al. 2009) and observations (Boon & Sharp 2003; Das et al. 2008) have established that surface meltwater can penetrate thick subfreezing ice to provide lubrication to the glacier bed. Das et al. (2008) recorded short-lived vertical and horizontal accelerations at the glacier surface coinciding with a supraglacial lake drainage event in Greenland.

It is worth noting, however, that not all or even most of the recent large-scale acceleration of outlet glaciers in western Greenland (Rignot & Kanagaratnam 2006) can be attributed to meltwater-induced speed-up (e.g. Joughin et al. 2008). Other suggested mechanisms include a reduction in back-stress through the decrease in floating ice in front of the glacier (Thomas 2004; Nick et al. 2009), triggered by the transport of relatively warm water into the fjord (Holland et al. 2008; Rignot et al. 2010; Straneo et al. 2010). Increased surface melting may play an ancillary role in this process by enhancing hydro-fracturing in water-filled crevasses and hence calving rates at the terminus (Sohn et al. 1998).

3. Model

The hydrologically coupled higher order ice dynamics model used in this work is described in detail by Pimentel et al. (submitted) with the important extension that our treatment of the basal drainage system now includes both slow and fast components as described in Flowers (2008) (see also Hewitt & Fowler (2008) for a similar approach). Additional extensions to the previously published work are the use of elastic beam theory to model uplift and the inclusion of vertical fracture propagation. We consider a fixed geometry and do not solve for surface evolution. The core equations of the model system are presented below.

(a) Ice dynamics

Our model is based on the force-balance approximations of Blatter (1995) that are suitable for glacier flow regimes where both rapid and slow slip can occur (Schoof & Hindmarsh 2010). They include a multi-layer longitudinal stress scheme that has been given the classification LMLa among the so-called higher order models (Hindmarsh 2004; Pattyn et al. 2008). Longitudinal stress gradients can play a significant role, particularly when glaciers are subject to considerable slip (Price et al. 2002; Pattyn et al. 2008). The model is two dimensional (x-axis horizontal and orientated along the flowline and the z-axis vertical through the ice column), but aspects of three-dimensional flow have been parametrized; for example, lateral shear stress is parametrized by use of a flow-unit half-width, W. Glen’s flow law (Glen 1955; Nye 1957), describing the nonlinear rheological behaviour of glacier ice, is used in conjunction with the Stokes flow approximation (Blatter 1995) to formulate the horizontal stress balance as a second-order differential equation for velocity:Embedded Image3.1where the ice viscosity ν=ν(x,z) is given byEmbedded Image3.2where u=u(x,z) denotes the horizontal velocity, s=s(x) the surface elevation, uL=uL(x,z) the lateral velocity along the side wall of the glacier (assumed zero for this study), ρ the ice density, g the gravitational acceleration, A and n the parameters from Glen’s flow law (Glen 1955) and Embedded Image a small number to prevent a singularity in the event of vanishing strain rate (Meier 1958). At the ice–air interface, we assume a stress-free surface. At the ice–bed interface, we apply a basal friction law.

(b) Basal friction law

Modelling the basal sliding of glaciers has been a major area of research for several decades (e.g. Nye 1969, 1970; Kamb 1970; Iken 1981; Fowler 1986). More recent theoretical developments in glacier sliding theory have evinced that the commonly used sliding law, Embedded Images (with As,p,q being positive constants), relating basal drag τb to sliding velocity ub, fails to predict bounded basal drag for a given effective pressure N (Schoof 2005). The effective pressure, Embedded Image, is the difference between basal ice pressure and subglacial water pressure. Basal sliding with cavitation is therefore modelled using a regularized Coulomb friction law as proposed by Schoof (2005) and similarly by Gagliardini et al. (2007):Embedded Image3.3where C is a constant, Embedded Image a dimensional wavelength of the dominant bumps in the underlying bedrock and Embedded Image the maximum slope of bedrock obstacles. Sliding is strongly influenced by effective pressure, which itself is modulated by spatial and temporal variations in subglacial water pressure. Equation (3.3) is treated as a nonlinear Robin-type boundary condition for the problem, which because of its nonlinearity, forms part of an iterative solution procedure. Large amounts of slip can have a leading-order effect on the mechanics of the flow through the deviatoric normal stress gradients (Schoof & Hindmarsh 2010). An application of this friction law can also be used to estimate slip along the side walls, uL (see Pimentel et al. (submitted) for further details).

(c) Subglacial water sheet

The subglacial component of the hydrology model of Flowers & Clarke (2002) is used here, in one dimension, to describe the slow drainage system. Conservation of mass in a subglacial water sheet with a local subglacial water thickness, hs, can be writtenEmbedded Image3.4where qs denotes the vertically integrated subglacial water flux, QG is the geothermal flux, L is the latent heat of fusion of ice and Embedded Image a water source term for the direct input of surface melt in the absence of englacial transport. Adopting a Darcian framework to describe the flux,Embedded Image3.5with fluid potential ψs=Psw+ρwgb; K denotes the effective subglacial hydraulic conductivity and b the basal elevation. The subglacial water pressure, Psw, is ascertained by the empirically based expressionEmbedded Image3.6where hsc is the critical water thickness such that Psw=Pi (Flowers & Clarke 2002). hsc can be viewed as the areally averaged volume of a water-filled cavity system or as the product of porosity and undilated layer thickness of a macroporous till sheet separating the ice and bed. Although this empirical relation was derived from data specific to a soft-bedded alpine glacier (Flowers 2000), it serves the general purpose of describing a nonlinear increase in water pressure with water volume. A more general approach would avoid the prescription of a Embedded Image relation altogether (e.g. Hewitt & Fowler 2008).

(d) Subglacial conduits

Comprehensive treatment of the basal drainage system requires a representation of fast or channelized drainage, as well as the capability of dynamic evolution between fast and slow drainage systems. Although channels occupy only a small fraction of the glacier bed, they can regulate the pressure of their hydraulically connected surroundings (e.g. Hubbard et al. 1995). The physics of water flow through ice-walled conduits was first described by Röthlisberger (1972) and Nye (1976). The mathematical treatment of Spring & Hutter (1981, 1982) and the more recent update of Clarke (2003) are the foundation of our formulation here. We use a conduit density parameter to extend the treatment of an isolated conduit to a conduit system.

The discharge, Qc, within a saturated semicircular bed-floored conduit is expressed asEmbedded Image3.7where S is the conduit cross-sectional area and ψc=Pcw+ρwgb the fluid potential in the conduit system, with Pcw denoting the fluid pressure within the conduits. With radius R, we have a wetted perimeter Pwet=(π+2)R, and following Clarke (2003), a Darcy–Weisbach roughness fR=8gn2/R1/3H with prescribed Manning roughness n′ and hydraulic radius RH=πR/(2π+4).

Flowers (2008) parametrized the rate of water exchange between the sheet and the conduit systems asEmbedded Image3.8where hs:c is obtained from equation (3.6) using Embedded Image. We use dc to denote the latent conduit spacing and χs:c an exchange coefficient.

By considering mass conservation of ice and water within the conduit system, the following coupled partial differential equations arise:Embedded Image3.9andEmbedded Image3.10where A and n are Glen’s flow-law parameters, cp is the heat capacity of water, Φ is the Clausius–Clapeyron gradient and Embedded Image is a water source term to the conduit system. Numerical stiffness is encountered because the conservation equations admit both fast pressure-wave solutions as well as the desired more slowly varying water-flow solution. In order to suppress this stiffness and slow the pressure waves, a slight amount of conduit distensibility is introduced through the parameter γ (Clarke 2003).

Application of this mixed system of sheet and conduit interaction has been demonstrated by Flowers (2008). What is original in this deployment is the integrated approach of ice dynamics and hydrology, including both sheet and conduit systems.

(e) Uplift

When large amounts of water impinge on the glacier bed, high water pressures are generated and cause flexure of the overlying ice (e.g. Jóhannesson 2002). We attempt to model this effect through use of an elastic beam equation, which has a relatively long history of application in glaciology (e.g. Holdsworth 1969). Following Holdsworth (1969) and Vaughan (1995), we treat the glacier as a uniform static beam and use the Euler–Bernoulli beam equation,Embedded Image3.11where E is the elastic modulus, Embedded Image is the second moment of area for a glacier of rectangular cross section with Embedded Image a representative ice thickness; ω=ω(x) is the vertical deflection, and Embedded Image is the pressure load. Analytic solutions to equation (3.11) can be obtained under a variety of boundary conditions (e.g. Griffel 1970).

We assume the beam is simply supported at both ends, leading to the following symmetric boundary conditions: (d2ω/dx2)(L)=0, (d2ω/dx2)(0)=0, ω(L)=0 and ω(0)=0; L is the glacier length. The deflection or uplift is then computed asEmbedded Image3.12where a is the location of each point load. Other boundary conditions can easily be implemented, but are not essential when considering point loads on a long beam. We parametrize the increase in basal drainage system volume owing to uplift by making a correction to equation (3.6) and modifying the critical water sheet thickness, hsc(x)=hsc(x)+ω(x). In lieu of visco-elastic effects, we prescribe a Maxwell time of 6 h; considering a shear modulus μ of approximately 3.5×109Pa (Hobbs 1974, p. 258) and calculated effective viscosity ν between 106 and 107 Pa a, the estimated relaxation time τ=ν/μ is between 2.5 and 25 h (cf. note in Price et al. (2008)).

(f) Fracture mechanics

Crevasses are fractures in ice that form primarily as a result of tension. The influx of water to a crevasse can drive a fracture to propagate vertically through the glacier (e.g. Weertman 1973). A crevasse will penetrate to a depth at which the stress intensity factor at the crevasse tip is equal to the fracture toughness of the glacier ice. While the far-field tensile stress and fracture toughness determine where crevasses can form, the propagation of these crevasses is primarily governed by the rate of water injection. The theory of linear elastic fracture mechanics (LEFM), based upon the seminal work of Griffith (1921), has been applied to the fracture of ice (e.g. Weertman 1971a,b, 1973; Alley et al. 2005; van der Veen 2007; Krawczynski et al. 2009). Here, we use the LEFM model of van der Veen (2007) to model crack propagation. In this method, the stress intensity factor is used to describe elastic stresses near the tip of a crack.

The net stress intensity factor is estimated byEmbedded Image3.13where dcrev is the depth of the crevasse and bw is the water level within the crevasse. The first term is associated with crevasse opening (closing) owing to far-field tensile (compressive) normal stress Rxx, the second term corresponds to the closing effect caused by the lithostatic stress and the third term represents crevasse opening for a crevasse filled with water. The normal stress, Rxx, is defined as the total stress less the weight-induced lithostatic stress (ice overburden pressure):Embedded Image3.14The deviatoric stresses σxx′ and σyy′ are calculated by the ice-dynamics model (§3a). The water level relative to the bottom of the crevasse bw can be related to a water infilling rate Qin as bw=Qint. For our application, we assume the crack is kept water-filled; the crack geometry and lake volume constraints necessary for this assumption have been studied by Krawczynski et al. (2009). For a given fracture toughness, KI=KIC, equation (3.13) can be solved to find the depth, dcrev, of the crevasse. In this study, the method is used to determine whether the stresses and water inputs are sufficient to drive a crevasse through the entire ice thickness needed to initiate supraglacial lake drainage (§4c, experiment 3). This method does not attempt to predict crack growth rate (cf. Tsai & Rice 2010).

4. Simulation design

We describe three tests of the model: experiment 1, an ordinary Alpine-type glacier seasonal cycle; experiment 2, the response of a Greenland-type outlet glacier to summer drainage; and experiment 3, a subglacial flood resulting from the drainage of a supraglacial lake.

(a) Experiment 1: melt season transition of an Alpine-type glacier

An initial simulation is conducted to demonstrate that the hydrologically coupled ice-dynamics model reproduces the qualitative features of glacier seasonal transitions. This would include enhanced basal flow, vertical uplift and the emergence of a channelized drainage system (e.g. Iken et al. 1983; Mair et al. 2003; Anderson et al. 2004). A synthetic air-temperature time series (figure 1) covering a hypothetical melt season is applied to an idealized mountain glacier profile (figure 1). The time series represents temperature at the mean glacier-surface elevation. Using an atmospheric lapse rate of −5.9°C km−1 and a degree day factor for ice of 6.2 mm w.e.°C−1 d−1 (within the typical range found in mountain areas (cf. table 1 in Hock 2003), surface melt rates are determined along the glacier flowline; this water is routed directly to the bed through the term Embedded Image in equation (3.10).

Figure 1.

(a) Alpine-type glacier geometry (width fixed at 1000 m) and (b) synthetic air-temperature time series generated by superimposing a diurnal sinusoid of 5°C amplitude on an annual sinusoid of 15°C amplitude and a mean of −10°C.

View this table:
Table 1.

Model constants.

(b) Experiment 2: summer drainage at a Greenland-type outlet glacier

In this example, we draw upon a recent observational study conducted at Russell Glacier in southwest Greenland (Shepherd et al. 2009). Using measurements of air temperature and surface displacement at three locations along the glacier flowline, Shepherd et al. (2009) documented a link between surface melting and enhanced flow speeds over a 6 day observation period. Taking estimates of daily meltwater production from Shepherd et al. (2009), we endeavour to replicate the reported fluctuations in ice speed and height change of late summer (Shepherd et al. 2009). At 37 km from the glacier terminus, ice speeds determined from interferometric synthetic aperture radar (InSAR) are shown to increase by 150 m a−1 from the winter baseline to the summer peak (Shepherd et al. 2009). Furthermore, diurnal variations in ice speed at this location reached a maximum amplitude of 110 per cent of the baseline value and height changes of up to 4 cm were recorded (Shepherd et al. 2009). The glacier geometry is not known, so we construct an idealized glacier profile (figure 2a) that is based upon reported ice thicknesses of 890, 1050 and 1120 m at the three global positioning system (GPS) measurement locations sited 37, 53 and 72 km inland of the ice margin, respectively (Shepherd et al. 2009).

Figure 2.

Synthetic surface and bed elevations for (a) the land-terminating Russell Glacier simulations and (b) the tidewater-terminating Greenland outlet glacier simulations.

Shepherd et al. (2009) used a positive degree-day model to obtain estimates of meltwater production between 0.02 and 0.04 km3 d−1 during the observed summer period. Over the 4100 km2 catchment, this amounts to supraglacial runoff in the range 0.5−1 cm d−1. In order to obtain approximate seasonal and diurnal forcing, we idealize the situation by applying a seasonal sinusoid such that these derived values represent the peak summer melt rates at the upper and lower limits of the glacier flowline. The daily meltwater amounts are then distributed over a prescribed sinusoidal diurnal cycle. We assume that water is routed directly to the bed (by means of source term Embedded Image in equation (3.10)). Shepherd et al. (2009) note that the subglacial drainage network received water, up to approximately 0.05 km3, on an episodic basis; this water originates from meltwater stored in supraglacial lakes predominantly at high elevations. For simplicity, we have neglected this component of the drainage and will treat supraglacial lake drainage events independently in the next example.

(c) Experiment 3: supraglacial lake drainage in western Greenland

Finally, we explore the simulated dynamical response of a Greenland outlet glacier to subglacial flooding, as documented by Das et al. (2008). The rapid supraglacial lake drainage event produced increased seismicity, transient acceleration, rapid ice-sheet uplift (1.2 m) and horizontal displacement (0.8 m) followed by subsidence and deceleration over 24 h (Das et al. 2008). The outlet glacier in question terminates in tidewater, but the floodpath geometry is not known. We assume a simple geometry where the glacier rests on a horizontal bed and increases in thickness from 500 m at the terminus to 1200 m at a distance of 50 km (figure 2b). The boundary condition at the ice–ocean interface is given in Paterson (1994) as a balance of the longitudinal stress gradient and the hydrostatic pressure of the ocean water:Embedded Image4.1Also at the terminus we have Psw=Pcw=Pi.

We simulate a supraglacial lake drainage event according to the specifications given by Das et al. (2008): the lake of volume 0.044 km3 drains through 980 m of ice and empties in 1.4 h. The LEFM model (§3f) is used to propagate a crack beneath the lake until a surface to bed connection is achieved and lake drainage is initiated. Assuming a sinusoidal shape for the lake drainage hydrograph, the rate of lake-water input to the glacier is prescribed asEmbedded Image4.2where hξf is the horizontal grid cell length. We prescribe this rate of input for 0≤t≤5040 s at the location along the flowline where the ice is 980 m thick. We restrict ourselves to simulating conditions downstream of the water influx location, although the response is not limited to this area. An equilibrium model state is obtained before the flood is initiated. The simulation is then run for a total of 30 days to monitor transients in the aftermath of the flood.

Das et al. (2008) hypothesize that the short-lived ice-dynamic response to the flood may have been due to the presence of an efficient subglacial hydrological system. In order to make an inference about the configuration of the basal hydrological system prior to the drainage event, we first include only the slow (distributed) drainage system in the model. For this simulation, lake-water input is prescribed in the same manner as equation (4.2), but through Embedded Image (equation (3.4)) rather than Embedded Image (equation (3.10)). We then add the fast (channelized) drainage system and experiment with the existence of steady-state channels of various sizes prior to the flood. Assuming a steady supply of background water to the bed, arising either from regular surface water drainage (as, for example, indicated in the previous section) or from multiple previous lake drainage events, we prescribe a continuous background water discharge into the conduit system of Qcin=10 and 100 m3 s−1. We also model the situation in which conduits are not pre-developed, using Qcin=0. Qcin is applied as an upstream flux boundary condition. As before, the flood is not initiated until an equilibrium model state is obtained.

(d) Parameter choices

Constants and reference parameters common to all simulations are given in table 1, while those relevant to individual simulations are provided in table 2. An explanation of the key choices follows; for more detailed discussion of parameter selection and sensitivity, see Pimentel et al. (submitted) and Flowers (2008).

View this table:
Table 2.

Model parameters.

The rate factor in Glen’s flow law, A, is associated with the rheological softness of the ice and is dependent on the ice temperature. Softer (warmer) ice results in greater deformation rates and sliding enhancement in the model. Using a modified Arrhenius relation (Paterson & Budd 1982), we choose a value corresponding to an ice temperature of −2°C for our Alpine-type example (e.g. Paterson 1994) and −5°C for our Greenland-type examples.

Transmissivity of the subglacial drainage system is a product of the effective subglacial hydraulic conductivity K and the critical sheet thickness hsc. To accommodate the scale dependence of subglacial transmissivity, we hold K fixed and adjust the critical sheet thickness for different simulations (hsc=0.15 m for the Alpine-type glacier and 5 m for the Greenland-type outlet glaciers, increasing according to the scale of the system). In our experience, this adjustment to drainage system capacity becomes necessary to simulate realistic subglacial water pressures under ice sheets. The hydraulic conductivity sets the flow speed for a given fluid potential gradient, and a constant value of 0.1 m s−1 is chosen which is within the range of other numerical studies (e.g. Flowers 2008). Increasing K and or hsc would prolong the transition from sheet- to conduit-dominated discharge because this effectively increases the capacity of the slow drainage system. The Manning roughness n′, which governs turbulent heat transfer, can also be increased to inhibit conduit growth and diminish the role of channelized drainage in the subglacial system (cf. Clarke 2003). The initial condition for conduit cross-sectional area S(t=0)=Sϵ represents a vestigial channel (i.e. S is not permitted to reduce below this value). We choose a larger value for the Greenland-type examples (10−3 m2 compared with 10−4 m2). In all cases, we specify the latent conduit spacing, dc.

Some degree of uncertainty surrounds the appropriate value of the elastic (Young’s) modulus, E, used in the beam equation (3.11). Values derived from tidal deflection records of ice shelves are as low as 0.9 GPa, whereas laboratory experiments with ice give E=9.3 GPa (Vaughan 1995; Reeh et al. 2003). We select a value of E=3.9 GPa for our simulations. Determining a relaxation time scale is also problematic. In both respects, a visco-elastic uplift model may therefore be more appropriate. This indeed is the reasoning of Reeh et al. (2003) when examining the tidal flexure of Nioghalvfjerdsfjorden glacier in northeast Greenland, and Evatt & Fowler (2007) when investigating cauldron subsidence on the surface of Vatnajökull ice cap in Iceland.

5. Results

(a) Experiment 1: melt season transition at an Alpine-type glacier

Evolution of the subglacial drainage system and glacier response are modelled over a complete melt season (figure 3). Discharge from the subglacial drainage system is separated into contributions from the sheet (figure 3a) and the conduit (figure 3b) systems. As would be expected early in the melt season, discharge from the distributed system (sheet) dominates, but as we enter the spring season, large amounts of water give rise to channelized flow and a transition occurs such that conduit discharge becomes dominant. This switch, entirely emergent from subglacial processes, occurs around day 152 in the simulation and reflects the commonly observed spring transition in the glacier drainage network where the fast (conduit) system grows at the expense of the slow (sheet) system as meltwater flux to the glacier bed increases (e.g. Nienow et al. 1998). The conduits enlarge and grow headward until the peak of the melt season, after which creep closure overtakes melt-opening and the conduit cross-sectional area decreases over most of the glacier. The cross-sectional area reaches a peak of 2 m2 near the terminus (figure 3c). In figure 3d, we see a maximum in basal water pressures whereby sections of the glacier, between days 145 and 165, approach or exceed flotation. Here, the water sheet exceeds saturation and water is diverted into the developing conduits. Once the conduits have grown to a sufficient size, such that they discharge most of the subglacial water (by approx. day 152), basal water pressures return to previous values.

Figure 3.

Position–time diagrams of key variables for simulation of a seasonal transition at an Alpine-type glacier: (a) sheet discharge, Qs; (b) conduit system discharge, Qc; (c) conduit cross-sectional area, S; (d) subglacial water pressure, Psw; (e) basal ice-flow speed, ub; ( f) basal shear stress, τb; (g) basal longitudinal stress, σxx(b); and (h) vertical uplift, ω. Diurnal variability is evident in these figures as horizontal striping.

While water pressures are increasing, the glacier markedly accelerates from its quiescent baseline (figure 3e), with basal speeds topping 50 m a−1 at their peak. The model therefore qualitatively captures the seasonal speed-up, with hydrologically driven accelerations clearly manifested. The effects of diurnal forcing on the glacier are also present, with diurnal modes of variability evident in the output fields. A clear correspondence can be identified between the zone of high basal water pressures (figure 3d) and low basal shear stress (figure 3f). When flotation occurs, basal traction is lost and free slip ensues. Vanishing basal shear stress is compensated by larger longitudinal stresses (figure 3g) and increased lateral drag (not shown), demonstrating the importance of longitudinal and lateral stresses in the force balance. Finally, hydraulic uplift occurs in the model when flotation is reached (figure 3h). In this example, uplift reaches a maximum value of approximately 1.3 cm, providing a small increase in the volume of the distributed drainage system when water pressures are greatest. This magnitude can be controlled to some extent by the value of the elastic modulus.

(b) Experiment 2: summer drainage at a Greenland-type outlet glacier

Prescribed meltwater production rates, with diurnal and seasonal variability, are used as influx to the modelled sheet-conduit subglacial system (figure 4a). Our simulations of the seasonal cycle show a speed-up, from a winter baseline (mean winter surface speeds) to peak summer values, of 160 m a−1 at 37 km from the terminus (figure 4b). This compares well with the results of Shepherd et al. (2009) who use InSAR to derive seasonal fluctuations in ice speed of around 150 m a−1 close to site 1 (37 km from the terminus). During the late summer season, covered by the GPS observation period (days 201–206), Shepherd et al. (2009) observe an average ice speed of 166 m a−1 at site 1; this compares well with our simulated speed of 162 m a−1 at the same position and over the same time period. Our numerical simulations produced negligible diurnal variability in glacier speed at measurement sites 2 and 3 (53 and 72 km from the terminus, respectively) and only minor diurnal variability less than (10 m a−1) at site 1 (37 km from the terminus). A stronger simulated diurnal cycle emerges 10 km further downglacier (figure 4c), but still not of the scale measured at site 1 where daily speed-up averaged 55 per cent of baseline values (figure 2a in Shepherd et al. 2009). After correcting for downslope motion, we calculate diurnal vertical deflections of the ice surface (magnitude approx. 2 cm, see figure 4d) that are consistent with the observed range. The measurement sites exhibited a 1–4 cm diurnal uplift of the ice surface (fig. 2b in Shepherd et al. 2009). Note that we are also assuming here that our calculated basal deflection is transmitted to the ice surface without damping or spreading.

Figure 4.

(a) Total (glacier wide) subglacial water influx, showing diurnal variations. (b) Seasonal speed-up from winter baseline using daily mean surface glacier speeds at 37 km from the terminus. (c) Glacier-surface speed at 27 km from the terminus. (d) Vertical motion at 37 km from the terminus.

(c) Experiment 3: supraglacial lake drainage in western Greenland

Simulations of this tidewater-terminating glacier that employ only the slow (distributed) drainage system produce unrealistically high water pressures resulting from the subglacial floodwater prescribed at the glacier bed (figure 5a). In this case, the subglacial deluge causes a spike in water pressure, reaching almost eight times the ice overburden at the source. It takes 9 days for water pressures along the flowline to subside below flotation. The further decline in pressure exhibits a long tail, with water pressures remaining slightly higher than their pre-flood values even after 30 days. Including the fast (channelized) drainage system and injecting water into the conduit system at the bed creates a different dynamic response (figure 5b). As would be expected, the maximum water pressures are reduced with the introduction of an efficient drainage system. However, with a vestigial conduit cross-sectional area of 10−3 m2 specified in the model prior to the flood, maximum water pressures still exceed twice overburden pressure in the early stages of the flood. Although the maximum pressures are reduced in this simulation, the disturbance to the basal drainage system remains long-lived. The observations of Das et al. (2008) suggest a much shorter lived glacier response to the lake drainage (the major transient response occurred within a day), raising the question of whether the flood may have encountered a pre-existing efficient drainage system at the glacier bed.

Steady-state channels of various sizes were introduced prior to the modelled flood by applying a background discharge into the conduit system (§4c). When pre-existing conduits are sufficiently large, they substantially reduce the duration of the flood response (figure 5c,d). The greater the prescribed conduit discharge, the larger the equilibrium cross-sectional area of the conduits prior to flood, hence the lower the peak water pressure and faster the egress of water from the flood-initiation point. Background forcing of Qcin=100 m3 s−1 produces a conduit cross-sectional area, averaged along the flowline, of 22 m2 with which to initiate the flood simulation. This scenario most closely resembles the situation described in Das et al. (2008) in terms of the duration of the flood response (figure 5d). The existence of conduits of this size seems plausible, given our simulations from §5b where conduit cross-sectional areas greater than 15 m2 extended 28 km up-glacier at the peak of the melt season. However, conduit cross-sectional area depends on the choice of the latent conduit spacing dc.

Figure 5.

Evolution of modelled subglacial water pressures (given as a flotation fraction Psw/Pi) plotted at 5 km intervals (spike is closest to flood injection point). (a) Sheet drainage only, (b) mixed (sheet-conduit) drainage, (c) mixed drainage with a prescribed background discharge of Qcin= 10 m3 s−1 and (d) Qcin=100 m3 s−1. Note the difference in vertical scales between the panels.

Using the case of Qcin=100 m3 s−1, we now look in detail at the modelled dynamic response of the ice during the 24 h period following the flood (figure 6). The large volume of water impinging on the bed requires the efficient drainage network to dominate proceedings, with discharge in the conduit system dwarfing that in the sheet (compare figure 6a,b). Conduits in the vicinity of the water source grow rapidly from an initial area of 25 m2 to exceeding 300 m2 in cross section (figure 6d). As the subglacial system adapts and water is evacuated, subglacial pressures are steadily reduced. Only a small fraction of the glacier bed remains at a high flotation fraction after 24 h (figure 6c). A consequence of the excess pressure beneath the glacier is an upward deflection of the ice. Vertical uplift of the subglacial roof reaches 0.5 m in this simulation (figure 6e), less than the 1.2 m surface deflection recorded at the lake shoreline (Das et al. 2008). Low basal resistance (figure 6g) arising from increased lubrication results in increased basal sliding speeds (figure 6f). A transition zone is created between the area of zero traction upstream (where fast flow occurs) and the area of increased basal resistance downstream; during this phase, a ridge of compressive stresses recedes up-glacier (figure 6h). Longitudinal (tensile) stresses are high near the terminus throughout the simulation (figure 6h). Increased ice speed generates greater shearing and the parametrized lateral drag contributes significantly to the force balance during the flood (figure 6i). Persistent features visible near the terminus (figure 6) are related to our tidewater boundary conditions.

Figure 6.

Position–time diagrams, displayed with distance downglacier from flood source and time elapsed since flood initiation, of key model variables: (a) Water sheet discharge, Qs; (b) conduit system discharge, Qc; (c) subglacial water pressure, Psw; (d) conduit cross-sectional area, S; (e) vertical uplift, ω; ( f) basal speed, ub; (g) basal shear stress, τb; (h) vertically averaged longitudinal stress, Embedded Image; and (i) vertically averaged lateral shear stress, Embedded Image.

6. Discussion

Ice motion on Alpine glaciers is known to be influenced by subglacial drainage system transitions over the melt season, and numerous observational studies have documented hydromechanical coupling over diurnal and seasonal time scales (e.g. Iken et al. 1983; Mair et al. 2003; Anderson et al. 2004). The switch from a system of low hydraulic capacity and efficiency to one of high capacity and efficiency has been qualitatively replicated by the model (§5a). The simulations successfully capture the conceptual features of subglacial drainage evolution and hydrologically modulated mountain glacier flow. Parallels between Alpine glaciers and Greenland in the way that large amounts of meltwater facilitate basal motion have been recently reported by Bartholomew et al. (2010), and suggest a unified treatment of basal processes across a range of scales (Parizek 2010).

The remainder of our discussion focuses on aspects of hydrologically controlled ice dynamics occurring on the Greenland ice sheet (§5b,c). Without detailed input data, such as an accurate record of the spatial and temporal variability of meltwater production and infiltration over the drainage catchment, a quantitative comparison between our simulations and the findings reported in Shepherd et al. (2009) is difficult. By making simplistic assumptions about the distribution and injection of reported meltwater volumes, a good deal of broad agreement between the simulation and observations is achieved. For example, seasonal accelerations are captured to within 10 m a−1. Diurnal accelerations of the magnitude reported in Shepherd et al. (2009) are not replicated. Such strong fluctuations can be produced by the model (see, for example, §5a), suggesting that the weakness of the diurnal variations in §5b is not necessarily a result of a fundamental model deficiency, but rather a product of parameter choices or the uncertainty in the forcing. Improvements in the comparability of the simulations should result from using measured air temperatures, observed lapse rates and site-specific degree-day factors, as well as a more precise description of the glacier geometry. The absence of the episodic lake drainage in the model is also a contributing factor. It is clear from our lake drainage simulation (§5c) that this would have an effect in priming the channels and increasing drainage efficiency.

The lake drainage event reported by Das et al. (2008) is clearly three dimensional in nature: it is described as moving the ice surface up and away from the lake centre. Horizontal northward motion of 0.8 m was measured at the lake shoreline, and this rapid northward motion was orthogonal to the regular westerly glacier velocity that remained relatively stable during the flood (fig. 2 in Das et al. 2008). The two-dimensional nature of our model is therefore limited in its ability to replicate the complex dynamics involved. The likely multi-directional flow of the floodwater as it impinged on the glacier bed (contrasted to our one-directional assumption) also makes a detailed quantitative comparison impossible. Our simulated vertical deflection was just over a factor of 2 smaller than the observed surface uplift (1.2 m; Das et al. 2008). More realistic solutions could be obtained from a visco-elastic plate rather than the less well-suited elastic beam implemented here. In another modelling study involving this drainage event, Tsai & Rice (2010) consider how the turbulent flow of the subglacial flood water opens a basal crack between ice and bedrock. Their modelling approach to this flood also underestimated surface displacements (likewise by a factor of 2–3). Although we are limited in the conclusions that can be drawn from our modelling exercise, our results do support the role of a pre-existing channel network in dissipating the response to the flood as quickly as was observed. One possible interpretation concerning this requirement is that the well-established channelized network hinders our ability to simulate the observed magnitudes of horizontal and vertical displacement. If this is the case, then incorporating horizontal turbulent hydraulic fracture for basal crack propagation (Tsai & Rice 2010) could be the mechanism that creates rapid initial channel growth without undermining the simulated surface displacements.

7. Conclusions

Satellite and direct observations of supraglacial lakes forming and draining seasonally over the Greenland ice sheet have been documented (Box & Ski 2007; McMillan et al. 2007) and close monitoring at particular sites (e.g. Das et al. 2008; Shepherd et al. 2009) reveals a direct link between drainage and glacier acceleration and uplift. The numerical model described in this paper endeavours to qualitatively replicate such behaviour by coupling higher order ice dynamics to a subglacial drainage system that enables switching between a slow distributed network and fast channelized water flow. Qualitative features of hydrologically regulated ice flow are reproduced in our numerical simulations on several spatial and temporal scales. Key features of the classical seasonal transition of an Alpine-type glacier are captured, including the diurnal and seasonal fluctuations of glacier flow speed, as well as the so-called spring transition whereby the fast conduit drainage system grows at the expense of the slow distributed drainage system as meltwater flux to the glacier bed increases (e.g. Nienow et al. 1998). Simulations of summer drainage at a Greenland-type outlet glacier based upon observations from Shepherd et al. (2009) produce modelled flow speeds that compare well with the measured seasonal speed-up. Both the mean glacier-surface speeds and the diurnal surface height changes of late summer are simulated with appropriate magnitudes, but the model fails to produce the strong diurnal cyclicity observed in the horizontal flow speed. This is probably a result of ill-constrained spatial and temporal variability in subglacial meltwater input. Simulations of a subglacial flood resulting from a supraglacial lake drainage event of the scale recorded by Das et al. (2008) produce features (rapid horizontal and vertical accelerations) broadly inline with the observations insofar as two-dimensional model results can be compared with observations of three-dimensional processes. The modelled glacier response to the flood can be confined to the short time scale observed (approx. 24 h) when the prescribed flood impinges on a pre-existing channel system. These subglacial channels are likely developed earlier in the season through both regular meltwater drainage (a view consistent with results from §5b) and previous lake drainage events within the catchment.


This work has been funded by Natural Sciences and Engineering Research Council (NSERC) of Canada as part of the Canadian contribution to International Polar Year project GLACIODYN (IPY30). The Canadian Foundation for Climate and Atmospheric Science (CFCAS) provided additional funding for S.P. through the Polar Climate Stability Network (PCSN) while working at the University of British Columbia. The authors thank Christian Schoof for his advice on sliding laws and parametrizing lateral shear stress.


  • Present address: Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada.

  • Received April 20, 2010.
  • Accepted July 5, 2010.


View Abstract