Elevation change, mass balance, dynamics and surging of Langjökull, Iceland from 1997 to 2007

Abstract

Glaciers and ice caps around the world are changing quickly, with surge-type behaviour superimposed upon climatic forcing. Here, we study Iceland's second largest ice cap, Langjökull, which has both surge- and non-surge-type outlets. By differencing elevation change with surface mass balance, we estimate the contribution of ice dynamics to elevation change. We use DEMs, in situ stake measurements, regional reanalyses and a mass-balance model to calculate the vertical ice velocity. Thus, we not only compare the geodetic, modelled and glaciological mass balances, but also map spatial variations in glacier dynamics. Maps of emergence and submergence velocity successfully highlight the 1998 surge and subsequent quiescence of one of Langjökull's outlets by visualizing both source and sink areas. In addition to observing the extent of traditional surge behaviour (i.e. mass transfer from the accumulation area to the ablation area followed by recharge of the source area), we see peripheral areas where the surge impinged upon an adjacent ridge and subsequently retreated. While mass balances are largely in good agreement, discrepancies between modelled and geodetic mass balance may be explained by inaccurate estimates of precipitation, saturated adiabatic lapse rate or degree-day factors. Nevertheless, the study was ultimately able to investigate dynamic surge behaviour in the absence of in situ measurements during the surge.

1. INTRODUCTION

Glaciers and ice caps are important components of the world's hydrological cycle. Their mass loss has contributed 226 ± 135 Gt a⁻¹ (0.62 ± 0.37 mm sea level equivalent (s.l.e.) a⁻¹) from 1971 to 2009, and 275 ± 135 Gt a⁻¹ (0.76 ± 0.37 mm s.l.e. a⁻¹) from 1993 to 2009, to the global oceans (IPCC, 2013). Monitoring their changes is key to understanding the impacts of global and regional climate change (e.g. Oerlemans, 1994), the implications for global sea-level rise (e.g. Radić and Hock, 2011; Marzeion and others, 2012; Chen and others, 2013) and the effects on regional and local hydrology, including flooding (Dahlke and others, 2012), river biodiversity (Jacobsen and others, 2012) and water supply to large populations (Hopkinson and Demuth, 2006; Björnsson and Pálsson, 2008; Baraer and others, 2011; Barry, 2011; Bolch and others, 2012).

Glacier and ice-cap mass balance is the key variable that must be monitored because it is changing mass balance that affects global sea levels and regional and local river regimes. Glacier and ice-cap mass balance has traditionally been monitored using the glaciological method based on point stake, pit and probe measurements. Increasingly, glacier and ice-cap mass balance is being assessed using the geodetic approach involving measurement of surface elevation and therefore volume change. Full details on both methods are available in the Glossary of Glacier Mass Balance (Cogley and others, 2011).

The glaciological and the geodetic approaches both have inherent uncertainties (e.g. spatial extrapolation in the former and density assumptions in the latter). The few studies that have been done comparing the two methods for specific glaciers and ice caps for common time periods, have shown differences that are both positive and negative (e.g. Andreassen, 1999; Cogley, 2009; Zemp and others, 2013; Wang and others, 2014). There is a need, therefore, to understand the biases in the two techniques and what controls them, in an attempt to reconcile the two types of measurement (Geist and others, 2005; Hagen and others, 2005; Fischer, 2011).

On glaciers and ice caps, surface elevation changes are not only a function of surface mass-balance changes, but also dynamic changes associated with glacier flow, referred to as flux divergence (Hubbard and others, 2000; Nuth and others, 2012). For example, both negative mass balance and flow acceleration contribute to glacier thinning. Conversely, positive mass balance and flow deceleration both contribute to glacier thickening. The observed surface elevation change signal is a combination of multiple factors (e.g. a thinning surface can in fact result from decelerating ice and a very negative surface mass balance).

Mass balance, velocity and surface elevation are linked through the concept of a glacier's balance flux. The balance flux is defined as ‘The hypothetical horizontal mass flux (dimension [M T⁻¹]) through a vertical cross section that would be equal to the mass balance (usually the climatic mass balance) over the region up-glacier from the cross section’ (Cogley and others, 2011). The balance flux can be divided by the thickness of the glacier to yield a balance velocity, which can then be compared with actual glacier velocity as an indication of the health of a glacier.

Quantifying the various contributions is important for understanding the drivers of elevation change (i.e. surface accumulation/ablation or flow acceleration/deceleration). Dynamic thinning has been identified on outlet glaciers from the Greenland ice sheet and across outlet glaciers and ice streams of Antarctica by comparing measurements of elevation change with surface mass-balance calculations (e.g. Zwally and others, 2005; Pritchard and others, 2009). However, identifying the surface mass balance vs dynamic contributions to elevation changes across glaciers and ice caps has, so far, been attempted by only a handful of studies (e.g. Fischer, 2011; Nuth and others, 2012).

A case where dynamic changes make an important contribution to elevation change is that of surging glaciers. Glacier surges cause elevation changes that are very different from surface mass-balance changes, deflating high-elevation source regions and uplifting low-elevation sink regions as mass is transferred at high velocities over several years and, in the case of ice caps, altering the location of ice flow divides (Björnsson and others, 2003; McMillan and others, 2014). Using elevation change to better define surge extent can aid in a deeper understanding of the other processes (e.g. basal hydrology) important in non-steady glacier flow.

Surge type glaciers are concentrated in particular regions of the globe including Svalbard, Canadian Arctic islands, Alaska, the Karakoram and Iceland (Jiskoot and others, 2000; Sevestre and Benn, 2015). Those in Iceland have been particularly well documented over many decades, largely on the basis of ground-based observations of terminus advance (Björnsson and others, 2003). More recently, remote sensing imagery has been used to identify surge-type glaciers over large regions (e.g. Grant and others, 2009; Paul, 2015; Sevestre and Benn, 2015) and to measure velocity increases associated with particular surges (e.g. Pritchard and others, 2005; Mansell and others, 2012). These last examples use feature tracking or InSAR to resolve horizontal variations in velocity but do not provide a full 3-D description of surge behaviour, including vertical changes. Minchew and others (2015) also measured horizontal flow velocities on Langjökull using InSAR (up to 75 m a⁻¹), but reliable vertical velocity estimates were confounded by moisture-induced phase offsets.

2. OBJECTIVE

The overall objective of this study is to compare measurements of surface elevation change with calculations of surface mass-balance change over a large ice cap, Langjökull, which contains both surge-type and non-surge type outlet glaciers. The study enables us not only to compare the geodetic mass-balance calculations with calculations of the surface mass balance for all the main outlet glaciers as well as the ice cap as a whole, but also to map spatial variations in the dynamic component of surface elevation change, specifically patterns of vertical submergence and emergence. Comparing such patterns across the ice cap enables us to locate and quantify the source and sink areas of the ice cap associated with glacier surge activity.

Using three DEMs generated for the entire ice cap for 1997, 2004 and 2007, we calculate elevation change for two epochs (20 April 1997–15 August 2004 and 15 August 2004–2 August 2007). For the same time periods, we use a degree-day approach to model the climatically-driven surface mass balance across the ice cap. The model is constrained by measurements of summer, winter and net mass balance made at a network of stakes across the ice cap.

3. STUDY SITE

Langjökull (‘Long Glacier’) is Iceland's second largest ice cap (Fig. 1). It is oriented SW–NE in central western Iceland with an area of ~900 km², an elevation range of ~460 to ~1440 m a.s.l., and a volume of ~190 km³, equivalent to 0.5 mm of eustatic sea-level rise (Pálsson and others, 2012). It has an equilibrium line altitude at ~1100 m a.s.l. up to ~1200 m a.s.l. at northern outlets, surface slopes ranging from ~5° on outlet glaciers to ~1° near the summit, and an average slope of 3.4° (Pope and others, 2013). It has a mean ice thickness of ~210 m and a maximum of ~650 m (Björnsson and Pálsson, 2008). Evidence suggests that Langjökull is completely temperate, and the widespread presence of moulins implies that meltwater is freely able to reach the ice-cap bed (Eyre and others, 2005). Precipitation is thought to exert a strong influence over the form and flow of the ice cap, with a steeply decreasing precipitation gradient from south to north (Palmer and others, 2009). This results in fast, steep south-flowing glaciers reaching down to elevations of <600 m and shallower, north-flowing glaciers terminating in broad fronts at higher altitudes.

Fig. 1. Map of Langjökull showing its 2004 topography (WGS84), ice divides, major outlets and the nearby Hveravellir weather station, underlain by a false-color Landsat 8 image from 7 June 2014. Also shown are the mass balance stakes, where winter and summer mass balance have been measured since 1996. The inset is a map of Iceland showing the country's major ice caps including Langjökull.

Termini of Langjökull have been surveyed since 1933. As with other ice masses across Iceland, there has been a significant retreat of Langjökull's main outlet glaciers through much of the 20th century, with an average mass balance of over 50 cm w.e. a⁻¹ over the century (Pálsson and others, 2012). The ice cap is predicted to disappear completely by ~2140 for the A1B scenario of ~2°C (100 a)⁻¹ warming (Björnsson and Pálsson, 2008) due to low-snow accumulation and high-annual temperatures (Björnsson and others, 2002). The three surge-type glaciers (Vestari-Hagafellsjökull, Eystri-Hagafellsjökull and Suðurjökull) have undergone periodic advance and retreat, superimposed on the overall climatically-driven retreat (Sigurðsson, 1998; Björnsson and others, 2003; Palmer and others, 2009). During surges, ice velocities increase over areas of up to 200 km². As Langjökull is temperate, evidence suggests that surges are due to the subglacial hydrological system intermittently switching from a channelized to a distributed system, and thereby allowing faster glacier flow (Björnsson and others, 2003).

In the south, the two major outlets, Vestari-Hagafellsjökull and Eystri-Hagafellsjökull, are separated by the Hagafell Ridge (Fig. 1). Vestari-Hagafellsjökull is ~7 km wide, 25 km long, and bounded on the east by the Hagafell Ridge. Eystri-Hagafellsjökull is ~4 km wide and feeds proglacial lake Hagavatn. It is constrained to the east by the Jarlhettur volcanic ridge, which it overtops in places forming small piedmont lobes in the Jarlhettukvísl Valley (Bennett and others, 2005). Vestari surged in 1971 and 1980 and Eystri surged in 1974 and 1980; whilst observations are based on increased turbidity of outflow and terminus advance, the true duration of the surge behaviour is unknown (Björnsson and others, 2003). In 1998, a surge began in the upper parts of both glaciers, stopping above the terminus of Vestari, but continuing all the way to the snout on Eystri-Hagafellsjökull (Björnsson and others, 2003). No previous surges are documented for either glacier (Sigurðsson, 1998; Bennett and others, 2005; Björnsson and Pálsson, 2008). Past surges of Eystri-Hagafellsjökull formed complex sediment deformation structures in and around lake Hagavatn (Bennett and others, 2000). In addition, Palmer and others (2009) report evidence for elevated velocity and therefore surges of the eastern outlet Suðurjökull (Fig. 1) in 1994 and 1999. This 15 km long, 3 km wide outlet terminates close to proglacial lake Hvítárvatn. Varved lake sediments suggest eight surges occurred between 1828 and 1930 AD with a periodicity of 14 ± 4 a when the glacier terminated in the lake (Larsen and others, 2013). During each surge, the terminus advanced up to ~1.6 km in <2 a.

4. DATA

4.1. Digital elevation models

Three DEMs are used to calculate changing elevation across Langjökull. The first covers the period from April 1997, interpolated from dGPS vehicle tracks across the ice cap spaced ~1 km apart; the average elevation accuracy is estimated by the authors to be <2 m (Pálsson and others, 2012; their Fig. 3). The second DEM is from August 2004, the result of three SPOT5 stereo pairs; elevation accuracy is estimated by the authors to be ~1 m (Pálsson and others, 2012; their Fig. 5). The third is based on airborne lidar data collected in August 2007. Full coverage of the ice cap was not achieved; every other flight line was flown, leaving ~2 km gaps between adjacent flight lines (Pope and others, 2013; their Fig. 3). A new 2007 DEM, gridded to 3 m, was generated by interpolating between lidar flight lines using the corresponding slope and aspect data from the spot-derived 2004 DEM. The elevations of the interpolated areas were adjusted (shifted and tilted) so that they match the 2007 elevation data at the edges of the lidar swaths. This new 2007 DEM has accuracies of 0.25 m in lidar areas and ~3 m in interpolated areas, prompting the authors to estimate an overall accuracy of ~2 m (Jóhannesson and others, 2013). All DEMs were re-gridded to 30 m resolution with common grid positions using bilinear interpolation.

4.2. Mass-balance measurements

In situ winter and summer mass-balance measurements have been conducted on Langjökull since 1996 at 22 stakes distributed across the ice cap (Björnsson and others, 2002; Pálsson and others, 2012) (Fig. 1). Manually produced contours and a 200 m grid are used to calculate the mass balance for the entire ice cap, with error estimated at 5–15% (full details in Pálsson and others, 2012; their Fig. 4). Here we use the point measurements of winter, summer and net mass-balance measurements and the ice-cap average winter, summer and net mass-balance estimates to parameterize a distributed model of surface mass balance for Langjökull.

4.3. Temperature data

Temperature data are used to drive a distributed surface mass-balance model for Langjökull.

Long-term temperature monitoring is undertaken across Iceland by the Icelandic Meteorological Office. The closest weather station to Langjökull is at Hveravellir (64°52.010′N; 19°33.733′W; 641 m a.s.l.), ~10 km from the north-east ice cap margin (Fig. 1). A lapse rate-based extrapolation of temperature from this one point across the ice cap is likely to produce large errors since previous studies have shown the relationship to be non-linear and variable in time and space (Guðmundsson and others, 2003; Flowers and others, 2007).

For this reason, we use a daily 1 km gridded product based on interpolation of temperature across the entire network of Iceland weather stations (Crochet and Jóhannesson, 2011). This product is based on adjusting the weather station temperatures to equivalent sea-level temperatures using the known heights of the stations and a saturated adiabatic lapse rate (SALR) of 6.5°C km⁻¹, and then interpolating across a 1 km grid of Iceland using the grid elevations and the SALR. This 1 km dataset is within 1°C of 60–80%, and within 2°C of 90–95% of independent station measurements in 1995–2010, depending on the month of the year (Crochet and Jóhannesson, 2011).

We downscaled the 1 km temperature data across Langjökull to the 30 m grid of our DEMs in two steps: (1) the 1 km grid was re-gridded to 30 m using bilinear interpolation; (2) for each 1 km grid cell, the elevation difference between the 1 km cell and each of the 30 m cells contained within it was determined, and the temperature was adjusted from the 1 km data using the 6.5°C km⁻¹ SALR used by Crochet and Jóhannesson (2011).

4.4. Precipitation data

Precipitation data are also needed to drive a distributed surface mass-balance model for Langjökull. Like the temperature data, precipitation measured at Hveravellir has been shown to have a complex and variable relationship with precipitation and accumulation measured across Langjökull (Björnsson and others, 2002). Thus, we use a 1 km gridded product of daily precipitation derived from ERA-40 reanalysis data and a linear model of orographic precipitation (Crochet and others, 2007). This product has been validated against accumulation measurements on Langjökull from 1997 to 2002 and shows good agreement: R ² = 0.4, slope = 0.95, intercept = 1.14 mm, mean relative error = 1.3%, median relative error = −4.9%, RMSE = 2.9 mm, mean error = −0.7 mm and median error = −0.3 mm. The 1 km precipitation data are downscaled to the 30 m grid of our DEMs using bilinear interpolation.

Unfortunately, the precipitation product is only available up to January 2007, and so for February–August 2007, another means of generating 30 m precipitation data was necessary. We use, therefore, the daily data collected at Hvervellir. Neither a precipitation lapse rate nor a linear interpolation between the Hvervellir data and data from other Icelandic weather stations are appropriate means of generating the precipitation field across Langjökull and so we use instead a quantile-mapping technique (e.g. Panofsky and Brier, 1968; Maurer and others, 2010; Rye and others, 2010). The daily precipitation measurements at Hveravellir from January 1996 to January 2007 were broken into 1000 bins. Similarly, for each 1 km grid cell across Langjökull, the daily precipitation data from the Crochet and others (2007) model were also broken into 1000 bins. Lookup tables were produced to extrapolate daily precipitation amounts at Hveravellir to daily precipitation in each 1 km grid across Langjökull. The daily precipitation amounts across Langjökull for February–August 2007 were then derived using the appropriate lookup tables from the February–August 2007 Hveravellir data. As above, the 1 km scale precipitation data were downscaled to 30 m resolution using bilinear interpolation.

5. METHODS

5.1. Elevation change

Elevation change is calculated across two epochs, 1997–2004 and 2004–07. Following previous naming conventions, H indicates an elevation measurement at a given time and dH is a measured change in height over an epoch. Because coincident on-ice and off-ice data were not available across all three time periods, it was not possible to develop a model that considers both registration and slope bias (e.g. Kohler and others, 2007; Rees and Arnold, 2007; Nuth and others, 2012). Although we do not expect such errors to be large, it is possible that some may persist despite the common 30 m grid. Such errors will not significantly influence volume change calculations for the entire ice cap, although the possibility of such errors must be considered when investigating elevation changes at individual points.

In elevation change calculations, errors will be both random and spatially autocorrelated, and these are taken into account using a geostatistical model (Rolstad and others, 2009). We assume a spherical semivariogram (as in Barrand and others, 2010) and a correlation length of 740 m (as suggested by Rolstad and others, 2009). Thus, the correlated area is a circle with radius of 740 m, and the accuracy of elevation change is:

(1) $${E_{{\rm geostat}}} = {\left( {\displaystyle{{0.2 E_{{\rm dH}}^2 \times {A_{{\rm cor}}}} \over {{A_{{\rm total}}}}}} \right)^{0.5}}$$

where E _geostat is the error taking into account a geostatistical model (i.e. spatial autocorrelation), A _cor is the area of spatial autocorrelation (in this case, πr ² ≈ 1.72 km²), A _total is the total surface area of the ice cap. E _dh is the error in elevation difference at each independent point given by:

(2) $${E_{{\rm dH}}} = {\left( {E_{{\rm H}1}^2 + E_{{\rm H}2}^2} \right)^{0.5}}$$

where E _H1 and E _H2 indicate the uncertainty in elevation of the DEMs at the beginning and end of the epoch. The uncertainty used is specific to each DEM, i.e. 2, 1 and 2 m for 1997, 2004 and 2007, respectively.

5.2. Surface mass-balance model

Surface ablation is calculated at a daily time step using a simple degree-day model (Hock, 2005). For all daily temperatures, T (°C) above zero, melt, M (mm) is modelled as:

(3) $$M = DDF \times T$$

where DDF is a degree-day factor specific to either snow or ice; the surface used is based upon a snow mask initialized with the first year's precipitation and subsequent tracking of the relative accumulation and ablation at each cell. We use values of 5.8 and 7.0 mm °C⁻¹ respectively, since these have been optimized in previous degree-day modelling studies on Langjökull (Guðmundsson and others, 2003). They are comparable with similarly optimized values for nearby Hofsjökull, ~30 km east of Langjökull (5.3 and 7.3 mm °C⁻¹, respectively) (Aðalgeirsdóttir and others, 2006) and close to observed lapse rates for Langjökull (6.0 mm °C⁻¹, Guðmundsson and others, 2009; 5.8 mm °C⁻¹, Hodgkins and others, 2013). Following earlier studies, a refreezing ratio of melt on snow of 0.07 is also included (Jóhannesson and others, 1995; Aðalgeirsdóttir and others, 2006). Surface accumulation is driven by the gridded precipitation field as described above with a precipitation threshold above/below which rain/snow falls, treated as a tuneable parameter varying from 0°C to 3°C at 0.5°C increments (cf. Jóhannesson and others, 1995; De Woul and others, 2006).

The model was run from September 1996 (the beginning of the 1996/97 mass-balance year) to August 2007 (the end of the 2006/07 year) with gridded output generated for the dates of each DEM collection, 20 April 1997, 15 August 2004 and 2 August 2007. This allows for surface mass balance differencing coincident with the epochs between the DEMs.

5.3. Dynamic change

Surface elevation change with time (dH/dt) is the result of contributions from multiple sources:

(4) $$\displaystyle{{{\rm d}H} \over {{\rm d}t}} = \displaystyle{A \over {{\rho _{{\rm sf}}}}} - {V_{{\rm fc}}} - \displaystyle{{{A_{\rm b}}} \over {{\rho _{\rm i}}}} - {V_{{\rm ice}}} + \displaystyle{{{\rm d}B} \over {{\rm d}t}}$$

where A is the snow accumulation rate, ρ _sf is the density of near-surface firn, V _fc is the velocity of firn compaction at the surface, A _b is the ablation rate, ρ _i is the density of ice, V _ice is the vertical velocity of the ice at the firn/ice transition and dB/dt is the vertical bedrock motion (e.g. Li and Zwally, 2011). As commonly defined, parameters H, A and dB/dt are positive upward and V _fc, A _b and V _ice are positive downward.

We ignore dB/dt, which is currently only ~0.02 m a⁻¹ measured at GPS stations around Langjökull (Árnadóttir and others, 2009; Geirsson and others, 2010; Compton and others, 2015). We also remove V _fc, which is only relevant in the accumulation area. There are very few temperate firn compaction models, and the surface processes they describe are quite complex, such as water percolation and refreezing (Vimeux and others, 2009; Huss, 2013). Thus, we do not attempt to model firn compaction explicitly. Since the epochs are at least 3 a long, and because the climate is largely stable over the complete time period (according to weather station data), we expect V _fc to be small. Removing these two terms and rearranging Eqn (4), we calculate V _ice from:

(5) $${V_{{\rm ice}}} = \displaystyle{A \over {{\rho _{{\rm sf}}}}} - \displaystyle{{{A_{\rm b}}} \over {{\rho _{\rm i}}}} - \displaystyle{{{\rm d}H} \over {{\rm d}t}}.$$

Since the mass balance model outputs A and A _b are in m w.e., we divide the net surface balance calculations by ρ _sf, where our model predicts net accumulation and by ρ _I, where it predicts net ablation. This is largely valid for end-of-seasons displacements, where the net mass loss is assumed to be ice and the net gain to be firn. However, for a spring-to-summer comparison (as in the first epoch), we modify Eqn (5) to account for the removal of seasonal snow:

(6) $${V_{{\rm ice}}} = \displaystyle{A \over {{\rho _{{\rm sf}}}}} - \displaystyle{{{A_{{\rm spring}}}} \over {{\rho _{{\rm snow}}}}} - \displaystyle{{{A_{{\rm b\; net}}}} \over {{\rho _{\rm i}}}} - \displaystyle{{{\rm d}H} \over {{\rm d}t}}$$

where A _spring is springtime snow accumulation as modelled on the day of DEM data collection and A _{b net} is net ablation after removal of that springtime snow (i.e. A _{b net} = A _b − A _spring, again with both accumulation and ablation treated as positive quantities). Where there is springtime ablation (i.e. no net accumulation), we revert to using Eqn (5).

We use recommended values of ρ _i = 850 ± 60 kg m⁻³ (Huss, 2013), ρ _sf = 600 ± 100 kg m⁻³ and ρ _snow = 400 ± 100 kg m⁻³ (Cuffey and Patterson, 2010; Huss, 2013). V _ice will be positive (i.e. downward) where dH/dt ≤ (A/ρ _sf − A _b/ρ _i). This will typically occur in the accumulation area, where the downward velocity is known as a submergence velocity. Conversely, V _ice will be negative (i.e. upward) where dH/dt ≥ (A/ρ _sf − A _b/ρ _i). This will typically occur in the ablation area, where the upward velocity is known as an emergence velocity.

6. RESULTS AND DISCUSSION

6.1. Elevation change

Elevation change in the two epochs, 1997–2004 and 2004–07, is shown in Figure 2, and the mean elevation changes for the main outlet glacier catchments and for the entire ice cap are shown in Table 1. Between 1997 and 2004, the dominant trend across the ice cap is one of surface lowering (Fig. 2a), with an average height change of −1.83 ± 0.01 m a⁻¹ (Table 1). More specifically, the area below 1100 m (ablation area) changed elevation by −1.99 ± 0.01 m a⁻¹ and the accumulation above 1100 m by −1.73 ± 0.01 m a⁻¹. The surface lowering in the ablation area was particularly marked for southern outlet glaciers and less so for those in the north. This thinning also includes the removal of seasonal snow cover in 1997, which although averaged across the entire epoch of over 7 a, can be up to 7 m in this region (Berthier and others, 2014) and would also be expected to have a stronger influence on the southern outlet glaciers and lower elevations.

Fig. 2. Elevation change derived from DEM differencing: (a) 1997–2004; (b) 2004–2007; and (c) elevation change as a function of elevation for both epochs.

Table 1. Elevation change and modelled mass balance across Langjökull and its outlet glaciers (listed from north, proceeding clockwise)

Outlet names in quotes are not officially accepted but are used for convenience based on nearby geographical features.

The obvious exception to the general trend is the surface height gain along the southern and southeastern margins of Eystri-Hagafellsjökull and ‘Jarlhettur.’ Here, surface elevation rose by up to 10 m between 1997 and 2004 across an area of ~16 km²; this is associated with the 1998 surge of this outlet glacier. If anything, we expect the surging signal to have been dampened by the inclusion of seasonal snow in the first DEM collection, and it is, nevertheless, quite strong. In addition, the magnitude of the surface elevation gain and the area experiencing the gain would likely have been larger had measurements been made immediately following the surge.

By contrast, between 2004 and 2007, the areas above ~1100 m typically experienced an elevation gain, while areas below ~1100 m largely experienced elevation loss (Fig. 2b). The average height change across the ice cap was therefore lower than in the earlier epoch at −0.41 ± 0.01 m a⁻¹ (Table 1). The area below 1100 m lost elevation (−1.55 ± 0.02 m a⁻¹), while the area above 1100 m gained elevation (0.21 ± 0.02 m a⁻¹). As before, some of the largest surface height changes are associated with the 1998 surge of Eystri-Hagafellsjökull as the outlet enters a renewed phase of quiescence and retreat, possibly building towards a future surge. In this later epoch, the largest elevation gains on the ice cap occurred in the lower accumulation area of Eystri-Hagafellsjökull, whereas some of the greatest elevation drops occurred in the glacier's ablation area, most notably in the upper part, in the west, and at the margin.

The effect of the surge can also be seen by plotting elevation change against elevation (Fig. 2c). Across the first epoch (surge, 1997–2004), the ablation zone is expected to thin due to seasonal snow loss but thicken from the surge itself. The accumulation zone, however, is expected to thin due to both seasonal snow cover and the surge. Accordingly, Figure 2c shows some thickening at the lowest elevations and increased thinning at higher elevations. By contrast, across the second epoch (recovery, 2004–07), the ablation zone is expected to strongly thin due to decreased mass flux and negative surface mass balance, while the accumulation zone increases should thicken as ice flux has waned. As expected, Figure 2c shows strong thinning at low elevations transitioning to modest thickening at high elevations (especially above ~1200 m). By comparing elevation change with elevation itself, despite the fact that the surge only majorly impacts one outlet, we see that the effects of the 1998 surge of Eystri-Hagafellsjökull and subsequent quiescence dominate the elevation change signal.

6.2. Surface mass balance

As mentioned above, we fix the parameters of temperature lapse rate, DDFs for snow and ice, and the refreezing ratio, as these have been relatively well constrained in previous modelling studies for Langjökull. We run the model seven times for the different temperature thresholds (for snow vs rain) and evaluate model performance against the measured mass balance data (Table 2). In particular, we evaluate it against: (1) the RMSEs of the modelled vs measured winter, summer and net surface mass balance for the individual stake data (b _w, b _s and b _n, respectively); (2) the coefficients of determination (R ² values) of regression equations fitted to modelled vs measured winter, summer and net mass balances at the individual stakes; (3) the RMSEs of the modelled vs measured winter, summer and net surface mass balance averaged across the ice cap (B _w, B _s and B _n, respectively); (4) the R ² values of regression equations fitted to modelled vs measured winter, summer and net mass balances averaged across the ice cap; (5) differences between icecap wide modelled and measured cumulative winter, summer and net mass balance at the end of the 10 a.

Table 2. Model evaluation against measurements for different precipitation thresholds (0–3°C), above/below which precipitation is treated as rain/snow

b _w, b _s and b _n are winter, summer and net stake mass balance, respectively. B _w, B _s and B _n are winter, summer and net ice cap wide mass balance, respectively. ‘Cum Diff’ is the difference between the cumulative mass balances after 10 a (m w.e.). The optimal threshold evaluated against the different criteria is highlighted in bold.

The results are shown in Table 2. Taking all measures into account, the optimal results are produced for a precipitation threshold of 2°C and hereafter model results are shown and discussed only for this run. We arrived at this decision as follows: for each evaluation criterion, we ranked the criteria across the precipitation thresholds (where 1 is the ‘best’, i.e. minimum RMSE or maximum R ² and 7 is the ‘worst’). We summed the ranks and chose 2°C as it had the lowest sum of ranks (Table 2).

Figure 3 shows the modelled outputs together with the stake measurements made between September 1996 and August 2007. For b _n, the model predicts small negative balances at a few stakes in some years where the measurements suggest small positive balances; other than this, the modelled b _n agree well with the measurements, but with a degree of scatter (Fig. 3a). It is possible that local snow drifting could explain this, but stakes are located to minimize such effects. For b _s, the model successfully predicts net ablation at all stakes and the overall fit is good, but with a tendency to under predict ablation by several metres at some stakes in some years (Fig. 3b). For b _n, the fit is generally good across all stakes in both the accumulation and ablation areas (Fig. 3c). Overall, the model successfully reproduces the magnitudes of b _w, b _s and b _n.

Fig. 3. Modelled versus measured surface mass balance (b) at the stakes for the mass balance years 1996/97–2006/07: (a) winter, (b) summer and (c) net. Also shown are the 1:1 line (solid) and the regression line (dashed). ‘+’ symbols indicate stakes below 1100 m and ‘x’ symbols indicate stakes above 1100 m; the red symbols indicate the ten largest outliers in each figure, further detailed in Table 3.

The ten biggest outliers for the winter balance all lie above the respective regression line showing that there are a few stakes in a few years where much more snowfall is measured than predicted (Fig. 3a). Conversely, the ten biggest outliers for the summer balance all lie below the respective regression line (Fig. 3b) showing that there are places and times where much more melt is observed than predicted. The ten biggest outliers for the net balance lay either side of the line, seven below (measured net balance is more negative than modelled) and three above (measured net balance is less negative or more positive than modelled) (Fig. 3c).

To detect any spatial or temporal patterns in these extreme biases, the stakes and dates of the 10 biggest outliers are shown in Figure 3 and Table 3. Biases in modelled b _w occur low down on Vestari-Hagafellsjökull (L01), in the col to the north east of Geitlandsjökull (L10), and along the high points of the ice cap (L07, L16 and L18) (Fig. 3a; Table 3; Fig. 1 for stake locations). The precipitation patterns across Langjökull reflect the passage of low pressure cyclones from the south-west. L10 is in the lee of Geitlandsjökull and may receive more snowfall than modelled due to the forced ascent and high deposition there (Fig. 1). L06, L07, L16 and L18 high on the ice cap may receive higher than modelled snowfall due to the precipitation model occasionally underestimating the orographic enhancement of precipitation at high elevations (field observations support this theory), or it may also be due to the use of a precipitation threshold for snow that is too low at these elevations. We use a threshold of 2°C as this best matches the overall winter, summer and net balances, although b _w alone is best modelled with a threshold of 3°C (Table 2). L01 may receive more snowfall than modelled as it is the first part of the ice cap to be reached by precipitation-bearing cyclones from the southwest and may therefore receive greater than average precipitation for its elevation and position.

Table 3. Summary of the 10 largest outliers in bias between measured and modelled stake mass balance

Stake locations are shown in Figure 1.

Biases in modelled b _s, with one exception, all occur at low elevations around the ice cap (L01 and L02 in the south, L11 in the west and L20 and L22 in the northwest and north east respectively; Fig. 3b; Table 3; Fig. 1). This may be due to the lapse rate of 6.5°C km⁻¹ being too low at these low elevations. Air temperatures will likely be enhanced around the margin due to enhanced heating of air above the fore field and the advection of the air onto the ice cap. Furthermore, the constant DDF for ice (7.0 mm °C⁻¹) may underestimate melt around the ice cap margin where albedos are noticeably lower than at higher elevations within the ablation area. Multiple studies have indeed shown that melt modelling at Langjökull is sensitive to variable lapse rates and DDFs (Guðmundsson and others, 2009; Hodgkins and others, 2012, 2013; Matthews and others, 2015).

Most of the biases in modelled b _n are due to the biases in either modelled b _s or b _w (Fig. 3; Table 3). These are L01 in 1999 and L06 in 2003, where b _n is less negative than modelled, due to the higher measured than modelled winter accumulation; and L01 in 2007, L11 in 2002 and 2007, and L20 and L22 in 1997, where b _n is more negative than modelled, due to the higher measured than modelled summer ablation. Some of the biases in modelled b _n are due to a combination of moderate biases in both the winter and summer balances (Fig. 3c; Table 3). These are L11 in 2000 and L15 in 2003, where b _n is more negative than modelled, and L15 in 2000, where b _n is more positive than modelled.

Figure 4 shows the cumulative modelled B _w, B _s and B _n for the entire ice cap between 1997 and 2007 compared with those derived by interpolation/extrapolation from the stake measurements. It confirms that modelled winter accumulation is slightly underestimated and summer ablation is slightly overestimated with respect to the measured balances (Figs 4a, b). The result is a systematic cumulative overestimation of modelled negative net mass balance compared with the stake-derived values, although by the end of the 10 a period, the cumulative B _n converges (Fig. 4c). Across the 10 a, modelled and measured B _w, B _s and B _n are in good agreement given the errors of both techniques.

Fig. 4. Modelled vs measured cumulative ice cap wide surface mass balance (B), 1996/97–2006/07: (a) winter, (b) summer and (c) net.

Patterns of net mass balance across the ice cap for the two epochs are shown in Figure 5. Patterns are similar for the two epochs with strong vertical gradients ranging from high negative values of up to −9.2 m w.e. a⁻¹ at the termini of the outlet glaciers in the south to small positive values of up to +1.4 m w.e. a⁻¹ at high elevations in the interior. There is also a significant south-east to north-west gradient with the tongues of both Eystri- and Vestari-Hagafellsjökull to the south, Leiðarjökull in the east, and Suðurjökull in the south-east having more negative mass balances than the tongues of glaciers at similar elevations further north and west, for example, Þrístapajökull. Although patterns are similar between the two epochs, magnitudes of net balance vary. Net balance is everywhere less positive (e.g. accumulation area) or more negative (e.g. ablation area) in 1997–2004 than in 2004–07, which can be seen both as a climate signal as well resulting from the timing of the epochs (i.e. the earlier epoch has a larger proportion of melt season). The net balance for the entire ice cap and for each outlet glacier is given in Table 1, confirming that in both epochs, the catchments with the most negative balances are in the east and south, especially ‘Jarlhettur’ and ‘Skálpanes’, while those with the least negative balance are in the north and west, including Baldjökull, Flosakarðsjökular and Geitlandsjökull. The net balance for the entire ice cap, using an uncertainty estimated at 15%, is −1.46 ± 0.22 m w.e. a⁻¹ for 1997–2004 and −0.86 ± 0.13 m w.e. a⁻¹ for 2004–07.

Fig. 5. Modelled ice cap wide net surface mass balance: (a) 1997–2004 and (b) 2004–07.

6.3. Ice dynamics

At each grid cell, the difference between the surface mass balance (Fig. 5) corrected for density and the surface height change (Fig. 2) is the vertical component of ice velocity (Eqns (5) and (6)) and is shown in Figure 6. As discussed above, the surface mass balance is first divided by snow, surface firn or ice density (for springtime accumulation, net accumulation at the end of the year and ablation areas, respectively) before dividing by the length of the epoch to convert to m a⁻¹. The 1997–2004 epoch shows predominantly submergence velocities in the accumulation area above ~1100 m elevation (lower in some places) and emergence velocities in the ablation area. The 2004–07 epoch shows a similar pattern, but shifted to slightly higher elevations, with submergence velocities confined to the upper accumulation area above ~1200 m elevation.

Fig. 6. Dynamic component of elevation change: (a) 1997–2004 and (b) 2004–07.

Showing by far the biggest difference between epochs is the Eystri-Hagafellsjökull and ‘Jarlhettur’ catchments, due to the 1998 surge (Fig. 7). In the 1997–2004 epoch, there is a very large submergence velocity zone, in which velocities are particularly high between the elevations of ~800–1100 m (Fig. 8a, labelled A). Conversely, there is a zone of particularly high emergence velocities at the terminus and along the eastern margin, below elevations of ~800 m (Fig. 8a, labelled B). Less obvious is a small area of emergence velocities at elevations near 900 m at the western margin where the ice comes up against the bedrock spur, Hagafell, separating Vestari- and Eystri-Hagafellsjökull (Fig. 8, labelled C). The patterns show the effects of the 1998 surge and the transfer of mass and associated extension and vertical subsidence from a reservoir area between ~800 and 1100 m to the receiving area at lower elevations, where enhanced compression and vertical uplift occur.

Fig. 7. Difference in the dynamic component of elevation change between the 2004–2007 epoch and the 1997–2004 epoch. Positive areas (in red) indicate decreased emergence (or increased submergence) in the more recent epoch, while negative areas (in blue) indicate decreased submergence (or increased emergence) in the more recent epoch.

Fig. 8. Dynamic component of elevation change showing the effect of the 1998 surge on Vestari-Hagafellsjökull: (a) 1997–2004 and (b) 2004–07.

In the 2004–07 epoch, the vertical velocity patterns across Eystri-Hagafellsjökull are more similar to those of the non-surge glaciers with submergence velocities above ~1200 m and emergence velocities at lower elevations. The key exception is the zone of higher emergence velocities at ~1100–1200 m, which we interpret as a zone of horizontal compression, vertical uplift and the build up of mass at the top of the reservoir area during the early stages of quiescence (Fig. 7b, labelled D). Conversely, we interpret the small area of high submergence velocities at the western margin, where the ice comes up against the bedrock spur, Hagafell, as a zone of vertical subsidence associated with the wastage and deformation of ice that has become isolated from the main flow of Eystri-Hagafellsjökull during the quiescent phase (Fig. 7b, labelled E). A similar pattern appears to dominate the terminus of Eystri-Hagafellsjökull and ‘Jarlhettur,’ although it is possible this is a data artefact.

The effects of the 1998 surge of the upper part of Vestari-Hagafellsjökull that failed to propagate beyond ~4 km from the terminus (Björnsson and others, 2003) are not apparent in the vertical velocity field for 1997–2004 (Fig. 6a), although it could perhaps be argued that there appears to be some emergence at similar elevations to its neighbouring surged outlet. This suggests that the effects of the surge on the vertical velocity field were relatively small, and were subsumed by method uncertainties and the effects of non-surge activity in the subsequent years.

Björnsson and others (2003) stated that although the surges of the southern outlets of Langjökull are well documented, there has been no reported surging of the steep western and eastern outlets. However, subsequent analysis of horizontal velocity patterns derived from interferometric synthetic aperture radar and dGPS led Palmer and others (2009) to conclude that Suðurjökull is a surge-type outlet, and experienced a surge between 1999 and 2004. Our data support this, showing that compared with the rest of Langjökull's outlets, Suðurjökull's terminus experienced high emergence velocities during the 1997–2004 epoch of up to 11 m a⁻¹ (Fig. 6a). This outlet glacier experienced more typical submergence and emergence velocities during the 2004–07 epoch (Fig. 6b). It appears that the neighbouring Norðurjökull also experienced sustained high emergence velocities (~8–10 m a⁻¹) in the earlier epoch (Figs 6, 7). Although no known surges have been reported for this outlet glacier, perhaps this is evidence that this outlet also experiences periods of slower and faster movement. It is also possible that the small, more northerly outlets of the Flosakarðsjöklar show evidence of increased emergence in the more recent epoch. Due to their small size and lack of notable advance, this also may be due to DEM artefacts (Fig. 7).

It is important to consider potential limitations of this method and therefore caution some interpretations. We have good confidence in the DEMs used to calculate elevation change, but the density corrections that have been applied have high uncertainties. Similarly, the assumption that firn compaction is not significant for Langjökull also contributes to high uncertainties. These assumptions were applied based on a general understanding of snow properties and compaction, in the absence of specific data from Langjökull in particular. Improvements in vertical velocity fields can thus come from a more complete understanding of snow compaction on temperate glaciers. The density assumption for ice is more robust, but the study would have been improved with two equal length epochs rather than one short and one long. This would also increase the signal-to-noise ratio for understanding surge behaviour.

In addition to the surface elevation change uncertainties, errors also come from the modelled surface mass balance. Some of these derive from the model itself (i.e. using a degree day model rather than a surface energy balance model), but this was justified by the lack of detailed climate data in the region. More complete, validated temperature and precipitation data will also remove uncertainty and the inherent dependency of mass balance on elevation. Ultimately, these errors are not expected to be systematic, and therefore this study's results are still valid, but lower uncertainty in surface mass balance would make the method more robust and increase confidence in interpreting small vertical velocity signals.

6.4. Mass balance comparisons: glaciological, geodetic and modelled

Mass balance is measured in many different units; to be able to compare glaciological, geodetic and modelled mass balance these must be reconciled. Above, elevation change measured between DEM collection dates is discussed. To convert volume to mass when integrating change over large areas (the entire ice cap and the individual outlet glaciers), the average density of the material gained or lost (ice, firn, snow) needs to be known.

For 2004–07, because both DEMs were collected in late summer, it is reasonable to assume steady-state firn compaction, negligible vertical bedrock motion and a density of 850 ± 60 kg m⁻³ (i.e. Sorge's Law; Huss, 2013). Thus, we can calculate the geodetic mass balance for the entire ice cap, of −0.35 ± 0.03 m w.e. a⁻¹. The glaciological mass balance (B _n) for the same time period is −1.13 ± 0.17 m w.e. a⁻¹, while modelled B _n falls in between −0.86 ± 0.13 m w.e. a⁻¹. It appears that the geodetic balance is therefore less negative than the other methods. However, the epoch is fairly short and so the uncertainty in density (Huss, 2013) may be insufficient.

The 1997–2004 epoch is more complex because the 1997 DEM was collected at the end of spring and the 2004 DEM was collected in late summer. Pálsson and others (2012) used the measured summer balance from 1997 to bring the two DEMs into line with each other, yielding −1.27 ± 0.15 m w.e. a⁻¹ (note that we use a density assumption 850 ± 60 kg m⁻³ (Huss, 2013) while Pálsson and others (2012) used 900 kg m⁻³). Using the same correction, the adjusted modelled balance is −1.07 ± 0.16 m w.e. a⁻¹. The glaciological mass balance (B _n) for the same time period is −1.36 ± 0.20 m w.e. a⁻¹. This epoch shows much stronger agreement between methods, agreeing with Pálsson and others’ (2012) conclusion that there does not, at least from this limited dataset, appear to be a bias between measurement methods.

For 1997–2004, our modelled mass balance is −1.07 ± 0.16 m w.e. a⁻¹ but the geodetic mass balance is −1.27 ± 0.15 m w.e. a⁻¹. Thus, our modelled surface mass balance is less negative than that derived from the surface elevation change, although the two are comparable given the uncertainties. The model under-predicts geodetic loss estimation by 16%. Conversely, during the 2004–07 epoch, the modelled mass balance is −0.86 ± 0.13 m w.e. a⁻¹ but the geodetic mass balance is −0.35 ± 0.03 m w.e. a⁻¹. In this epoch, the modelled surface mass balance is more negative than that derived from the surface elevation change. While the defined error bars do not overlap, it is possible that the short epoch may render uncertainties taken from Huss (2013) to be insufficient. Nevertheless, the model over-predicts the geodetic loss by 1.45 times. However, over the entire period 1997–2007, modelled surface mass balance is −1.01 ± 0.15 m w.e. a⁻¹ and the geodetic balance is −0.99 ± 0.11 m w.e. a⁻¹. Thus, across the whole time period, agreement is very good and there is minimal bias between methods.

Several factors may explain the difference between the modelled and geodetic mass balance. The under-prediction of mass loss for the 1997–2004 epoch may be due to overestimated precipitation inputs (including snowfall, due to the choice of precipitation threshold), underestimated air temperatures resulting from an incorrect SALR, or underestimated DDFs. Conversely, the over-prediction of mass loss for the 2004–07 epoch may be due to underestimated precipitation inputs (and/or associated precipitation threshold), overestimated air temperatures resulting from an incorrect SALR, or overestimated DDFs. In this study, we used the best available knowledge and datasets from published studies for driving our model. We also used locally derived parameters and physically downscaled temperature and precipitation fields. Effectively, this study necessitated a trade off between resolution and local accuracy in lapse rate. The highly local and temporal variability in all of these factors makes them very difficult to diagnose, in particular due to the limited datasets available. Ultimately, that the model is capable of both over- and under-predicting mass loss compared with the geodetic calculations suggests the variability associated with the model is random rather than systematic.

In a study similar to this one, Tennant and others (2012) compared geodetic balances with modelled mass balances over two glaciers for various epochs between 1949 and 2009. Similar to our results for Langjökull, over the entire period the model over-predicted volume loss, but within particular epochs the model both over- and under-predicted volume loss compared with the geodetic calculations. Tennant and others (2012) discussed possible reasons for these discrepancies, including variable debris cover, the assumption of static ice, inappropriate or variable DDFs, and ‘errors inherent in the low-resolution dynamical downscaling.’ Again, these results match our results for Langjökull, and future studies should directly address these uncertainties and those stated above. Apart from the Tennant and others (2012) study and our study, there have been relatively few such studies comparing geodetic and modelled mass balances for glaciers and ice caps.

However, there have been many more studies comparing geodetic balances (GeB) with direct glaciological balances (GlB). Cogley (2009) compared 105 coincident GeB and GlB measurements covering 29 glaciers and found no systematic difference between the two types of measurement, with an average difference of just −0.07 m w.e. a⁻¹ (GeB slightly more negative than GlB). The spread was large, however, showing that for particular glaciers in particular years the GeB estimate was larger or smaller than the GlB. The RMSE of the coincident measurements was 0.38 m w.e. a⁻¹ although only 3% of the differences, all with the GeB measurements more negative than the GlB measurements, were statistically significant. Similarly, Zemp and others (2013) compared ~50 coincident GeB and GlB measurements from 12 glaciers after first correcting for biases in each measurement and generic differences between them. They too found a small average difference of just 0.12 m w.e. a⁻¹ (again with GeB slightly more negative than GlB) and an RMSE of 0.23 m w.e. a⁻¹. Conversely, in this study, the GlB (−1.29 ± 0.19 m w.e. a⁻¹) is more negative than the GeB (−0.99 ± 0.11 m w.e. a⁻¹), although just barely within the error bounds. This is consistent with other studies and adds another data point to our understanding of comparing GlB with GeB.

7. CONCLUSIONS

In this study, we used in situ mass-balance survey data series, multiple DEMs and a surface mass-balance model to investigate the dynamic component of elevation change of Langjökull, Iceland from 1997 through 2007. The largest uncertainties were introduced by DEMs being created at different times of year (i.e. April vs August) and the consequent density and surface mass-balance corrections, which needed to be applied. Nevertheless, the mass-balance model successfully reproduced Langjökull's surface mass balance. Therefore, the maps of emergence and submergence velocity highlighted the immediate effects of a surge of one of Langjökull's outlets in 1998 (i.e. high elevation submergence and low elevation emergence), as well as the return to a quiescent state (i.e. refilling of a source zone, wastage of surged extent). We were thus able to document dynamic surge behaviour in the absence of in situ measurements during the surge. In addition, glaciological, geodetic and modelled mass balances were in good agreement over the whole ice cap and over the full time period, giving confidence in the study's results.

Ultimately, this study demonstrates a method for deriving vertical velocity fields for surge and non-surge glaciers. This enables spatial dimensions of source and sink areas to be quantified. This information could be used to better understand the basal conditions that control glacier surging behaviour. Improvements to the method will come from using a more detailed surface energy balance model to calculate surface mass balance, as well as using a model of surface density changes and firn compaction, calibrated against in-situ measurements.