Localized uplift of Vatnajökull, Iceland: subglacial water accumulation deduced from InSAR and GPS observations

Abstract

We report on satellite and ground-based observations that link glacier motion with subglacial hydrology beneath Skeiðarárjökull, an outlet glacier of Vatnajökull, Iceland. We have developed a technique that uses interferometric synthetic aperture radar (InSAR) data, from the European Remote-sensing Satellite (ERS-1/-2) tandem mission (1995–2000), to detect localized anomalies in vertical ice motion. Applying this technique we identify an area of the glacier where these anomalies are frequent: above the subglacial course of the river Skeiðará, where we observed uplift of 0.15–0.20 m d⁻¹ during a rainstorm and a jökulhlaup, and subsidence at a slower rate subsequent to rainstorms. A similar pattern of motion is apparent from continuous GPS measurements obtained at this location in 2006/07. We argue that transient uplift of the ice surface is caused by water accumulating at the glacier base upstream of an adverse bed slope where the overburden pressure decreases significantly over a short distance. Most of the frictional energy of the flowing water is therefore needed to maintain water temperature at the pressure-melting point. Hence, little energy is available to enlarge water channels sufficiently by melting to accommodate sudden influxes of water to the base. This causes water pressure to exceed the overburden pressure, enabling uplift to occur.

Introduction

In recent decades, various studies have described the physical properties of subglacial water flow (e.g. Röthlisberger, 1972; Walder, 1986) and how basal sliding of glaciers is influenced by the configuration of subglacial drainage (e.g. Kamb, 1987) and its water pressure (e.g. Weertman, 1957). In order to understand recent observations of accelerated mass loss of arctic glaciers through changes in ice flow (e.g. Rignot and Kanagaratnam, 2006), the role of hydraulic processes in glacier motion is currently under revision (Van de Wal and others, 2008; Shepherd and others, 2009; Bartholomew and others, 2010; Schoof, 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. Iken and Bindschadler, 1986; Fudge 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. Iken and others, 1983; Sugiyama 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. Fatland and Lingle, 2002; Gray 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 (Schoof, 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 (Fricker and others, 2007; Bell, 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 (Iken and others, 1983; Anderson 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 (Bacher 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

Research area

Vatnajökull (Fig. 1) is Iceland’s largest ice cap, covering an area of 8100 km², with an average ice thickness of 380 m (Björnsson and Pálsson, 2008). Excluding a thin surface layer where the temperature varies seasonally, the ice cap is at the pressure-melting temperature, and thus is classed as a temperate or warm-based glacier (Björnsson and Pálsson, 2008). During the 20th century, Vatnajökull decreased by ∼10% (∼300 km³), corresponding to ∼1 mm global sea-level rise (Björnsson and Pálsson, 2008). Systematic mass balance and meteorological measurements have been conducted on Vatnajökull since 1991. Between 1991 and 1995 the balance was positive, but since 1995 it has been negative, with average balance corresponding to ∼−1 m a⁻¹ (Björnsson and Pálsson, 2008).

The main focus of this study is the ablation area of the outlet glacier Skeiðarárjökull, which is a surge-type glacier with seven surges reported since the late 18th century (Björnsson and others, 2003). The glacier terminus retreated 0.5–3 km in the 20th century (Generalstabens topografiske Afdeling, 1905; Bacher and others, 2001) and has been retreating continuously since its last surge in 1991 (Björnsson, 1998). Radio-echo soundings were carried out on Skeiðarárjökull in 1993–98 (Fig. 1). The survey revealed complex subglacial topography, with ice up to 940 m thick in the accumulation area (unpublished data from the Institute of Earth Science, University of Iceland). Deep troughs reaching far below sea level, with a ridge between them, exist underneath the eastern and the western part of the glacier tongue (the eastern one reaches 240 m below sea level).

Significant temporal variations have been observed in the surface velocity field of Skeiðarárjökull. Short-term variations related to intense rainfall, glacier melting and jökulhlaups (glacier outburst floods) lasting for hours or days have been reported (Magnússon and others, 2006, 2007), as well as deceleration linked to long-term changes in subglacial water flow affecting the surface velocity for several years (Magnússon and others, 2010).

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 (Hanssen, 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 (Hanssen, 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 (Bacher and others, 2001) and 2004 (Berthier, 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 (Guðmundsson and others, 2003, 2009). 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 (Guð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. Fatland 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 (Magnússon and others, 2006, 2007, 2010) 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 Reeh and others (2003) (for details see also Magnússon and others, 2007; Magnú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. Paterson, 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. Paterson, 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⁻¹ (Bjö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 (Bjö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 (Magnú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. Hooke 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. Truffer 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 (Truffer and others, 1999), which is comparable to grain-size distribution of till recovered from the forefield of Skeiðarárjökull (Waller 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; Truffer 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 (Iken and others, 1983; Anderson 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 Bartholomaus 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 Flowers (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 (Röthlisberger, 1972; Nye, 1976; Spring 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 (Clarke, 2005). The subglacial water potential, Φ, is given by (e.g. Paterson, 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. Sugiyama 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

(Clarke, 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 (Alley and others, 2003) and drainage of jökulhlaups (Roberts 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 (Walder, 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 Rö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; Fricker and others, 2007; Bell, 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.