Introduction

In recent decades, various studies have described the physical properties of subglacial water flow (e.g. Reference RöthlisbergerRöthlisberger, 1972; Reference WalderWalder, 1986) and how basal sliding of glaciers is influenced by the configuration of subglacial drainage (e.g. Reference KambKamb, 1987) and its water pressure (e.g. Reference WeertmanWeertman, 1957). In order to understand recent observations of accelerated mass loss of arctic glaciers through changes in ice flow (e.g. Reference Rignot and KanagaratnamRignot and Kanagaratnam, 2006), the role of hydraulic processes in glacier motion is currently under revision (Reference Van de WalVan de Wal and others, 2008; Reference Shepherd, Hubbard, Nienow, McMillan and JoughinShepherd and others, 2009; Reference Bartholomew, Nienow, Mair, Hubbard, King and SoleBartholomew and others, 2010; Reference SchoofSchoof, 2010). Studies linking basal water-pressure measurements and glacier velocities are commonly based on observations from discrete boreholes and stake velocities on the ice surface (e.g. Reference Iken and BindschadlerIken and Bindschadler, 1986; Reference Fudge, Harper, Humphrey and PfefferFudge and others, 2009). Point measurements of vertical surface displacement have shown that subglacial water accumulation is commonly related to periods of intense rainfall or heightened melt rates (e.g. Reference Iken, Röthlisberger, Flotron and HaeberliIken and others, 1983; Reference Sugiyama and GudmundssonSugiyama and Gudmundsson, 2004), often resulting in enhanced horizontal motion.

Beyond discrete, point measurements of surface motion, interferometric synthetic aperture radar (InSAR) provides a large-scale spatial overview of glacier surface displacement. Large spatial phase anomalies in InSAR data, often referred to as ‘bull’s-eye’ patterns, have been interpreted as signs of localized water accumulation or depletion (e.g. Reference Fatland and LingleFatland and Lingle, 2002; Reference Gray, Joughin, Tulaczyk, Spikes, Bindschadler and JezekGray and others, 2005). Understanding how, why and when water storage takes place will provide further knowledge of the mechanisms producing high basal water pressure. Recent work indicating that increased variability in subglacial water flow, rather than increased volume of subglacial water flow, is speeding up arctic glaciers (Reference SchoofSchoof, 2010) highlights the need for further investigations of water-pressure transients and consequent speed-up events. Studying ‘bull’s-eye’ patterns caused by short-lived events (days) in the subglacial hydrology of ice caps and mountain glaciers may also help to explain temporal formation of subglacial lakes beneath ice sheets (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference BellBell, 2008), despite different time and length scales.

In this paper, InSAR displacements over 24 hour time-spans (European Remote-sensing Satellite (ERS-1/-2) tandem data from 1995–2000) are utilized to study bed separation due to subglacial water accumulation. We introduce a method to identify anomalies in the vertical ice motion using InSAR data of glaciers with variable ice flux, and use it to document episodes of ice-surface uplift and subsidence at Skeiðarárjökull, an outlet glacier of Vatnajökull, Iceland (Fig. 1). We identify a region where the most pronounced and frequent uplift and subsidence occurs. The surface motion pattern was also observed by differential GPS (DGPS) measurements in 2006/7. We argue that this localized uplift and subsidence cannot be explained by the conventional model of subglacial water cavities enlarged by increased sliding rate (Reference Iken, Röthlisberger, Flotron and HaeberliIken and others, 1983; Reference AndersonAnderson and others, 2004). Instead we propose that the uplift takes place due to hydraulic jacking on an adverse bed slope, where the gradient in water pressure is high and most of the frictional energy of water flow is used to maintain the water at the pressure-melting point. This slows expansion of water-filled ice channels, causing restriction of the water flow.

Fig. 1. The study area, Skeiðarárjökull outlet at Vatnajökull, Iceland. The triangles marked F and S show the locations of the weather stations Faghólsmýri and Skaftafell, respectively. The color map shows the bedrock elevation (previously unpublished data from the Institute of Earth Sciences) based on radio-echo sounding in 1994, 1997 and 1998 (red points and profiles). The contours show surface elevation (Reference Bacher, Bludovsky, Dorrer, Münzer, Oltan and GründigBacher and others, 2001) in m a.s.l. Solid triangles signify the locations of the GPS stations (see Fig. 4a). The light blue areas at the eastern side of the glacier indicate marginal lakes

Skeiđarárjökull Insar Data

Data processing

An InSAR image reveals the difference in distance of a radar beam in the radar line of sight (LOS) between two radar images, acquired from two similar positions (Reference HanssenHanssen, 2001). The phase difference of the radar beam in the two images reveals the relative change in LOS distance across the image. To use repeat-pass satellite InSAR images for motion analyses, the geometric effect due to a slight difference in orbits, usually referred to as the topographic phase, must first be removed using orbit parameters and a digital elevation model (DEM) of the survey area (Reference HanssenHanssen, 2001). In our study with the ERS-1/-2 tandem data (24 hour observation repeat time) from 1995 to 2000, this was carried out using DEMs from 1997 (Reference Bacher, Bludovsky, Dorrer, Münzer, Oltan and GründigBacher and others, 2001) and 2004 (Reference BerthierBerthier, 2005). Topographic correction of InSAR scenes between 1997 and 2000 was achieved by interpolating the DEMs, whereas topographic correction of scenes before 1997 was extrapolated from the models. The surface velocity in the LOS direction, v _LOS, over the time period, Δt, between two radar observations, is derived from the residual phase difference, φ, by using

(1)

where λ is the radar wavelength (5.6 cm for the ERS). InSAR scenes obtained during the summer were adjusted using a degree-day ablation model to compensate for surface lowering due to melting during the observation repeat time (24 hours). Studies on Icelandic ice caps during the last two decades show that degree-days inferred from meteorological stations at some distance from the glacier yield better estimations of the ablation than degree-days observed close to or within the boundary layer of the glacier (Reference Guðmundsson, Björnsson, Pálsson and HaraldssonGuðmundsson and others, 2003, Reference Guðmundsson, Björnsson, Pálsson and Haraldsson2009). Temperature observations (24 hour average) from Fagurhólsmýri meteorological station, T _F (data from the Icelandic Meteorological Office), were therefore used rather than those from the Skaftafell meteorological station closer to the glacier (Fig. 1) in the degree-day model:

(2)

where the value for the ablation model coefficient, ddf, of 6 mm ^°C⁻¹ and the temperature lapse rate, γ, of 0.006°C m⁻¹ were based on observations from Breiðamerkurjökull outlet glacier, ∼40 km east of Skeiðarárjökull (Reference Guðmundsson, Björnsson, Pálsson and HaraldssonGuðmundsson and others, 2003). Z ₀ is the surface elevation at the point of estimated melting, and Z _F is the elevation of Fagurhólsmýri weather station.

The residual phase term representing the motion of the ice surface in the LOS direction may still include errors (in general, <1.4cm d⁻¹) due to inaccurate melting corrections and temporal changes of atmospheric propagation properties, mainly related to water vapor. These errors are, however, expected to vary smoothly in space at regional scales and are therefore of little significance when extracting localized spatial anomalies from InSAR data. We therefore cautiously estimate an uncertainty of 1.4 cm d⁻¹ for the localized spatial anomalies in LOS velocity (corresponding to ∼1.5 cm d⁻¹ if interpreted as vertical motion) due to errors in the InSAR data.

Identifying vertical motion anomalies

Localized anomalies in InSAR data, appearing as bull’s-eye patterns and commonly interpreted as localized vertical motion, have previously been detected by calculating the difference between two interferograms from the same orbit (e.g. Reference Fatland and LingleFatland and Lingle, 2002). This approach is possible only when the two interferograms show approximately the same ice motion outside the bull’s-eye pattern. We adopt a new approach here to identify spatial anomalies in the glacier flow field that can be applied to glaciers with variable ice flux like Skeiðarárjökull (Reference Magnússon, Rott, Björnsson, Roberts, Berthier, Pálsson and LacosteMagnússon and others, 2006, Reference Magnússon, Rott, Björnsson and Pálsson2007, Reference Magnússon, Björnsson, Rott and Pálsson2010) and that enables the use of data from different orbits. Instead of studying the direct difference between two InSAR observations, we take into account velocity changes and derive the LOS velocity anomaly terms at time t and location (x,y) as

(3)

where the ice flow rate factor, α, which we define later, describes the spatially smooth acceleration/deceleration of the flow field relative to our reference LOS velocity, .

Instead of using a measured interferogram from the same orbit as , we first calculate a three-dimensional (3-D) reference velocity field for the glacier surface, V_ref (Fig. 2a), and derive by taking the scalar projection of V_ref onto a given LOS. To obtain V_ref we follow the approach of Reference Reeh, Mohr, Madsen, Oerter and GundestrupReeh and others (2003) (for details see also Reference Magnússon, Rott, Björnsson and PálssonMagnússon and others, 2007; Reference MagnússonMagnússon, 2008), combining LOS velocities from N orbits with the mass continuity equation to obtain the 3-D V_ref (Fig. 2a) for the glacier surface. If the InSAR observation represents approximately the same velocity field and includes observations from both ascending and descending orbits, the three velocity components of V_ref = (v_x, v_y, v_z ) can be derived from the equations (as before, all variables are functions of t, x and y, but the attachment (t,x,y) is skipped for abbreviation)

(4)

(5)

where φ is the radar azimuth angle, measured relative to the x-axis in counterclockwise direction, θ is the radar incidence angle relative to a level surface, Z ₀ is the surface elevation, H is the glacier thickness and F is the ratio of the average velocity in a vertical ice column to the surface velocity, which in general should vary between 0.8 (pure ice deformation) and 1 (pure sliding) (e.g. Reference PatersonPaterson, 1994). Assuming equal contribution of basal sliding and ice deformation to the surface ice motion, F is approximated as 0.9. We apply this to Skeiðarárjökull during dry and cold weather conditions, when we do not expect significant temporal variation in the glacier flow field and the vertical motion caused by water accumulation or depletion should be absent. For this we use two ascending interferograms (obtained on 29–30 December 1995 and 2–3 February 1996) and two descending interferograms (obtained on 27–28 December 1995 and 31 January–1 February 1996). Despite the vertical component of V_ref being relatively inaccurate, the absolute value is generally <1 cm d⁻¹, corresponding to a small fraction of v _LOS for most of the study area. The left side of Equation (5) (with an opposite sign) is referred to as the emergence velocity, v _em (e.g. Reference PatersonPaterson, 1994). It reveals the thickening (or thinning referred to as submergence velocity) of the ice due to ice dynamics at a point fixed in space and should therefore equal the balance at each point with an opposite sign for a glacier in a steady state. v _em derived from V_ref (Fig. 2b) is positive for almost the entire study area, revealing values commonly between 0 and 2 cm d⁻¹ corresponding to 0–7 m a⁻¹. Only limited balance measurements have been carried out on Skeiðarárjökull, but the net balance in the ablation area of Vatnajökull generally varies between 0 and −12 m a⁻¹ (Reference Björnsson and PálssonBjörnsson and Pálsson, 2008). v _em on Skeiðarárjökull’s ablation zone should be somewhat lower in magnitude than the balance since it has been continually retreating since the last surge in 1991 (Reference BjörnssonBjörnsson, 1998). The emergence velocity, v _em, and hence also the vertical velocity component of V_ref is therefore more realistic than with the surface-parallel approach, which assumes v _em = 0 and is commonly used for estimating 3-D surface velocity fields of glaciers from InSAR data.

Fig. 2. The horizontal components of (a) the reference velocity field, V_ref, and (b) the corresponding emergence velocity, v _em, for Skeiðarárjökull outlet.

We now define the ice flow rate factor, α:

(6)

where A is a disk 6 km in diameter (∼tenfold the glacier thickness above the troughs shown in Fig. 1). Locations where or where the unfiltered ratio exceeded 10 (the latter criterion usually related to areas where the ice flow is nearly perpendicular to the radar LOS) are ignored within A. An ice flow rate factor α < 1 indicates that the glacier is moving slower than in V_ref, while α > 1 indicates the opposite. Due to the continuity of the ice flow, we expect the relative changes in the horizontal ice flow to vary evenly over the glacier surface. This means that the term (from Equation (3)) should compensate for most of the variation in the v _LOS caused by changes in the horizontal flow. The residual caused by horizonatal motion that does not account for should be particularly small on Skeiðarárjökull due to its southward ice motion. At the latitude of the study area (∼64° N), ERS InSAR data are about nine times more sensitive to vertical motion than to surface motion southwards. We therefore estimate vertical motion anomalies by projecting onto the vertical axis (θ is the incidence angle of the radar):

(7)

This approach assumes that the horizontal direction of the ice flow is constant. The errors produced by this approximation will be discussed in relation to flow direction variations observed in our GPS observations. In case of events where significant vertical motion occurs over areas with dimensions comparable to the disk, A, the flow rate factor, α, will be skewed in that area, which smooths the estimated anomalies causing the peak magnitude of the estimated vertical motion to be underestimated.

Observed Uplift and Subsidence

The method described above was applied to 35 interferograms from both ascending and descending orbits of Skeiðarárjökull outlet from the end of May 1995 to February 2000 (Fig. 3). At two locations, anomalies in vertical motion were observed repeatedly. One site is near a series of marginal lakes on the eastern side of Skeiðarárjökull (Figs 1 and 3), which appear to drain annually in jökulhlaups. In a map from the US Defense Mapping Agency and the Iceland Geodetic Survey (DMA and IGS, 1990), the lakes have identical surface elevations, suggesting they are subglacially connected, causing floating of the adjoining ice mass at high lake level. The other site of localized vertical motion is above the estimated subglacial path of the river Skeiðará (derived from water potential assuming full ice overburden pressure; explained in more detail in the Discussion), and covers an area across the entire eastern trough underneath the glacier (Fig. 1). This is the location where significant uplift (∼15 cm d⁻¹) was observed during a jökulhlaup in March 1996 (Reference Magnússon, Rott, Björnsson and PálssonMagnússon and others, 2007). We also observe at this location >20 cm d⁻¹ uplift in a scene from 17–18 January 1996. During these two days, rainfall of 8 and 3 mm d⁻¹ was measured, respectively, in Skaftafell ∼5 km east of Skeiðarárjökull (data from the Icelandic Meteorological Office). In both cases, the glacier flowed much faster than V_ref, reflected in the calculated flow rate factor, α, between 2 and 3. This corresponds to 24 hours of horizontal motion of 0.7–1.0 m d⁻¹ in the vicinity of the uplift area, whereas V_ref is ∼0.35 m d⁻¹. If all the uplift in the affected area (∼8 km²) in the latter case was due to water accumulation, it would correspond to ∼10 m³ s⁻¹ discharge averaged over 24 hours, between one-third and one-half the typical winter drainage in the river Skeiðará (personal communication from S. Zóphóníasson, 2006). We also observe slower but significant subsidence (<0.15 m over 24 hours) at the same location in six scenes during or following periods of significant rainfall (5–30 mm d⁻¹). For these scenes, α is generally close to 1. In some of these scenes the observed subsidence can be traced 6–8 km along the subglacial course of the river Skeiðará (Fig. 3d).

Fig. 3. Calculated in four InSAR scenes. (a) The local uplift, , during winter rainfall in the course of the river Skeiðará. (b) Uplift at the same spot during the beginning of a glacier outburst flood (jökulhlaup). (d, e) Significant subsidence 1 day after 10–30 mm d⁻¹ rain. (c) The average vertical motion at each location for all the 35 InSAR scenes where . (f) Same as (c), but for all scenes where .

In Situ Observations of Glacier Surface Displacement

Three GPS stations were mounted on Skeiðarárjökull in the spring of 2006 (Fig. 4). One of them (SKE1) operated until spring 2007, with a 4 month break during the winter. The location of the station was chosen based on the results of the InSAR data. SKE1 was placed within the area of frequent uplift and subsidence (Fig. 3). The other two stations (SKE2 and SKE3) were located outside this area and were operated until autumn 2006. The instruments were fixed on a low-lying platform on the glacier surface. The measured vertical motion therefore includes subsidence caused by melting. The stations recorded their position every 15 s and the data were corrected kinematically using three base stations at Skrokkalda (SKRO), Grímsfjall (GFUM), located on a nunatak in central Vatnajökull, and Höfn (HOFN) 40–95 km from the survey area (Fig. 4a; further information on these base stations is available at the Icelandic Meteorological Office web page, http://hraun.vedur.is/ja/englishweb/gps.html). The software used for the data processing was Trimble Geomatics Office v. 1.63. The GPS records were processed three times, each time using one of the three base stations independently. Due to the large distances between measurement stations and base stations, the solution did not remain stable at all times. All data were therefore manually reviewed, and unrealistic jumps in the solutions (which usually only appeared in one of the base station solutions at a time) were masked out. The solutions obtained with base stations Höfn and Skrokkalda were then shifted towards the solution obtained for Grímsfjall to compensate for any long-term difference between the solutions. At this stage, a single record was obtained from the three using a weighted average, where solutions from Grímsfjall had the weight 4, Höfn 2 and Skrokkalda 1 (the weight based on apparent solution quality).

Fig. 4. (a) The location of the GPS stations deployed on Skeiðarárjökull and the base stations used for the differential correction. The blue triangle, S, shows the location of the weather station Skaftafell. The red curve shows the estimated subglacial course of the river Skeiðará. (b) The horizontal ice motion at SKE1, SKE2 and SKE3 from April 2006 to May 2007. (c) The elevation change of the stations relative to initial location of each segment minus the vertical motion due to the local surface slope and melting. (d) The rainfall in Skaftafell, ∼5 km east of Skeiðarárjökull, during the GPS survey period as well as the estimated melting at SKE1 from a degree-day model.

The vertical motion due to the local surface slope, measured specifically along the path of each station, and melting estimated from a degree-day model was subtracted from the elevation curves in Figure 4c. Below we describe in detail the applied melting correction. The temperature data applied in the model (Equation (2)) were from the Fagurhólsmýri meteorological station, using a temperature lapse rate, γ, of 0.006°C m⁻¹.The ablation model coefficients, ddf, were derived from the melting measured over a 2 month period at SKE1 and SKE2 in summer 2007. The derived ablation model coefficients corresponded to 6.4 mm ^°C⁻¹ lowering over 24 hours at SKE1 and 5.7 mm ^°C⁻¹ at SKE2. For SKE3, located at a similar elevation to SKE2, the latter value was adopted. These values agree well with observations from Breiðamerkurjökull, mentioned previously, despite the short observation period which may affect the result to some degree. We estimate snow accumulation based on temperatures and precipitation in Skaftafell, but this never exceeded 10 cm w.e. Snow accumulation in the lower area of Skeiðarárjökull is low in general and the glacier is often snow-free during most of the winter. If our estimate resulted in snow accumulating on the ice surface, on which the GPS platforms were lying, then snowpack had to be melted (applying the same degree-day model) before the melting correction for the GPS station became active.

The derived melting correction was typically ∼5 cm d⁻¹ for all stations during summer 2006 and it never exceeded 7 cm d⁻¹. For other seasons the melting correction was generally 0–5 cm d⁻¹. Only a single melting correction value was calculated for each day using the average temperature over 24 hours at the meteorological station. The actual melt rate, however, varies diurnally, which in addition to the long baselines of the DGPS measurements may add significant noise to the derived elevation over short time periods. The elevation curves were therefore filtered using 6 hour running averages. The filtered elevation curve of SKE3 (Fig. 4c), which shows insignificant vertical motion, reveals noise with amplitude ∼2 cm, indicating the relative accuracy of the elevation curves. The horizontal positions were filtered with 3 hour running averages, since the horizontal component has less noise than the vertical component. We generally expect the horizontal displacement calculated over 24 hours (Fig. 4b) to be accurate to within 2 cm.

The observed displacement at SKE1 (Fig. 4b and c) shows correlation of both the vertical and the horizontal motion with increase in estimated water flow to the glacier base due to rainfall and/or changes in glacier melt rate (Fig. 4d). During these events we observe horizontal acceleration of the glacier surface from 0.3–0.4 m d⁻¹ to values often exceeding 1.0 m d⁻¹, coinciding with uplift of the glacier surface. This lasts for a few days, with total uplift often between 0.4 and 1.2 m, and exceeds 0.5 m over 24 hours in the extreme case (October 2006). The surface is commonly still rising after the glacier starts decelerating. This is followed by subsidence at a rate generally slower than the rise, but in some cases lasting for >1 week. During the subsidence period, approximately half of the total uplift during the event is reversed. For most of the subsidence period, the horizontal motion of the station is at normal rate (0.3–0.4 m d⁻¹). Most of these uplift events occur outside the main melting season. We only observe one summertime event, during intense rainfall in early July 2006. This is partly due to more extreme changes in water input from the surface for single events between autumn and spring, as demonstrated by Figure 4d which shows the estimated water input from the surface during the GPS observation period. A well-developed system of water channels during midsummer is also more likely to be able to accommodate changes in surface water injection.

At station SKE3, outside the area where anomalies in the vertical motion were frequently observed in the InSAR data, we note the effect of increased water input on horizontal motion, changing the 24 hour displacement from 0.15 m to 1 m in an extreme case (Fig. 4b). The subsequent 24 hour uplift (Fig. 4c) is only ∼5 cm and less than that for other speed-up events. The horizontal motion of SKE2 seems to show a pattern similar to that of SKE1 and SKE3. However, the record is incomplete due to power shortage at SKE2. Unfortunately the data from SKE2 cover only smaller uplift events in July and September; during these events no clear uplift is observed at this location.

The GPS data show some dependence of the flow direction on horizontal velocity, causing the ice flow to diverge less to the sides during speed-up events (Fig. 5a). This adds some uncertainty to the estimated uplift in Figure 3. At the location of SKE1 a deviation of 5° in the flow direction from the presumed flow angle (derived from Fig. 5a) would produce a phase signal in the interferograms which we would misinterpret as ∼3 cm d⁻¹ uplift in Figure 3a and b (compared to 15–20 cm d⁻¹ maximum uplift), but much less in Figure 3d and e. Only a minor fraction of the phase signal anomalies observed in the interferograms can therefore be caused by changes in ice flow direction.

Fig. 5. (a) The flow direction of SKE1 as a function of the horizontal velocity (velocity dependence also observed in SKE2 and SKE3). The red curve shows, for comparison, how the bearing would vary given a constant deformation velocity of 0.3 m d⁻¹ with 153° bearing while everything exceeding 0.3 m d⁻¹ is due to sliding with 161° bearing. (b) The strain rates (green diamonds) between SKE1 and SKE3 as a function of the horizontal velocity at SKE1 for the common observation period. The yellow diamonds show the strain rates between the same stations during a period of subsidence at SKE1 in September and October 2006.

Discussion

The InSAR observations (1995–2000) and the GPS data (2006–07) show frequent anomalies in the vertical motion in an area above the subglacial course of the river Skeiðará, ∼9 km from the glacier margin. In the following discussion, we argue that significant bed separation took place in this area. However, we are not able to directly separate these anomalies quantitatively into motion due to vertical strain and due to water accumulation and depletion. The GPS stations were too few and too far apart for such an analysis and we lack constraints for the InSAR data due to the high horizontal flow rate. The relation between the flow direction and the horizontal velocity observed in the GPS data (Fig. 5a) further complicates estimates of vertical motion due to strain.

To calculate uplift produced by strain, we take the horizontal component of V_ref linearly scaled to fit velocities observed at the GPS stations SKE1 and SKE3 during an event when the highest strain was observed between the two stations (Fig. 5b). From Equation (5) (with F = 1 for simplicity) we derive the emergence velocity revealing the ice thickness change caused by strain as

(8)

where V = (vx, vy) is the scaled velocity field (Fig. 6). Due to the crude estimate for V, we can only consider this as an order-of-magnitude estimate. Since the flow divergence (the right side of the equation) is independent of the chosen coordinate system, we can obtain the same result by choosing the x-axis in the direction of a flowline. The two terms of flow divergence would then reveal uplift caused by longitudinal (the x term) and transverse strain (the y term) as commonly estimated in glacier studies using GPS or other discrete velocity measurements (e.g. Reference Hooke, Calla, Holmlund, Nilsson and StroevenHooke and others, 1989).

Fig. 6. Order-of-magnitude estimate of the emergence velocity during speed-up event when strain rate between SKE1 and SKE3 is at maximum observed value.

This simple calculation indicates that uplift due to strain can generally not be excluded as a possible source for uplift of the same order of magnitude (>0.1 m d⁻¹) as observed at SKE1. However, if the uplift was caused entirely by vertical strain, it would not be limited to the vicinity of SKE1 alone as the InSAR data suggest (Fig. 3a–c) and it would be unlikely to produce a bell-shaped pattern across the trench underneath SKE1 (Fig. 1). We also observe from the GPS data at SKE1 that typically around half the vertical motion occurring during uplift events is reversed during the subsequent subsidence period (Fig. 4c). During the subsidence period, the horizontal motion rate at SKE1 is back to normal. Calculated strain rates between SKE1 and SKE3 during periods of subsidence at SKE1 (the data from SKE2 are too sparse for such calculation) are in general roughly the same as observed at other times (Fig. 5b). Extensile strain rates are therefore unlikely to cause the observed subsidence.

Expansion of subglacial till due to decrease in the effective basal water pressure is also unlikely as a source of the rapid uplift observed. Reference Truffer, Echelmeyer and HarrisonTruffer and others (2001) modeled 10 cm uplift due to till expansion at a borehole on Black Rapids Glacier, Alaska, USA. This high uplift was, however, derived assuming much higher till compressibility than expected considering the grain size distribution of the till recovered from the borehole (Reference Truffer, Motyka, Harrison, Echelmeyer, Fisk and TulaczykTruffer and others, 1999), which is comparable to grain-size distribution of till recovered from the forefield of Skeiðarárjökull (Reference Waller, van Dijk and KnudsenWaller and others, 2008). We expect the variations in the effective water pressure (the main variable in the till expansion) to be of similar order since the ice thicknesses at SKE1 (∼700 m; Fig. 1) and at the borehole on Black Rapids Glacier (∼500 m; Reference Truffer, Motyka, Harrison, Echelmeyer, Fisk and TulaczykTruffer and others, 1999) are comparable.

Based on the above arguments, we expect that most of the reversed vertical motion is due to depletion subsequent to water accumulation. The InSAR scenes where subsidence can be traced along the subglacial course of the river Skeiðará suggest water drainage towards the main channel of the river from the surrounding areas, further supporting this interpretation.

The magnitude of the uplift events and the lack of spatial correlation between the uplift and the sliding velocity observed in the InSAR data and further supported by the GPS data suggest that the localized uplift on Skeiðarárjökull cannot be explained with the conventional model of subglacial water cavities enlarged by increased sliding rate (Reference Iken, Röthlisberger, Flotron and HaeberliIken and others, 1983; Reference AndersonAnderson and others, 2004). Based on our observations, we suggest that outflow of water from channels to distributed flow is responsible for the uplift events. This is in better agreement with the observations of Reference Bartholomaus, Anderson and AndersonBartholomaus and others (2008).

Water passing the location of most frequent uplift (Figs 3 and 4) is drained from a 1070 km² catchment area upstream. Due to high discharge from such a large catchment area, we expect the hydrological system to be, in general, governed by channel flow at this location, at least during the summer. During uplift events observed at SKE1, the estimated total water input from the surface (precipitation plus melting from Fig. 4d) never exceeded 10 cm d⁻¹. We cannot see how water accumulation producing an uplift rate of decimeters per day during rainfall or an increase in melting can be explained without effective water transfer towards the area of uplift, even during winter. The large meltwater catchment upstream as well as the maritime conditions at southern Vatnajökull, causing a large part of the study area to be snow-free throughout most of the winter, may maintain channelized basal flow throughout the year.

Modulation of subglacial water flow, as described by Reference FlowersFlowers (2008), suggests that channel flow adapts ineffectively to increases in water flow on locally adverse bed slope, as found underneath SKE1. To elaborate how adverse bed slope may influence the channel flow during events of increased water flux, we consider the growth of a channel (semicircular) cross-sectional area, S, given by the differential equation (Reference RöthlisbergerRöthlisberger, 1972; Reference NyeNye, 1976; Reference Spring and HutterSpring and Hutter, 1981)

(9)

where Q _w is the water discharge in the channel, ρ_{i ce} the density of ice, L the latent heat of fusion, p _ice the ice overburden pressure, p _w the water pressure in the channel, B the ice-flow parameter and n the Glen’s flow law exponent. The variable, μ, we refer to as the melt-rate ability of the channel walls. It determines how much energy is available to melt the ice from the channel wall. If p _w approaches the ice overburden pressure, p _ice, the creep closure of the channel (the second term on the right side of Equation (9)) diminishes, and hence the channel water discharge, Q _w, and the melt-rate ability, μ, control how fast the size of tunnel can change. Assuming that heat is transferred instantaneously from the water to the ice walls, maintaining the water at its pressure-melting point, μ can be defined as

(10)

where s is the distance along the channel, C_t the pressure-melting coefficient and c _w the heat capacity of water (Reference ClarkeClarke, 2005). The subglacial water potential, Φ, is given by (e.g. Reference PatersonPaterson, 1994)

(11)

where Z _b is the elevation of the bedrock above some datum, g the acceleration due to gravity and ρ _w the density of water. In our discussion we assume that the basal water pressure, p _w, relates to the ice overburden pressure, p _ice, as

(12)

where the pressure ratio, K, varies both spatially and with time. Hence

(13)

Borehole studies from temperate glaciers show that during events of high motion triggered by increased water input, subglacial water may be close to the ice overburden pressure (e.g. Reference Sugiyama and GudmundssonSugiyama and Gudmundsson, 2004). Hence, we assume here that K = 1 is reasonable for the speed-up events in Figure 4, and from Equations (11– 13) we derive

(14)

By using c _w = 4220 J kg⁻¹ K⁻¹ and C_t = 7.42 × 10⁻⁸ K Pa⁻¹, μ will be negative for adverse bed slopes where

(Reference ClarkeClarke, 2005). Negative μ means that supercooled water freezes to the ice walls as it flows along the channel. The effect of glaciohydraulic supercooling is important for glacier erosion (Reference Alley, Lawson, Larson, Evenson and BakerAlley and others, 2003) and drainage of jökulhlaups (Reference RobertsRoberts and others, 2002).

We calculate μ from Equation (14) for a 100 m × 100 m grid using surface and bedrock DEMs of Skeiðarárjökull (Fig. 1). We assume that the water flow direction is controlled by the potential, Φ, derived from Equation (13) with K = 1, hence we calculate the derivatives ∂Z ₀/∂s and ∂Z _b/∂s in the direction of −∇Φ. The potential gradient was also used to map watercourses underneath Skeiðarárjökull (Fig. 7), and based on this the main course of the river Skeiðará was located.

Fig. 7. (a) The melt-rate ability, μ, derived from Equation (14) for Skeiðarárjökull (color map). The black curves show watercourses estimated from the potential given by Equation (13), with K = 1. (b) μ along the estimated main course of the river Skeiðará (shown with light blue line in (a)), as well as bedrock and surface elevation. Locations of GPS stations are shown in (a) and (b) for reference.

Even though the assumption K = 1 may not be entirely realistic, the model points at locations at the glacier bed where channel flow may adapt slowly to increased water input from the surface, hence limiting the capacity of the ice tunnels to accommodate for rapidly injected water and causing water to spread out. The estimated μ is close to zero, even slightly negative due to an adverse bed slope, downstream of the area of greatest uplift observed in the InSAR data (Fig. 4). There the water-pressure gradient along the water channels (the second derivative in Equation (10)) is expected to be high. Hence, most of the frictional energy of the water flow is needed to maintain the water at the pressure-melting point, leaving little for melting of the ice walls. This limits the enlargement of the channels in response to suddenly increased water input. Hence, the basal water pressure increases and bed separation takes place upstream of a constriction in the water channels, and pressurized water starts accumulating in a sheet or in cavities (Reference WalderWalder, 1986), in the vicinity of the channels. However, given time after the water flow rate from the surface stops increasing, the channels are expected to gradually expand and water flow through the channels starts behaving more like the steady-state channel flow described by Reference RöthlisbergerRöthlisberger (1972), hence the water pressure decays. Pressurized water in the vicinity of the channels will then start to drain back to the channels, causing the overlaying ice to subside slowly back to its former state as observed from both the InSAR data and the GPS data at SKE1.

In another environment of different scales in time and dimensions, Antarctica, the mechanism described above may help to understand the formation of short-lived subglacial lakes (Grey and others, 2005; Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference BellBell, 2008), at locations where static water pressure will not sustain a lake.

Conclusions

Our study demonstrates the potential of InSAR as a tool for studying links between subglacial hydrology and glacier surface motion. By using InSAR data from the ERS-1/-2 tandem mission, we identify an area of Skeiðarárjökull, an outlet of Vatnajökull, above the estimated subglacial course of its main river, where anomalies in vertical ice motion are frequent; these events take place during rapid injection of water to the glacier base. We observe 0.15–0.20 m d⁻¹ uplift during a rainstorm and one jökulhlaup, followed by a slower subsidence (<0.15 m d⁻¹). In some scenes the subsidence can be traced along the subglacial main course of the river, suggesting that water drains towards the main channel of the river from the surrounding areas. Continuous GPS records on the glacier surface, obtained in the area of the most pronounced and frequent vertical motion anomalies in 2006 and 2007, also show high vertical uplift (0.10–0.50 m d⁻¹) related to episodes of rainfall or increases in melting, followed by slower subsidence. In contrast, vertical motion of this magnitude was not observed in GPS data from two other neighboring stations on the glacier. We argue from the spatial as well as the temporal pattern that the uplift and subsequent subsidence were partly caused by hydraulic uplift of the glacier base followed by subsidence when the water pressure drops. We propose that during events of sudden water influxes to a base of adverse bed slope the subglacial channels conducting water may not expand sufficiently by frictional melting due to near supercooling conditions of the water, causing a constriction (bottleneck) in the channel flow and consequent water accumulation upstream of it. (Fig. 7). Given time, the channels adapt to the basal water flow by expansion, and the pressurized water accumulated upstream drains into the channels, resulting in decreased water pressure, allowing the overlying ice to slowly subside.