A spectrum of jökulhlaup dynamics revealed by GPS measurements of glacier surface motion

Abstract

GPS campaigns on glaciers during jökulhlaups show how subglacial floods affect glacier motion and shed light on the dynamics of such floods. Three such campaigns have been carried out on southern and western Vatnajökull, southeast Iceland, over known jökulhlaup paths. Two slowly rising jökulhlaups from Grímsvötn and two rapidly rising jökulhlaups from the western and eastern Skaftá cauldrons were captured in these campaigns, with maximum discharge ranging from 240 to 3300 m³ s⁻¹. Glacier surface movements measured in these campaigns are presented along with the corresponding discharge curves. The measurements are interpreted as indicating: (1) initiation of rapidly rising jökulhlaups with a propagating subglacial pressure wave, (2) decreased glacier basal friction during jökulhlaups, (3) subglacial accumulation of water in slowly rising jökulhlaups and (4) lifting of the glacier caused by subglacial water pressure exceeding overburden in both rapidly and slowly rising jökulhlaups. The latter two observations are inconsistent with assumptions that are typically made in theoretical and numerical modelling of jökulhlaups. Both viscous and elastic deformation of the glacier as well as turbulent hydraulic fracture at the ice/bedrock interface are important in the dynamics of the subglacial pressure wave at the front of rapidly rising jökulhlaups.

INTRODUCTION

Jökulhlaups (glacier outburst floods) from subglacial geothermal areas, marginal lakes and subglacial volcanic eruptions are common in Iceland (e.g. Björnsson, 1975, 1976, 2002). These floods overfill the subglacial hydrological system and have a substantial effect on the dynamics of the glaciers, where they occur (Anderson and others, 2005; Magnússon and others, 2007; Sugiyama and others, 2007; Bartholomaus and others, 2008, 2011; Stearns and others, 2008; Riesen and others, 2010). Jökulhlaups can be categorized into two main groups, rapidly and slowly rising (Björnsson, 2002), which are characterized by marked differences in the hydrographs.

The hydrographs of slowly rising jökulhlaups have been explained with a model of melt enlargement of an ice tunnel (Nye, 1976) or melt enlargement of a constriction (bottleneck) at the location of a seal near the source lake (Clarke, 1982; Björnsson, 1992). Clarke (2003) developed the theory of Nye (1976) further based on Spring and Hutter (1981, 1982) and showed that the effective pressure (ice overburden minus water pressure) in a flood through a subglacial tunnel may become negative in parts of the flood path. Negative effective pressure is expected to lead to water escaping the tunnel due to hydraulic fracturing at the ice/bedrock interface and lifting of the glacier, creating a wider and more complex flood path at the corresponding location. Glacier uplift and widespread velocity increase have been observed for slowly rising jökulhlaups from the subglacial lake Grímsvötn in 1996 (Magnússon and others, 2007) and the marginal lake Gornersee, Switzerland, in 2005 (Huss and others, 2007), indicating that the flood path of slowly rising jökulhlaups is not always restricted to one or more tunnels. Björnsson (2010) therefore states that ‘current models may reconstruct discharge curves while not describing all the factual hydraulic and glaciodynamic processes of each jökulhlaup’.

In contrast, it has been proposed that the flood path of rapidly rising jökulhlaups is created by the propagation of a subglacial pressure wave (Jóhannesson, 2002) at the front of sheet flow underneath the glacier (Björnsson, 2002; Flowers and others, 2004). During rapidly rising jökulhlaups from Grímsvötn and from the lake beneath the western Skaftá ice cauldron, the estimated water volume stored in the flood paths was an order of magnitude larger than could be stored in conduits formed by melting. Lifting and deformation of the overlying glacier was, therefore, important in their dynamics (Björnsson, 2002). Indications of subglacial water pressure higher than ice overburden, and associated glacier lifting have been observed for rapidly rising jökulhlaups from Grímsvötn in 1996 (Björnsson, 1997; Roberts and others, 2000) and for Gornersee in 2004 (Huss and others, 2007; Sugiyama and others, 2008). Indications of subglacial sheet flow in the early stages of a rapidly rising jökulhlaup from the eastern Skaftá cauldron have also been observed (Magnússon and others, 2007), but direct observations of the propagation of a subglacial pressure wave have not been reported.

The propagation mechanism of jökulhlaups at the glacier bed is poorly understood. Continuous GPS measurements on the glacier surface can provide valuable information on the effect of jökulhlaups on the motion of the glacier and shed light on jökulhlaup dynamics. Three GPS campaigns have been carried out on southern and western Vatnajökull, Iceland, over jökulhlaup paths from subglacial lakes at Grímsvötn and the Skaftá cauldrons. Glacier surface movements during jökulhlaups from these campaigns are presented and interpreted in this paper. Similarities and differences in the initiation of rapidly and slowly rising jökulhlaups are addressed based on observations of both types of jökulhlaups. Based on these observations, we argue that jökulhlaups show a spectrum of behaviour instead of being classifiable into two distinct categories, although, the traditional categories remain useful in discussions of flood characteristics. We identify physical processes that we infer to be acting at the glacier bed and within the glacier during the floods and compare our observations with the conceptual background of existing models.

Subglacial water flow and variations in subglacial water pressure have attracted increasing attention in recent years as a likely cause of large variations in ice-flow velocities that have been observed on the main outlets of the Greenland ice sheet and some of the ice streams of Antarctica (e.g. Rignot and Kanagaratnam, 2006; Fricker and others, 2007; Stearns and others, 2008; Doyle and others, 2015). The focus of this paper is on subglacial water flow in the context of jökulhlaup dynamics, but the subject has a wider relevance, as improved understanding of subglacial hydrology is important for understanding the response of glaciers and ice sheets to hydrological variations such as changes in meltwater input due to climate change.

THE STUDY SITES AND THE INVESTIGATED JÖKULHLAUPS

This study was conducted over two jökulhlaup paths in the Vatnajökull ice cap, at the Skeiðarárjökull and Skaftárjökull outlet glaciers in the southern and western parts of the ice cap, respectively (Fig. 1).

Fig. 1. Location of the subglacial lakes at Grímsvötn and the Skaftá cauldrons. Estimated flood paths, based on the gradient of the hydraulic potential, and location of GPS stations, discharge measurement sites and hydrometric stations.

Jökulhlaups from Skeiðarárjökull originate in the Grímsvötn subglacial lake (Fig. 1) (Björnsson, 1974; Thórarinsson, 1974). The lake is formed by geothermal melting in a ~10 km wide and ~300 m deep subglacial volcanic caldera (Björnsson, 1974). The lake accumulates basal meltwater produced by geothermal and volcanic activity in addition to surface meltwater and rain (Björnsson and Kristmannsdóttir, 1984). Most jökulhlaups from Grímsvötn are classified as slowly rising (Björnsson, 2002) and the exponential rise of the discharge during these jökulhlaups has been physically explained by invoking the conduit-melt–discharge feedback (Nye, 1976; Björnsson, 1992). A jökulhlaup in November 1996 that rose to a maximum discharge of ~50 000 m³ s⁻¹ at the terminus in a little over a day after outflow from Grímsvötn started (Björnsson, 1997, 2002; Snorrason and others, 1997) is, however, an example of a rapidly rising jökulhlaup (Björnsson, 2002; Jóhannesson, 2002). This jökulhlaup released ~3 km³ of meltwater from an eruption in Gjálp north of Grímsvötn (Guðmundsson and others, 1997) and damaged the seal of the Grímsvötn subglacial lake, altering the characteristics of jökulhlaups after 1996. Less water can now accumulate in Grímsvötn, hence the volume of jökulhlaups is smaller. At the same time, their duration became shorter. The peak flow of the two largest post-1996 jökulhlaups, in the autumn of 2004 and 2010, ~3000 m³ s⁻¹, is therefore, comparable with typical jökulhlaups prior to 1996. A GPS station was operated on Skeiðarárjökull during the 2004 and 2010 jökulhlaups. The station was located ~1 km up-glacier of the outlet of the ~50 km long flood path in 2004 (denoted by SKE2) and ~9 km up-glacier of the outlet in 2010 (denoted by SKE1) (Figs 1 and 2). Previous studies have shown water accumulation at SKE1 during a jökulhlaup and also during an intense rain event (Magnússon and others, 2011).

Fig. 2. Longitudinal profiles of bedrock and ice surface along jökulhlaup paths from Grímsvötn (a) and the western Skaftá cauldron (b) (see Fig. 1). The lower 32 km of the flood paths from the western and eastern Skaftá cauldrons are identical and the elevation profiles for the upper parts are so similar that they are nearly indistinguishable at the scale of this figure. Therefore, only the path from the western cauldron is shown in (b). The locations of the GPS stations are indicated. The elevation of the GPS stations corresponds to the glacier surface at the time of each measurement campaign. Bedrock elevation is based on radio-echo soundings (Björnsson, 1988; Björnsson and others, 1999; Magnússon, 2008; Magnússon and others, 2009).

Jökulhlaups in the Skaftá river from western Vatnajökull have been known since 1955 (Björnsson, 1977). These jökulhlaups originate from two subglacial lakes under 1–3 km wide and 50–150 m deep depressions in the glacier surface, commonly referred to as ice cauldrons, formed by geothermal melting. The average time interval between jökulhlaups from each cauldron is ~2 a. The jökulhlaups travel ~40 km under a glacier that reaches maximum thickness of ~750 m and emerge in the Skaftá river at the terminus of Skaftárjökull outlet glacier (Figs 1 and 2). A hydrometric station is located in the river by Sveinstindur, 25 km downstream of the glacier margin (Fig. 1). The hydrographs of all measured jökulhlaups from the Skaftá cauldrons have been rapidly rising (Zóphóníasson, 2002). Three continuous GPS stations were installed over the flood path of jökulhlaups from the Skaftá cauldrons in 2008, located 3, 8 and 15 km from the glacier margin and denoted D3, D8 and D15, respectively. Two jökulhlaups occurred during the GPS campaign, one with a maximum discharge of 240 m³ s⁻¹ from the western cauldron in August 2008 and the other with a maximum discharge of 1290 m³ s⁻¹ from the eastern cauldron in October 2008.

The jökulhlaup paths from Grímsvötn and the Skaftá cauldrons have similar length and ice thickness, whereas the vertical relief is larger for Grímsvötn (Fig. 2). The Grímsvötn ice-flow basin (~300 km²) is an order of magnitude larger than the ice-flow basins of each of the Skaftá cauldrons (~20 and ~30 km² for the western and eastern cauldron, respectively) and accumulates an order of magnitude more meltwater annually in the subglacial lake (Björnsson, 1974, 1977, 2002; Pálsson and others, 2006). The jökulhlaups from Grímsvötn have been more varied in size and characteristics than jökulhlaups from the Skaftá cauldrons, some released in connection with subglacial volcanic eruptions and others caused by the continuous accumulation of water in the lake. The jökulhlaups from the Skaftá cauldrons, on the other hand, vary in size over an order of magnitude and are all caused by slow accumulation of water from geothermal melting, surface melting and rain.

METHODS

Instrumentation

The GPS instruments used to monitor the two events from the Skaftá cauldrons and the 2010 event from Grímsvötn were installed on top of the glacier surface on quadropods, while the GPS instrument used in the 2004 event from Grímsvötn was installed on a pole drilled into the ice. The height of the instruments deployed on the quadropods is affected by surface melting as they are lowered with the surface. Observed elevation changes at all instruments are affected by the local surface slope and the down-glacier movement of the ice. These processes are not corrected for and need to be accounted for, when the elevation data are studied over periods of days or periods with large horizontal movements.

The GPS instruments used during recording of all events were Trimble NetRS (Trimble 5700 in 2004) dual-frequency receivers recording at 15 s intervals. At the time of October 2008 jökulhlaup, the D15 instrument at the highest elevation was snowed in and out of power. The D8 instrument at the second highest elevation lost power on 14 October. The maximum discharge during the October event had been reached at that time and recession had started but was not fully completed. SKE2, the station running in 2004 on Skeiðarárjökull, was taken down on 6 November, while the glacier was still subsiding after the 2004 jökulhlaup from Grímsvötn.

Data processing

The GPS data were processed with the GAMIT-Track utility (Herring and others, 2010) using a set up for long baselines. The August and October events on Skaftárjökull in 2008 were processed with respect to a base station at Skrokkalda (SKRO) in August and Ísakot (ISAK) in October, giving baselines of ~35 and ~85 km, respectively. The stations operated on Skeiðarárjökull were processed using Höfn (HOFN) as base station at ~100 km distance. The standard deviation of unfiltered positions around a daily mean during periods of slow motion and low melt is ~4 cm in the vertical and 2 cm in the horizontal coordinates and may be used as an indication of the precision of the GPS measurements. Horizontal velocities were calculated for 1 h periods using filtered positions (2 h wide triangular filter). During days prior to the jökulhlaups, when we expect the glacier motion to be approximately constant, the velocity records reveal standard deviation of 8, 7 and 11 cm d⁻¹ (3–5 mm h⁻¹) for the Skaftárjökull, 2004 Skeiðarárjökull and 2010 Skeiðarárjökull data, respectively.

The volume of water stored in lake Grímsvötn can be calculated from the mapped hypsometry of the caldera, where the lake is located and the thickness and surface elevation of the lake ice cover (Guðmundsson and others, 1995). The elevation is from a fixed station on the ice cover above the deepest part of the lake. In 2004, the elevation was measured with a navigation GPS (code-based GPS). In 2010, the elevation was calculated from barometric measurements at the station on the ice cover and at a nunatak at ~3 km distance using temperature measurements at the two locations. In both cases, the uncertainty of the daily average elevation is 2–4 m but the derived total elevation drop was 42 m in 2004 and 54 m in 2010. The hypsometry of the caldera is from Guðmundsson and others (1995), updated with post-1996 radio-echo sounding measurements (unpublished data of the Institute of Earth Sciences) to compensate for changes caused by the 1996 Gjálp eruption. The ice-cover thickness is based on radio-echo sounding measurements in 2000.

The amount of additional meltwater produced by friction during the Grímsvötn jökulhlaups can be calculated from the loss of potential energy in the 1280 m descent from the water level in the subglacial lake to the ice margin. We assume that the temperature of the water stored at Grímsvötn is at the pressure-melting point (Ágústsdóttir and Brantley, 1994).

To enable integration of the cumulative jökulhlaup discharge in Gígjukvísl in 2010, we assume that the river discharge increases linearly from a background value at 10:00 on 31 October, when river level starts rising significantly at the bridge, to the first discharge measurement 4 h later. We assume linear change in discharge between discharge measurements, except for the peak discharge. The peak discharge was estimated at ~3000 m³ s⁻¹ at 11:30 on 3 November based on the intersection of exponential approximations to the observed discharge before and after the maximum (unpublished data of the Icelandic Meteorological Office) (Fig. 3). The last discharge measurement was carried out in the afternoon on 4 November, but three additional values of the discharge were derived using the measured river level at the bridge and a linear approximation relating river level and discharge based on simultaneous measurements of the discharge and river level. The river is assumed to have reached background discharge at midnight on 11 November when the river level fully settled.

Fig. 3. (a) Discharge during the 2004 jökulhlaup from Grímsvötn (green), the 2010 jökulhlaup from Grímsvötn (black) and the October 2008 jökulhlaup from the eastern Skaftá cauldron (blue), (unpublished data from the discharge database of the Icelandic Meteorological Office). Discrete discharge measurements are indicated with filled squares and estimated discharge values are indicated with hollow squares. (b, c) Observed elevation changes (b) and horizontal velocities (c) of the glacier surface during the jökulhlaups. The 2004 results (green) are from SKE2, a location 1 km from the glacier margin of Skeiðarárjökull, while the 2010 results (black) are from SKE1, a location 9 km from the margin. Results at two different locations over the flood path of the October 2008 jökulhlaup from the eastern Skaftá cauldron at 3 km, D3 (red), and 8 km, D8 (blue), from the glacier margin are presented.

To estimate water storage in the subglacial flood path of the jökulhlaups, we calculate the difference between the integrated water added into the flood path (flow from Grímsvötn plus meltwater produced by friction) and the outflow of jökulhlaup water from the glacier margin estimated by shifting the Gígjukvísl discharge back 1 h (corresponding to estimated travel time from the glacier margin to the bridge) minus an estimated fixed background discharge of 29 m³ s⁻¹, which leads to equal volume released out of Grímsvötn and discharged into Gígjukvísl. This magnitude of the background discharge can be considered realistic for Gígjukvísl at this time of year.

Analysis of local conduit growth

Estimation of local conduit growth, assuming a steady supply of water from above through a sheet-like flood path formed by a subglacial pressure wave, can be used to investigate possible conduit development during the creation of the flood path at selected locations.

Here, the discharge, Q, along a subglacial conduit is calculated using the Gauckler-Manning-Strickler formula for the mean flow speed v̄, as is traditional in jökulhlaup modelling (e.g. Nye, 1976; Clarke, 2003). Using $Q = S\bar v$ , where S is the cross-sectional area of the conduit, and assuming semicircular channel geometry, the discharge may be expressed as

(1) $$Q = \displaystyle{{S^{4/3} \left \vert {\partial \phi /\partial s} \right \vert ^{1/2}} \over {n_{\rm m} \sqrt {\rho _{\rm w} g} (2/\pi )^{1/3} (2 + \pi )^{2/3}}}. $$

ρ _w is the density of flood water, n _m is Manning's friction coefficient for the conduit walls and bottom and g is the acceleration of gravity. An along-conduit horizontal coordinate is denoted by s, ϕ = ρ _w gz _b + p _w is the hydraulic potential, p _w is the water pressure, z _b is the elevation of the glacier bed and ∂ϕ/∂s is the hydraulic gradient in the conduit.

The time-dependent development of S and Q can be calculated as a function of the hydraulic gradient and an estimate for an initial conduit cross-sectional area, S ₀, derived for example from information about the subglacial discharge of water before the jökulhlaup, using Nye's (1976) theory for jökulhlaups. Creep closure of the conduit during the initial phase of the flood is neglected as the water pressure is assumed to be close to overburden, where this equation is applied. The time-dependent development of S will therefore only be controlled by melting of the conduit walls and is given by

(2) $$\displaystyle{{\partial S} \over {\partial t}} = \displaystyle{{m_{\rm i}} \over {\rho _{\rm i}}}, $$

where t is time, m _i is the melt rate of the channel walls and ρ _i is the density of ice. All initial heat of the flood water in the source lake is assumed to have been dissipated before the flood reaches the location of the GPS instruments, as they are tens of kilometres downstream from the source lake, and heat transfer in subglacial water flow towards the surrounding ice is very rapid (Björnsson, 1992; Jóhannesson, 2002). Due to the rapid heat flow, flood water temperature is assumed to be maintained at the pressure-melting point, so that all released frictional heat is spent locally for melting of ice. Under these assumptions the melt rate of the channel walls will be

(3) $$m_{\rm i} = \displaystyle{Q \over L}\left( {\displaystyle{{ - \partial \phi} \over {\partial s}} + c_{\rm t} \rho _{\rm w} c_{\rm w} \displaystyle{{\partial p_{\rm w}} \over {\partial s}}} \right), $$

where L, c _w and c _t are the latent heat of fusion, the heat capacity and the change of the melting point with pressure for water, respectively. The term in parentheses is the melt-rate ability of the conduit, denoted μ by Magnússon and others (2011). It takes into account the thermodynamic effect, first analysed by Röthlisberger (1972), of the variation of the melting point of water with pressure on the rate of melting of the conduit walls, assuming that energy is released or spent as needed to lower or raise the temperature to maintain it at the local pressure-melting point. Combining Eqns (1–3), the time-dependent development of S can be expressed as

(4) $$\displaystyle{{\partial S} \over {\partial t}} = \displaystyle{{S^{4/3} \left \vert {\partial \phi /\partial s} \right \vert ^{1/2} (( - \partial \phi /\partial s) + c_{\rm t} \rho _{\rm w} c_{\rm w} (\partial p_{\rm w} /\partial s))} \over {\rho _{\rm i} Ln_{\rm m} \sqrt {\rho _{\rm w} g} (2/\pi )^{1/3} (2 + \pi )^{2/3}}}. $$

The local hydraulic potential gradient ∂ϕ/∂s is estimated from the glacier bed and surface topographies, z _b and z _s. Assuming that the subglacial water pressure equals ice overburden as the subglacial hydraulic system is flooded with jökulhlaup water,

(5) $$\displaystyle{{\partial \phi} \over {\partial s}} = (\rho _{\rm w} - \rho _{\rm i} )g\displaystyle{{\partial z_{\rm b}} \over {\partial s}} + \rho _{\rm i} g\displaystyle{{\partial z_{\rm s}} \over {\partial s}}. $$

The bed and surface gradients ∂z _s/∂s and ∂z _b/∂s are estimated over a horizontal distance of one ice thickness at each location.

Analysis of sheet-like water flow and propagation of a subglacial pressure wave

Following Jóhannesson (2002), discharge per unit width, q, in a subglacial sheet, which is much wider than the sheet thickness, h _w, may be calculated from Manning's formula as

(6) $$q = \displaystyle{{h_{\rm w}^{5/3}} \over {2^{2/3} n_{\rm m}}} \sqrt {\displaystyle{{ - \partial \phi _{\rm s} /\partial s} \over {\rho _{\rm w} g}}}. $$

The notation is the same as above for channelized flow, except that the potential ϕ _s includes the sheet thickness h _w and is defined as ϕ _s = ρ _w gz _b + ρ _w gh _w + p _w, where p _w is defined at the top of the sheet. The sheet thickness is considered as an average, h _w = A _w/b _w, over an effective width, b _w, of the flood path with a total cross-sectional area, A _w, resulting in a Manning's hydraulic radius of h _w/2.

The difference, Δp = p _w − p _i, between the water pressure and ice overburden near the front of the subglacial pressure wave drives the uplift of the glacier at the flood front and thereby the propagation of the front. It may be expressed as

(7) $$\Delta p = p_{\rm w} - p_{\rm i} = \phi _{\rm s} - (\rho _{\rm w} - \rho _{\rm i} )gz_{\rm b} - \rho _{\rm i} gz_{\rm s} - \rho _{\rm w} gh_{\rm w}. $$

If the propagating flood front is thicker than the sheet that follows, as hypothesized by Jóhannesson (2002), then the potential gradient, ∂ϕ _s/∂s, within this thicker flow will be smaller than higher upstream. The longitudinal profile of ϕ _s will thus be flatter within the thicker flood front and the water pressure within it will approach a hydrostatic pressure distribution. A propagating front of this type implies a lifting of the glacier surface and a subsequent lowering as the pressure wave passes by, when viewed from a fixed position on the flood path. The water pressure at the base of the glacier can be assumed to be approximately equal to overburden, when lifting stops and the vertical displacement reaches maximum.

If the distance between the front of the flood and the position of maximum vertical displacement can be estimated, as well as the potential gradient within the pressure wave with greater sheet thickness, then Eqn (7) can be used to derive an estimate for the pressure difference, Δp, at the flood front. We assume that the geometry of the propagating pressure wave is approximately constant, so that the depth-averaged velocity of the flow of water is approximately the same within and behind the wave over the time span needed for the wave to pass, and that the discharge behind the pressure wave is driven by the local gradients of the bed and surface topographies. Then the potential gradients within and behind the pressure wave, ∂ϕ _wave/∂s and ∂ϕ _sb/∂s, according to Eqn (6), modified for flow speed, will satisfy

(8) $$\displaystyle{{\partial \phi _{{\rm wave}}} \over {\partial s}} = \left( {\displaystyle{{h_{{\rm sb}}} \over {h_{{\rm wave}}}}} \right)^{4/3} \displaystyle{{\partial \phi _{{\rm sb}}} \over {\partial s}}, $$

and the gradient within the pressure wave may be calculated from estimates of the sheet thickness within and behind the wave, h _wave and h _sb.

Glacier flow is typically dominated by viscous deformation (Cuffey and Paterson, 2010), but the propagation velocity of the flood front and the magnitude of the basal overpressure are such that elastic deformation of the overlying ice may also be important for the separation of the glacier from the bedrock during jökulhlaups. The relative importance of viscous and elastic deformation in a given situation may be assessed by the Maxwell time

(9) $$\tau _{\rm M} = \displaystyle{\eta \over G}, $$

where η is the effective viscosity of the material and G the shear modulus (Melosh, 2011). The effective viscosity of ice (and thereby the Maxwell time) depends on the deviatoric stress, τ _E , as ice is a shear thinning fluid with

(10) $$\eta = \displaystyle{{B^n} \over {2\tau _E^{n - 1}}}, $$

where B and n are parameters in Glen's flow law (Cuffey and Paterson, 2010). Both viscous and elastic deformation are important on timescales close to the the Maxwell time, but elastic deformation dominates on shorter timescales and viscous deformation on longer scales (Melosh, 2011). For typical values of shear and longitudinal stresses and strain rates on temperate glaciers, the Maxwell time is of the order of several hours. As the ice viscosity is reduced by the deformation induced by the propagating pressure wave, Eqn (10) predicts that the Maxwell time will be shorter in the vicinity of the wave than elsewhere.

Jóhannesson (2002) derived an approximate expression for the rate of viscous lifting, w _v, at a flood front of transverse width b _w with water pressure locally higher than the ice overburden

(11) $$w_{\rm v} = \displaystyle{{2h_{\rm i}} \over {n + 2}}\left( {\displaystyle{{b_{\rm w}} \over {2h_{\rm i}}}} \right)^{n + 1} \left( {\displaystyle{{\Delta p} \over B}} \right)^n, $$

where h _i = z _s − z _b is the ice thickness. In case the transverse width is not shorter than the length of the pressure wave as assumed by Jóhannesson, the equation may be assumed to provide an order-of-magnitude estimate for the rate of lifting, if b _w is replaced by the shortest horizontal dimension of the pressure wave over which vertical shear motion associated with the lifting takes place. This equation ignores the effect of strain softening due to the horizontal shearing in the downslope flow of the glacier, which may be estimated to be an order of magnitude smaller than the shear rate induced by the propagation of the flood fronts that are analysed here.

Equation (11) only considers lifting of the glacier by viscous deformation, but elastic deformation may also be expected to play a role in the dynamics of the pressure wave. The magnitude of elastic lifting, ε_e, of the overlying ice after the passage of the initial flood front may be crudely estimated with the theory of elastic plates (Pollard and Fletcher, 2005) as

(12) $${\epsilon} _{\rm e} = \displaystyle{{\Delta p} \over {cR}}a^4. $$

a is the horizontal dimension of the area affected by overpressure Δp acting at the bottom of the plate, $R = E{\kern 1pt} h_{\rm i}^3 /(12(1 - \nu ^2 ))$ is the flexure rigidity of the plate and E and ν are the Young's modulus and the Poisson ratio of the plate material, respectively. If the shortest horizontal dimension of the area, a, is much smaller than the longest horizontal dimension, the constant c = 384, but if the area is circular with diameter a, c = 1024. This equation provides an order-of-magnitude estimate for the elastic response of the glacier if a is the shortest dimension of the subglacial pressure wave.

Elastic lifting according to Eqn (12) approximately captures the effect of elasticity after the passage of the flood front, but not the possible influence of elastic deformation on the propagation of the subglacial front. Subglacial, turbulent, hydraulic fracture propagation in a medium with small fracture toughness has been analysed by Tsai and Rice (2010, 2012) based on linear elastic fracture mechanics (LEFM) for a horizontal fracture that grows parallel to a free surface, both for a short fracture and for a fracture that eventually becomes long compared with the ice thickness. They take into account the time needed for the water flow to fill the fracture as it propagates down-glacier and grows in thickness utilizing the Gauckler–Manning–Strickler approximation to represent turbulent friction in the fracture. The analysis is intended to shed light on the emptying of supraglacial lakes, and does not fully correspond to the situation analysed here, as Tsai and Rice assume maximum effective pressure at a fixed inflow location. This leads to a steadily growing separation between ice and bedrock all the way between the inflow location and tip of the propagating fracture at the ice/bedrock interface. The analysis of Tsai and Rice (2012) may be expected to approximately represent the tip of a propagating jökulhlaup, if the inflow location is assumed to correspond to a moving position at a distance L _Δp behind the tip although their solution regarding the variation of pressure and discharge along the flood front will not be fully consistent in this case. This leads to the following expression (in our notation) for propagation velocity of the crack tip:

(13) $$U_{{\rm tip}} = \theta \sqrt {\displaystyle{{\Delta p} \over {\rho _{\rm w}}}} \left( {\displaystyle{{\Delta p} \over {E^{\prime}}}} \right)^{2/3} \left( {\displaystyle{{L_{\Delta {\rm p}}} \over k}} \right)^{1/6}, $$

where θ is a non-dimensional multiplicative factor (denoted ϕ by Tsai and Rice), E′ = E/(1 − ν ²) and k is the Nikuradse roughness height. Tsai and Rice (2012) give an approximation for the factor θ ≈ 5.13 + 0.64(L _Δp/h _i) + 0.94(L _Δp/h _i)², valid for L _Δp/h _i ≤ 5, and the Nikuradse roughness height is related to the Manning's friction coefficient by $n_{\rm m} = (0.0380\;{\rm s}{\kern 1pt} {\rm m}^{ - 1/2} )\;k^{1/6} $ .

The analysis of Tsai and Rice also yields approximate solutions for the crack opening (again in our notation)

(14) $$h_{\rm w} = \displaystyle{{\hat w\Delta pL_{\Delta {\rm p}}} \over {E^{\prime}}}. $$

The time-varying spatial average of the non-dimensional crack opening, $\hat w$ , is approximately given as a function of the ratio of the crack length to the ice thickness $\hat w_{{\rm avg}} \approx 1.72 + 0.89(L_{\Delta {\rm p}} /h_{\rm i} )^2 $ , and the ratio of the maximum opening to the average opening may be read off figure 5 in Tsai and Rice (2012) for 0.02 ≤ L _Δp/h _i ≤ 5.

The vertical uplift rates predicted by viscous deformation (Eqn (11)), the lifting predicted by plate elasticity (Eqn (12)) and the crack opening given by the LEFM analysis of Tsai and Rice (Eqn (14)) are not directly comparable with our jökulhlaup observations, as the lifting shown by the GPS instruments will be the result of both elastic and viscous deformation. We will calculate the contributions of the different processes described above and compare them to our observations to investigate, which processes are likely to be important for the propagation of the subglacial pressure wave.

RESULTS

Observations of jökulhlaups from the Skaftá cauldrons

The October 2008 jökulhlaup in Skaftá

The October 2008 jökulhlaup from the eastern Skaftá cauldron was among the larger jökulhlaups in Skaftá, reaching maximum discharge of flood water at Sveinstindur of 1290 m³ s⁻¹ (Fig. 3). The total volume of flood water was 0.27 km³. A delay of 3 h is estimated for the propagation of the flood from the glacier margin to Sveinstindur, based on a simple calculation for wave propagation in open channels. The discharge at the glacier margin started rising at 02:45 on 11 October, reached a maximum ~26 h later at 05:00 on 12 October, and was back to expected base flow ~19 October.

On-site observations showed flow from large continuous segments of the ~5 km wide glacier margin of the Skaftárjökull outlet glacier during the beginning of the jökulhlaup. Around midday on 11 October, substantial flow was still observed from an overflight from large continuous parts of the glacier margin, but a few main outlets discharged most of the water. During a field trip late in the evening of 11 October, the jökulhlaup was observed to discharge from a few main outlets.

Indications of high water pressure near the terminus during the initiation of the jökulhlaup, such as water flowing out of moulins and a ~2 km long, transverse thrust fault formed through the glacier ~1 km from the ice margin, were observed from an overflight and during a field trip. The formation of the fault, which was ~2 m high, was associated with an M2 glacier originated earthquake at 02:05 on 11 October located near the glacier terminus by Iceland's seismic network (unpublished data of the Icelandic Meteorological Office). Breaking-up of the glacier snout was observed as well.

The arrival of the jökulhlaup at GPS station D8 at 22:30 on 10 October was associated with lifting of 1.06 m in 50 min followed by lowering of 0.41 m in 50 min (Fig. 4). A second phase of slower lifting then started, with total vertical displacement of 0.58 m in 11 h. The vertical position of the station then remained constant for ~14 h at a total displacement of 1.23 m, compared with the initial elevation before the jökulhlaup. The glacier surface started to fall at 01:30 on 12 October and had been lowered by 0.77 m when the station lost power on 14 October.

Fig. 4. Observed horizontal displacements (a, b) and elevation changes (c) during the October 2008 jökulhlaup from the eastern Skaftá cauldron at 3 km, D3 (red), and 8 km, D8 (blue), from the glacier margin. The onset of discharge increase at the glacier margin is indicated with a dashed vertical line and the timing of an M2 earthquake near the terminus with a dotted line.

The ice movement at station D3 at the beginning of the flood was similar to the movement at D8 with a delay of ~2.1–2.5 h, with lifting of 0.84 m in 40 min followed by lowering of 0.47 m in 60 min (Fig. 4). A second phase of slower lifting then started with total lifting of 0.12 m in 2 h followed by slow lowering over the next 5 d by 1.25 m in total, after which no further changes in elevation were observed. The station was thus 0.76 m lower after the jökulhlaup than before.

Substantial changes in horizontal velocity were observed at both D8 and D3 during the jökulhlaup (Fig. 3). The onset of horizontal motion was more gradual and started some tens of minutes before the vertical motion (Fig. 4). Rapid horizontal movements coincided with rapid vertical movements, and the horizontal velocity generally increased during lifting and started to decrease as soon as lifting stopped. Reversed loop motion (in a vertical section) was superimposed on the above-described general displacement pattern at both stations. The reverse motion occurs around the time of maximum uplift during the passage of the subglacial wave at the onset of the jökulhlaup at each location. The timing and relative magnitude of changes in the south–north and west–east velocity components often did not coincide during the jökulhlaup, causing temporary changes in the ice-flow direction by several tens of degrees.

Lowering due to surface melting is estimated at 2–3 mm d⁻¹ by simple degree-day calculations, both at D8 and D3, and can be ignored compared with the observed vertical motions. Lowering due to the movement of the instruments down the surface gradient of the glacier, with the horizontal flow, is 2–4 mm d⁻¹ at both stations for days before and after the jökulhlaup. During the jökulhlaup, bed slip is substantial, leading to a combination of normal ice flow following the surface gradient and block flow controlled by the bed gradient. Using the higher of the two gradients leads to estimated lowering of 25 mm d⁻¹ at D8 and 45 mm d⁻¹ at D3 at the time of maximum horizontal velocities during the jökulhlaup. This is one to three orders of magnitude smaller than the observed vertical velocities and is also ignored here.

The propagation velocity of the subglacial flood front is estimated as 0.5–0.6 m s⁻¹ between D8 and D3 and 0.4–0.5 m s⁻¹ between D3 and the glacier margin, based on the timing of the onset of lifting and the time of the beginning of discharge increase at the glacier margin. The propagation velocities allow an estimation of the length of the subglacial pressure wave (Jóhannesson, 2002), which is assumed here to cause the first phase of lifting and subsequent lowering, to be 3.0–3.5 km at D8 and 2.3–3.5 km at D3.

The overpressure, Δp, at the flood front may be calculated from Eqns (7) and (8). The water pressure at the base of the glacier is assumed to be equal to overburden when lifting stops and the vertical displacement reaches maximum at the distance L _Δp from the flood front. L _Δp is found to be 1.8–1.5 km at D8 and 0.9–1.3 km at D3. This leads to Δp = 0.16–0.18 MPa at D8 and Δp = 0.17–0.25 MPa at D3. The correction for the effect of friction in the water flow within the flood front expressed by Eqn (8) amounts to ~50% of the estimated overpressure.

The August 2008 jökulhlaup in Skaftá

The August 2008 jökulhlaup from the western Skaftá cauldron was a relatively small event, with a maximum discharge of flood water at Sveinstindur of 240 m³ s⁻¹ (Fig. 5) and total volume of flood water of 0.10 km³. A delay of 3.5 h is in this case estimated for the flood propagation from the glacier margin to Sveinstindur. The start of the flood at the glacier margin was obscured by the diurnal discharge variation in the base flow from the glacier. Conductivity started to rise at ~19:00 on 8 August, indicating some flood water from the geothermal source lake at the margin. Notable discharge increase was first observed around noon on 9 August and the discharge increased sharply until the evening of 10 August. The maximum in discharge, on 11 August, might have been caused by diurnal variation in base flow rather than a maximum in flood water discharge. The discharge receded back to expected base flow ~17–18 August. During a field trip on 10 August, the jökulhlaup was observed to flow out of a few main outlets and no evidence of a widespread outflow from the glacier margin or flow out of fractures through the glacier surface was observed.

Fig. 5. (a) Discharge during the August 2008 jökulhlaup from the western Skaftá cauldron (unpublished data from the database of the Icelandic Meteorological Office). (b, c) Observed elevation changes (b) and horizontal velocities (c) of the glacier surface at 3 km, D3 (red), 8 km, D8 (blue), and 15 km, D15 (black), from the glacier margin during the jökulhlaup.

Slow lifting started at D15 ~00:30 on 9 August (Fig. 5). About 3 h later, the rate of lifting suddenly increased by a factor of 4 and rapid lifting sustained for 2 h. The rate of lifting then decreased substantially and slow lifting took place for the next 11 h, reaching total vertical displacement of 0.15 m. A flat maximum in the displacement was then maintained for ~20 h. Slow lifting started at D8 ~07:30 on 9 August. About 7 h later, the rate of lifting was accelerated by a factor of 15–20 and rapid lifting took place for 1 h. This was followed by slower lifting for 3 h and maximum vertical displacement of 0.35 m was reached at 19:00. After the maximum was reached, slow lowering took place at both locations until ~15 August when the stations reached the rate of lowering expected from melting and down-glacier ice movement. The total vertical displacement of the instruments over the course of the jökulhlaup was equal to the expected lowering due to melting and downslope movement of the GPS instruments with the ice, within +/−5 cm and ±10 cm for D15 and D8, respectively.

Substantial changes in horizontal velocities were observed at both D15 and D8 during the jökulhlaup. The horizontal velocity generally increased during lifting and started to decrease as soon as lifting stopped. A short event of reversed motion in the south–north direction took place at D8 around the time of maximum vertical displacement. The timing and relative magnitude of changes in south–north and west–east velocity components did not coincide during the jökulhlaup, as for the October event, causing temporary changes in the ice-flow direction by tens of degrees.

No changes were observed at D3 during the jökulhlaup (Fig. 5). The vertical movement was as expected, due to surface melting and down-glacier ice movement. The horizontal velocities were approximately constant.

As for the October 2008 jökulhlaup, the effect of surface melting and down-glacier movement with the ice of the GPS instruments at D15 and D8 on the measured vertical velocities is small in comparison with the velocity changes induced by the jökulhlaups and can be ignored.

The propagation velocity of the subglacial flood front is estimated at ~0.2 m s⁻¹ between D15 and D8. The propagation velocity cannot be calculated between D8 and D3 as no signal of the flood front was observed at D3. The propagation velocity between D8 and the glacier snout is uncertain, as the start of flood at the ice margin is unclear, but it can be inferred to have been in the range 0.1–0.3 m s⁻¹.

The length of the subglacial pressure wave is harder to estimate than for the October event because the lowering in the wake of the pressure wave is hard to distinguish. The distance from the flood front to the point of maximum vertical displacement can be estimated at L _Δp = 1–2 km for D15 and L _Δp = 2–4 km at D8. The potential gradient that drives the sheet flow inside the pressure wave cannot be calculated from Eqn (8) because the sheet thickness behind the wave, h _sb, cannot be estimated. Assuming that the correction for head loss within the pressure wave is 50% of the estimated overpressure, similar to that for the October event, leads to Δp = 0.05–0.10 MPa at D15 and Δp = 0.15–0.31 MPa at D8.

Interpretation of observations of the Skaftá jökulhlaups

Observations of outflow from large continuous segments of the glacier margin during the beginning of the October 2008 jökulhlaup in Skaftá suggest drainage from a subglacial sheet-like flow in the wake of a pressure wave and lifting due to basal water pressure higher than ice overburden (Sugiyama and others, 2008). Such sheet-like flow is unstable to irregularities in sheet thickness, so conduits are expected to eventually grow by melt–discharge feedback (Walder, 1982) in areas where the melt-rate ability of the conduit, μ, is positive (Magnússon and others, 2011). The development of the subglacial flood path in the wake of the pressure wave is consequently dependent on the local value of μ. Contraction of the flood path into a few main conduits near the glacier margin is expected as μ is positive there and up to the location of D3.

The second phase of lifting at D3 and subsequent lowering, while the discharge of the October jökulhlaup was still increasing (Fig. 3), can be interpreted as an evolution of an initial sheet-like flood path into efficient conduits. Increasing discharge for a sheet-like subglacial flood path implies a thicker sheet and lifting of the glacier surface, whereas conduit growth can increase the efficiency of the flood path without increased bed separation. Draining of the water stored in the sheet into the growing conduits then leads to the observed lowering. Lower surface elevation at D3 at the end of the jökulhlaup may be due to the closure of such a conduit or erosion of bed sediments by the flood (Ng, 1998) but this cannot be ascertained on the basis of the available observations.

Sustained lifting while local inflow was increasing at D8 during the October jökulhlaup and at D15 during the August jökulhlaup implies increased flood path efficiency by thickening of the sheet during the whole event with little or no conduit growth. This is consistent with the adverse bed slope at D15, i.e. more than four times the magnitude of the surface slope (Fig. 2), which may be expected to cause glaciohydraulic supercooling at this location (Alley and others, 1998; Clarke, 2005). This leads to negative melt-rate ability μ that prevents melt enlargement of conduits. Adverse bed slope a short distance downstream of D8 may be expected to inhibit conduit growth there too as further discussed below.

A second phase of uplift was not detected at D8 during the August 2008 Skaftárjökull jökulhlaup, which had an order-of-magnitude smaller discharge and much smaller induced surface motions compared with the October flood. The reason for this difference in induced ice motion at D8 between the August and October events is not obvious but could be related to transverse variability across the subglacial path in the much smaller August flood that is not resolved by our point observations.

The difference between the ice motion at D3 and D8 induced by the October jökulhlaup and the difference between ice motion at D8 in the October and the August jökulhlaups demonstrates the importance of local surface and bed topographies for the path development and shows that detailed temporal and spatial modelling of the flood path development and ice dynamics is needed to understand rapidly rising jökulhlaups.

The subglacial propagation velocities for the two jökulhlaups in Skaftá are in the range 0.1–0.6 m s⁻¹, similar to ~0.5 m s⁻¹ estimated by Magnússon and others (2007) for a small jökulhlaup in Skaftá in 1995 and 0.2–0.4 m s⁻¹ for a jökulhlaup from the western cauldron in 2006 (unpublished data from the Icelandic Meteorological Office). These speeds are substantially lower than the mean propagation velocity of ~1.3 m s⁻¹ for the large jökulhlaup from Grímsvötn in November 1996 (Björnsson, 2002), indicating that large subglacial floods propagate faster along the glacier bed as would be expected from the larger spatial dimensions of the subglacial flood tongue and higher overpressure of the larger floods.

Dynamics of the flood front at Skaftárjökull

The uplift associated with the passage of a subglacial pressure wave can be caused by horizontal compression of the glacier due to longitudinal velocity gradients across the front as well as by viscous and elastic lifting of the glacier from the bed. The contribution due to horizontal compression for the two jökulhlaups can be estimated to be <5–15% based on the thickness of the glacier, L _Δp and horizontal motion accumulated between the arrival of the wave and maximum horizontal velocity and is therefore ignored here.

The Maxwell time for the glacier before the arrival of the flood front is of the order of several hours, which is longer than the observed 30–120 min duration of lifting associated with the passage of the flood front. The Maxwell time for the deforming ice above the pressure wave is considerably shorter. Using Eqns (9) and (10) and physical parameters from Table 1, it is found to be several minutes for overpressures of 0.15–0.35 MPa at D3 and D8 in both jökulhlaups and <2 h for an overpressure of 0.05–0.10 MPa, which was estimated at D15 during the August jökulhlaup. The Maxwell time thus indicates that the flood front propagates in glacier ice where elastic effects are initially dominant, before the deformation pattern of the flood front is established, but strain softening induced by the separation of the glacier from the bed eventually leads to dynamics that are dominated by viscous deformation.

Table 1. Physical parameters used in conduit growth calculations

The vertical ice velocity due to viscous deformation and the elastic lifting in connection with the passage of the flood front may be calculated with Eqns (11) and (12) from estimates of the overpressure and the shortest horizontal dimension of the pressure wave. In a similar manner, the propagation velocity of the flood front and the maximum crack opening predicted by the LEFM theory of Tsai and Rice (2010, 2012) can be found with Eqns (13) and (14). Because of the approximate nature of these analyses, the results is expected to be interpreted as order-of-magnitude estimates, intended to give an indication of the plausibility of the pressure wave as a conceptual model for the propagating subglacial flood front. The calculated vertical velocities and elastic lifting fall within wide ranges, w _v = 0.03–55 m d⁻¹ and ε_e = 0.05−0.3 m, for overpressures of 0.15–0.35 MPa acting over horizontal dimensions in the range 1–2 km, and using parameters from Table 1. The viscous vertical velocities correspond to total viscous lifting of 0.01–2 m during the passage of the wave. The lowest estimated value of 0.05 MPa for the overpressure at D15 in August leads to much lower velocities that cannot be considered realistic. More reasonable velocities can be obtained if a longer length scale is chosen as Eqn (11) is very sensitive to the choice of the length scale. The observed rates of vertical lifting, w _v = 1−30 m d⁻¹, and maximum lifting of 0.15–0.98 m, are on the same orders of magnitude but the wide ranges preclude detailed comparison for each location in the two jökulhlaups. The analysis indicates that the pressure wave has a lateral width greater than or equal to the along-flood-path length because a narrower wave, of the order of several hundred metres wide for example, would lead to much lower vertical velocities (of the order of mm d⁻¹) and maximum elastic lifting (of the order of mm) that cannot be reconciled with the observation of the propagating flood front.

The LEFM theory of Tsai and Rice (2012) leads to crack-tip propagating velocities U _tip = 0.02−9 m s⁻¹ for overpressures Δp = 0.05−0.35 MPa, L _Δp = 1−2 km, (L _Δp/h _i) = 2−5, Nikuradse roughness heights calculated from Manning's friction factor, n _m = 0.01−0.1 m^−1/3 s as reported for subglacial conditions by Cuffey and Paterson (2010), and physical parameters from Table 1. This range overlaps the observed flood front velocities of 0.1–0.6 m s⁻¹ for the jökulhlaups. The predicted elastic lifting at the point of maximum overpressure is found to be 0.02–3.6 m, again overlapping the observed lifting of 0.15–0.98 m. A narrower range of overpressures, 0.15–0.35 MPa and a middle-of-the-range value, n _m = 0.033 m^−1/3 s, for the Manning roughness, corresponding to natural river beds, leads to crack-tip velocities in the range 0.4–3 m s⁻¹ and maximum lifting of 0.08–3.6 m. These calculations indicate that turbulent, hydraulic fracture at the tip of the pressure wave is an important physical process in the propagation of the flood front and that both elastic and viscous deformations contribute to the lifting of the glacier.

Local conduit growth at Skaftárjökull

The field observations of the October 2008 jökulhlaup show development from sheet flow to conduits at the glacier margin over the course of the jökulhlaup, and the GPS measurements at D3 during the August jökulhlaup indicate that the subglacial water flow at D3 was confined to conduits for the whole jökulhlaup. Farther up-glacier, this jökulhlaup appears to have started with the propagation of a pressure wave implying sheet flow. Conditions for conduit development can be analysed in terms of the melt rate ability, μ = −(∂ϕ/∂s) + c _t ρ _w c _w(∂p _w/∂s), using Eqn (4) to estimate conduit enlargement based on the local flood path geometry, an estimate for an initial conduit size and physical parameters from Table 1. The modelled conduit development and estimated initial cross-sectional area is highly dependent on the selected Manning's coefficient for friction, n _m. We will consider the same range of n _m as above but use n _m = 0.033 m^−1/3 s, corresponding to natural river beds, when we need a single explicit value.

The melt-rate ability is negative at D15 and a short distance downstream from D8 so conduit growth is not expected at those locations as further discussed below. The melt-rate ability is positive at D3, so conduits are predicted to grow there. An initial cross section of 6–30 m² for a conduit at D3 in August 2008 is estimated for the range of n _m given above, if a single conduit that carries the measured pre-jökulhlaup river discharge is assumed (after subtraction of runoff from the watershed below this location (Kristinsson, 2005)). The increased carrying capacity of such a conduit in 2 d is 100–1000 m³ s⁻¹, and 720 m³ s⁻¹ for n _m = 0.033 m^−1/3 s. The observed discharge increase during the first 2 d of the jökulhlaup is 170 m³ s⁻¹. So this analysis shows that conduit growth starting from the late-summer subglacial hydrological conditions at D3 can accommodate the observed rise in the discharge of the August 2008 jökulhlaup, which is consistent with the GPS measurements at D3 that show no lifting of the glacier for this jökulhlaup.

Estimating a plausible initial cross-sectional area in August is harder at D8 than D3. Assuming an initial cross-sectional area of 10 m² as an upper bound results in conduit growth that is an order of magnitude too slow to accommodate the observed discharge increase. The initial rise in the discharge at this location is, thus, hard to explain without invoking sheet flow, which is consistent with the GPS measurements that show lifting of the glacier indicating sheet flow at the beginning of the jökulhlaup.

Similar calculations using the measured base discharge before the October 2008 jökulhlaup to estimate initial cross-sectional areas for a conduit at D3 and D8 indicate that conduit growth is two orders of magnitude too slow to accommodate the rapid initial discharge increase of the October flood.

Observations of jökulhlaups from Grímsvötn

The two jökulhlaups from Grímsvötn are similar in terms of duration, discharge and total volume (Fig. 3). The lowering of the Grímsvötn ice cover corresponds to the release of ~0.5 and 0.6 km³ during the 2004 and 2010 jökulhlaups, respectively. An additional ~0.1 km³ of meltwater was produced during the 2004 eruption in Grímsvötn (Guðmundsson and others, 2009), starting around the peak of the jökulhlaup. Up to ~0.2 km³ or ~1/3 of the flood water was stored in the flood path beneath the glacier (Fig. 6) during the 2010 jökulhlaup. Despite significant uncertainty in the ice-cover elevation at any given time, the peak subglacial water storage is fairly well established (uncertainty <25%) due to the time lag between the outflow from Grímsvötn and the observed river discharge. Analysis of the 2004 jökulhlaup (not shown) similar to that used to derive Figure 6 for the 2010 jökulhlaup shows that when the eruption started in 2004, ~0.2 km³ of water was already stored in the subglacial flood path. This flood water storage is an order of magnitude larger than the volume of a conduit formed by melting during these floods.

Fig. 6. Observed accumulated drainage out of Grímsvötn (blue) and accumulated discharge at the glacier margin in Gígjukvísl (black) during the 2010 jökulhlaup from Grímsvötn. Subglacial storage (green) is estimated as the difference of the two curves plus the calculated amount of melt due to friction in the flow (magenta). Estimated uplift due to water accumulation (red) at SKE1 is presented by the right y-axis.

Seismic tremor measured at Grímsfjall, ~1 km from the lake, and elevation measurements on the Grímsvötn ice cover suggest that water started to flow out of the lake no later than 26 and 27 October in 2004 and 2010, respectively. The influence of the flood water on glacier motion at the GPS station was in both cases observed ~3.5 d later, but increased conductivity was measured in Gígjukvísl on 29 October during the 2010 jökulhlaup, a day before acceleration of surface motion was detected. Timing of the onset of increased discharge at the glacier margin is uncertain in the two jökulhlaups as only water stage is measured in Gígjukvísl (in 2010) and Skeiðará (in 2004) and the calculation of discharge from water stage for braided glacial rivers on a sediment bed is uncertain. Clear signs of increased discharge are seen shortly after 01:30 on 30 October in 2004 and 10:00 on 31 October in 2010 but the discharge may have started to increase up to a day earlier in both cases. The delay due to the flow from the glacier margin to the water stage and discharge measurement sites may be estimated to be ~1 h.

The horizontal and vertical displacement records at the GPS stations do not indicate the propagation of a pressure wave at the beginning of the jökulhlaups as is observed on Skaftárjökull in August and October 2008. The effect on the glacier motion is quite different between the stations (Fig. 3). In 2010, gradual acceleration in horizontal motion probably starts ~18:00 on 30 October at SKE1 located 9 km from the glacier margin, where the glacier is ~600 m thick (Fig. 2). The velocity increases from 0.4 m d⁻¹ between 17:00 and 18:00 to 2.1 m d⁻¹ around midnight and continues to accelerate to a peak velocity of 4.0 m d⁻¹ between 05:00 and 06:00 the following morning. Uplift of the station starts shortly before midnight, ~6 h later than the horizontal acceleration. The station rises by 0.4 m until 06:00 the following morning, when the rate of rise slows down by half. The vertical rise remains constant at ~0.8 m d⁻¹ and the horizontal velocity ~3 m d⁻¹ for 3 d until 06:00 on 3 November. The rise then stopped and the station remained at approximately constant elevation for 6 h, while the peak discharge of the jökulhlaup passed the station. While the GPS station remained at this elevation, the horizontal velocity decreased from 3.1 to 1.1 m d⁻¹. The station subsided for ~1 week after the discharge peaked, by 0.9 m during the first 24 h, but then at a gradually slower rate and the horizontal velocity decreased gradually from 1.1 m d⁻¹ to a typical winter velocity value of 0.3 m d⁻¹.

In the morning of 30 October 2004, strong acceleration from ~0.15 m d⁻¹ between 06:00 and 07:00 to ~1.0 m d⁻¹ between 07:00 and 08:00 was observed at the GPS station SKE2, 800 m from the glacier margin where the ice thickness is ~200 m (Fig. 2). The start of the rise at the station is slow. The timing of the start is therefore inaccurate but seems to be ~3 h later. The station rises ~0.20 m during the next 2 d and the horizontal velocity fluctuates ~1 m d⁻¹. Around 17:00 on 1 November, the station starts accelerating steadily from 1.1 m d⁻¹ between 17:00 and 18:00 to 2.2 m d⁻¹ between 22:00 and 23:00. The rise also becomes faster ~22:00; the station rose by 0.15 m during the following 12 h. An eruption started in Grímsvötn in the evening of 1 November. A seismic tremor measured at Grímsfjall indicates the onset of the eruption at 21:50, a few hours after the horizontal velocity at SKE2 started increasing and at approximately the same time as the uplift rate increased. The horizontal velocity decreases gradually from 2.2 m d⁻¹ at midnight to 1.4 m d⁻¹ between 12:00 and 13:00 the following day, around the time when the maximum vertical displacement has been reached. Then the horizontal velocity was reduced by half to 0.7 m d⁻¹ in 1 h after which the glacier gradually decelerated. The station stayed at the peak elevation for ~8 h, while the peak discharge of the jökulhlaup passed this location, and then started subsiding slowly. When the station was taken down 4 d later, it had subsided by 0.11 m and the subsidence was still ongoing. The daily horizontal displacement of the station at that time was 0.15 m, whereas the daily displacement prior to the jökulhlaup was ~0.08 m.

The total vertical displacement at both stations on Skeiðarárjökull is positive, most likely due to adverse bed slope below the stations (Fig. 2). During rapid basal sliding, the glacier may be assumed to move effectively as a solid block along the bed, causing the bed slope to control the direction of ice motion at the surface. It is difficult to estimate the bed slope underneath the station accurately based on the existing bedrock data (Björnsson and others, 1999; Magnússon, 2008; Magnússon and others, 2009). The slope of the 3-D total displacement vector during the jökulhlaups, 2.6° at SKE1 in 2010 and 3.6^° at SKE2 in 2004, is within the likely range of adverse bed slope at both locations. If we compensate for vertical motion due to block movement along the bed using the above slope values, the residual vertical displacement may be assumed to be caused by water accumulation and depletion. It reaches a maximum of 2.5 m at SKE1 in 2010 (Fig. 6) and 0.17 cm at SKE2 in 2004. Uplift due to water accumulation may be somewhat underestimated at SKE2 because the station was still subsiding when it was taken down.

As mentioned above, the damage done to the seal at Grímsvötn by the catastrophic jökulhlaup in November 1996 altered the characteristics of slowly rising jökulhlaups from the subglacial lake. Figure 7 compares the rise of the 2004 and 2010 jökulhlaups studied here with seven earlier slowly rising jökulhlaups for which discharge measurements are available (Rist, 1955, 1973; unpublished data from the Icelandic Meteorological Office). The timescale needed for the flood discharge to rise by a factor of e = 2.7 was shortened from 3–6 d to slightly over a day in 1996. This distinct change indicates that the timescale of the rise in discharge is controlled by conditions at a particular location of the flood path rather than the large-scale path geometry, such as the vertical relief, the total path length or the overall scale for the ice thickness or potential gradient, as these quantities did not change fundamentally in 1996.

Fig. 7. Discharge variation with time for nine slowly rising jökulhlaups from Grímsvötn. The dashed lines show a log–linear fit, Q = ae ^bt , to the rising limb of each curve, where t is time in days since the start of the rise in discharge and a and b are statistical parameters. The slope of the least-squares line, b (d⁻¹), for each jökulhlaup is given in the legend. The timescale needed for the exponential approximation to the discharge to rise by a factor of e = 2.7 is given by 1/b. It was shortened by a factor of ~3–4 after the catastrophic jökulhlaup in 1996.

Interpretation of observations of the Grímsvötn jökulhlaups

The lifting at SKE1 on Skeiðarárjökull during the 2010 jökulhlaup shares a number of characteristics with the second phase of lifting at D8 on Skaftárjökull during the October 2008 jökulhlaup, such as close resemblance to the variation in discharge (see Figs 3 and 6). The geometry of the flood path at the two locations is similar, with an adverse bed slope, steep enough to prevent the development of conduits, less than one ice thickness farther down-glacier (Fig. 2). We interpret the rise in the ice surface as an indication of widespread subglacial accumulation of water that we expect to occur upstream of locations with adverse bed slope and conditions for supercooling that prevent the formation of efficient conduits. This accumulation is driven by two processes. Firstly, the rising discharge carried by sheet flow over the adverse bed slope requires a correspondingly greater sheet thickness, h _w, that may be expected to reach some distance up-glacier from the start of the back slope, due to rigidity of the ice on timescales of hours, further enhanced by block flow of the glacier up the adverse bed slope. Secondly, the smaller efficiency of the path without conduits on the adverse bed slope will lead to greater hydraulic gradient there, which leads to comparatively high subglacial water pressures, lifting and accumulation of water upstream (Magnússon and others, 2011).

Uplift and subsequent subsidence over the duration of the jökulhlaup is also observed at SKE2 in 2004, probably due to water accumulation upstream of adverse bed slope, but an order of magnitude smaller than at SKE1 in 2010 (Fig. 3). The reason for this difference in magnitude is not clear but might be related to the short distance (800 m) to the glacier margin where the water pressure becomes equal to atmospheric pressure and remnants of conduits from the previous summer may still be present due to slow creep closure under thin ice.

Subglacial accumulation of water of the magnitude encountered at both Skaftárjökull and Skeiðarárjökull during the jökulhlaups considered here, can be expected to greatly reduce basal drag over large areas by drowning bed roughness elements and decreasing effective pressure at the bed (Iken and Bindschadler, 1986). This change in the force balance of the glacier, from local control to global control (Cuffey and Paterson, 2010), is a probable explanation for the observed changes in horizontal velocities and ice-flow direction during the jökulhlaups. Similar increase in horizontal velocity and changes in flow direction in connection with enhanced basal slip during jökulhlaups or intense rainfall events have been observed before for Skeiðarárjökull (Magnússon and others, 2007, 2011) and observed and modelled for Gornergletscher, Switzerland (Sugiyama and others, 2007; Riesen and others, 2010).

DISCUSSION

The GPS measurements of ice motion during the four jökulhlaups described above, show many features that shed light on the dynamics of the subglacial hydrological system during rapid changes in water input. The dynamics are apparently complex and depend among other things on the state of the subglacial hydrological system before the flood and the local geometry of the glacier surface and the bed.

Grímsvötn is well known in the glaciological literature, and the jökulhlaup in 1972 (Rist, 1973) was used by Nye (1976) in his development of a dynamical theory for slowly rising jökulhlaups where the conduit-melt–discharge feedback was first used as an explanation for the rise in jökulhlaup discharge with time. In its original form, this theory assumes that the rise of the flood is driven by the average potential gradient from the source lake to the outlet. The change in the character of slowly rising jökulhlaups from Grímsvötn in 1996, when the timescale of the rise in discharge was shortened by a factor of ~3–4 (Fig. 7), indicates that the dynamics of slowly rising jökulhlaups is more complex than can be reconciled with Nye's theory and points to local control or bottlenecks in the path (Clarke, 1982, 2003) that determine the rise in discharge. It appears that the damage done to the seal in 1996 moved the control of the flood to a different location where conditions allowed a much more rapid rise in discharge. This implies that Manning's friction coefficient, determined from Nye's theory by statistical fitting to the variation of the discharge with time using a single along-path average of the potential gradient, must be considered as a phenomenological parameter that does not have a clear meaning in terms of conduit roughness or other physical properties of the path. The same applies to similar theories based on along-path averages of the potential gradient and conduit cross section (e.g. Clarke, 1982).

The lifting of the glacier surface shown by the GPS measurements and the subglacial accumulation of water in the 2004 and 2010 jökulhlaups reveal another interesting aspect of slowly rising jökulhlaups that cannot be accounted for by Nye's theory or its later developments (e.g. Spring and Hutter, 1981, 1982; Clarke, 2003). Subglacial water accumulation at Skeiðarárjökull has also been observed before the change of the timescale of discharge rise in 1996, for a jökulhlaup under normal conditions (Magnússon and others, 2007) and for a jökulhlaup during a surge (Björnsson, 1998). The subglacial accumulation of 0.2 km³ of flood water in the 2010 jökulhlaup corresponds to a 2 m thick and 2 km wide tongue of water extending along the entire 50 km long flood path. Of course, the water stored subglacially does not have this simple geometry. Most likely, it was accumulated in a few areas along the flood path of greater width and water depth where subglacial ponding is favoured by the local glacier surface and bed topographies. The discharge from glacier margin during the jökulhlaup is a combination of outflow from the source lake, with a delay because of the time needed for the water to flow down-glacier, and the release of ponded flood water stored subglacially. In spite of the fact that Nye's theory does not account for this subglacial ponding, the rise in discharge can be reproduced surprisingly well by the theory with parameters statistically determined to fit the observed discharge curve as mentioned above. This points to local control or bottlenecks in the path that determine the rise in discharge by conduit-melt–discharge feedback. These bottlenecks would not necessarily be located at a seal below the thickest ice as envisioned by Clarke (1982). They could be located in regions of channel growth between areas of subglacial storage and glacier lifting. Melting at such bottlenecks could be driven by a potential gradient that differs substantially from both the average potential gradient for the whole flood path and the gradient predicted by modelling assuming conduit flow along the path (Spring and Hutter, 1981, 1982; Clarke, 2003). This again indicates that such theories should be considered phenomenological rather than physical so that the parameters, such as the Manning's coefficient, n _m, and the variables, such as the channel cross section, S, may not have a clear physical meaning.

Jökulhlaups in Skaftá appear to be fundamentally different from jökulhlaups from Grímsvötn in having a much more rapid rise in discharge that does not have the near-exponential form so characteristic of most floods from Grímsvötn (Fig. 7). The GPS measurements show that the initiation of these rapidly rising floods is associated with a wave of surface lifting and subsequent subsidence moving down-glacier, in agreement with the concept of a localized, subglacial pressure wave suggested by Jóhannesson (2002). A several-km wide pressure wave followed by sheet flow is, furthermore, consistent with observations that indicate reduced basal friction over a several-km wide flood path for a jökulhlaup at Skaftárjökull in 1995 (Magnússon and others, 2007). Order-of-magnitude agreement between theoretical predictions and the observed ice velocities, lifting of the glacier and the down-glacier propagating velocities of the jökulhlaup wave indicates that viscous and elastic deformation as well as turbulent hydraulic fracture at the tip of the flood front are all important processes in the propagation of the pressure wave.

We have not identified quantitative conditions that determine whether a jökulhlaup develops as a rapidly rising flood, with an initial subglacial pressure wave, or slowly, with a near-exponential rise in the discharge. This must be considered as the main unresolved question in the study of jökulhlaups. Both these types of jökulhlaups appear from our observations to be associated with widespread subglacial ponding of flood water at preferred locations along the flood path where a substantial part of the flood volume is temporarily accumulated and released. This implies considerable challenges for quantitative understanding and modelling of jökulhlaups as detailed modelling of widespread lifting, and related dynamic effects on the movement of the overlying ice, requires an integrated modelling of subglacial water flow and ice movement that has not been attempted so far in jökulhlaup studies.

CONCLUSIONS

We have presented several different types of movement of the overlying ice in response to jökulhlaups. They include rapid vertical uplift, often followed by (partial) subsidence, during the initiation of the flood, uplift and subsequent subsidence that takes place during the entire flood, horizontal acceleration, temporary change in horizontal ice-flow direction by up to tens of degrees and temporary reverse ice motion during uplift and subsidence creating backward loop motion. These observations can be interpreted in terms of the development of the jökulhlaup at the bed of the glacier. The observations of jökulhlaups at Skeiðarárjökull in 2004 and 2010 indicate subglacial water pressure higher than overburden, lifting and lateral water storage in typical slowly rising jökulhlaups. The rapid uplift at the beginning of jökulhlaups in Skaftá indicates the passage of a subglacial pressure wave that forms the initial flood path of rapidly rising jökulhlaups. We, furthermore, interpret observations from a rapidly rising jökulhlaup in Skaftá in August 2008 such that parts of the initial flood path were formed by conduit-melt–discharge feedback as is considered typical for slowly rising jökulhlaups. The classification of jökulhlaups into two distinct categories, rapidly and slowly rising, therefore, seems to be too simple. Jökulhlaup behaviour appears to be more of a spectrum. Extremely rapidly rising jökulhlaups, such as the large November 1996 jökulhlaup from Grímsvötn and large jökulhlaups from the eastern Skaftá cauldron, where most of the initial flood path is formed by propagation of a subglacial pressure wave (Jóhannesson, 2002), will then form one end of the spectrum. Slowly rising jökulhlaups, confined in a conduit enlarged by melt–discharge feedback for the whole flood path, are then at the other end of the spectrum. The slowly rising jökulhlaups from Grímsvötn after the November 1996 event, which rise much faster than jökulhlaups before November 1996 as discussed above, and rapidly rising jökulhlaups from the western Skaftá cauldron that rise much more slowly than large jökulhlaups from the eastern cauldron, are examples of events that seem to be of an intermediate type between the two extremes in the spectrum of rapidly and slowly rising jökulhlaups.