1. INTRODUCTION
Glaciers and ice caps are important components of the world's hydrological cycle. Their mass loss has contributed 226 ± 135 Gt a−1 (0.62 ± 0.37 mm sea level equivalent (s.l.e.) a−1) from 1971 to 2009, and 275 ± 135 Gt a−1 (0.76 ± 0.37 mm s.l.e. a−1) from 1993 to 2009, to the global oceans (IPCC, Reference Stocker2013). Monitoring their changes is key to understanding the impacts of global and regional climate change (e.g. Oerlemans, Reference Oerlemans1994), the implications for global sea-level rise (e.g. Radić and Hock, Reference Radić and Hock2011; Marzeion and others, Reference Marzeion, Jarosch and Hofer2012; Chen and others, Reference Chen, Wilson and Tapley2013) and the effects on regional and local hydrology, including flooding (Dahlke and others, Reference Dahlke, Lyon, Stedinger, Rosqvist and Jansson2012), river biodiversity (Jacobsen and others, Reference Jacobsen, Milner, Brown and Dangles2012) and water supply to large populations (Hopkinson and Demuth, Reference Hopkinson and Demuth2006; Björnsson and Pálsson, Reference Björnsson and Pálsson2008; Baraer and others, Reference Baraer2011; Barry, Reference Barry2011; Bolch and others, Reference Bolch2012).
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, Reference Cogley2011).
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, Reference Andreassen1999; Cogley, Reference Cogley2009; Zemp and others, Reference Zemp2013; Wang and others, Reference Wang, Li, Li, Wang and Yao2014). 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, Reference Geist, Elvehøy, Jackson and Stötter2005; Hagen and others, Reference Hagen, Eiken, Kohler and Melvold2005; Fischer, Reference Fischer2011).
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, Reference Hubbard2000; Nuth and others, Reference Nuth, Schuler, Kohler, Altena and Hagen2012). 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−1]) 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, Reference Cogley2011). 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, Reference Zwally2005; Pritchard and others, Reference Pritchard, Arthern, Vaughan and Edwards2009). 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, Reference Fischer2011; Nuth and others, Reference Nuth, Schuler, Kohler, Altena and Hagen2012).
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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003; McMillan and others, Reference McMillan2014). 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, Reference Jiskoot, Murray and Boyle2000; Sevestre and Benn, Reference Sevestre and Benn2015). 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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003). More recently, remote sensing imagery has been used to identify surge-type glaciers over large regions (e.g. Grant and others, Reference Grant, Stokes and Evans2009; Paul, Reference Paul2015; Sevestre and Benn, Reference Sevestre and Benn2015) and to measure velocity increases associated with particular surges (e.g. Pritchard and others, Reference Pritchard, Murray and Luckman2005; Mansell and others, Reference Mansell, Luckman and Murray2012). 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 (Reference Minchew, Simons, Hensley, Björnsson and Pálsson2015) also measured horizontal flow velocities on Langjökull using InSAR (up to 75 m a−1), 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 km2, an elevation range of ~460 to ~1440 m a.s.l., and a volume of ~190 km3, equivalent to 0.5 mm of eustatic sea-level rise (Pálsson and others, Reference Pálsson2012). 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, Reference Pope, Willis, Rees, Arnold and Pálsson2013). It has a mean ice thickness of ~210 m and a maximum of ~650 m (Björnsson and Pálsson, Reference Björnsson and Pálsson2008). 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, Reference Eyre, Payne, Baldwin and Björnsson2005). 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, Reference Palmer, Shepherd, Björnsson and Pálsson2009). 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.
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−1 over the century (Pálsson and others, Reference Pálsson2012). The ice cap is predicted to disappear completely by ~2140 for the A1B scenario of ~2°C (100 a)−1 warming (Björnsson and Pálsson, Reference Björnsson and Pálsson2008) due to low-snow accumulation and high-annual temperatures (Björnsson and others, Reference Björnsson, Pálsson and Haraldsson2002). 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, Reference Sigurðsson1998; Björnsson and others, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003; Palmer and others, Reference Palmer, Shepherd, Björnsson and Pálsson2009). During surges, ice velocities increase over areas of up to 200 km2. 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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003).
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, Reference Bennett, Huddart and Waller2005). 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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003). 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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003). No previous surges are documented for either glacier (Sigurðsson, Reference Sigurðsson1998; Bennett and others, Reference Bennett, Huddart and Waller2005; Björnsson and Pálsson, Reference Björnsson and Pálsson2008). Past surges of Eystri-Hagafellsjökull formed complex sediment deformation structures in and around lake Hagavatn (Bennett and others, Reference Bennett, Huddart and McCormick2000). In addition, Palmer and others (Reference Palmer, Shepherd, Björnsson and Pálsson2009) 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, Reference Larsen, Miller and Geirsdóttir2013). 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, Reference Pálsson2012; 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, Reference Pálsson2012; 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, Reference Pope, Willis, Rees, Arnold and Pálsson2013; 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, Reference Jóhannesson2013). 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, Reference Björnsson, Pálsson and Haraldsson2002; Pálsson and others, Reference Pálsson2012) (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, Reference Pálsson2012; 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, Reference Guðmundsson, Björnsson, Pálsson and Haraldsson2003; Flowers and others, Reference Flowers, Björnsson, Geirsdóttir, Miller and Clarke2007).
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, Reference Crochet and Jóhannesson2011). 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−1, 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, Reference Crochet and Jóhannesson2011).
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−1 SALR used by Crochet and Jóhannesson (Reference Crochet and Jóhannesson2011).
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, Reference Björnsson, Pálsson and Haraldsson2002). 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, Reference Crochet2007). This product has been validated against accumulation measurements on Langjökull from 1997 to 2002 and shows good agreement: R 2 = 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, Reference Panofsky and Brier1968; Maurer and others, Reference Maurer, Hidalgo, Das, Dettinger and Cayan2010; Rye and others, Reference Rye, Arnold, Willis and Kohler2010). 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 (Reference Crochet2007) 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, Reference Kohler2007; Rees and Arnold, Reference Rees and Arnold2007; Nuth and others, Reference Nuth, Schuler, Kohler, Altena and Hagen2012). 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, Reference Rolstad, Haug and Denby2009). We assume a spherical semivariogram (as in Barrand and others, Reference Barrand, James and Murray2010) and a correlation length of 740 m (as suggested by Rolstad and others, Reference Rolstad, Haug and Denby2009). Thus, the correlated area is a circle with radius of 740 m, and the accuracy of elevation change is:
$${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 2 ≈ 1.72 km2), 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:
$${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, Reference Hock2005). For all daily temperatures, T (°C) above zero, melt, M (mm) is modelled as:
$$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−1 respectively, since these have been optimized in previous degree-day modelling studies on Langjökull (Guðmundsson and others, Reference Guðmundsson, Björnsson, Pálsson and Haraldsson2003). They are comparable with similarly optimized values for nearby Hofsjökull, ~30 km east of Langjökull (5.3 and 7.3 mm °C−1, respectively) (Aðalgeirsdóttir and others, Reference Aðalgeirsdóttir, Jóhannesson, Pálsson, Björnsson and Sigurðsson2006) and close to observed lapse rates for Langjökull (6.0 mm °C−1, Guðmundsson and others, Reference Guðmundsson, Björnsson, Pálsson and Haraldsson2009; 5.8 mm °C−1, Hodgkins and others, Reference Hodgkins, Carr, Pálsson, Guðmundsson and Björnsson2013). Following earlier studies, a refreezing ratio of melt on snow of 0.07 is also included (Jóhannesson and others, Reference Jóhannesson, Sigurdsson, Laumann and Kennett1995; Aðalgeirsdóttir and others, Reference Aðalgeirsdóttir, Jóhannesson, Pálsson, Björnsson and Sigurðsson2006). 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, Reference Jóhannesson, Sigurdsson, Laumann and Kennett1995; De Woul and others, Reference De Woul2006).
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:
$$\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, Reference Li and Zwally2011). 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−1 measured at GPS stations around Langjökull (Árnadóttir and others, Reference Árnadóttir2009; Geirsson and others, Reference Geirsson2010; Compton and others, Reference Compton, Bennett and Hreinsdóttir2015). 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, Reference Vimeux2009; Huss, Reference Huss2013). 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:
$${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:
$${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−3 (Huss, Reference Huss2013), ρ sf = 600 ± 100 kg m−3 and ρ snow = 400 ± 100 kg m−3 (Cuffey and Patterson, Reference Cuffey and Patterson2010; Huss, Reference Huss2013). 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−1 (Table 1). More specifically, the area below 1100 m (ablation area) changed elevation by −1.99 ± 0.01 m a−1 and the accumulation above 1100 m by −1.73 ± 0.01 m a−1. 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, Reference Berthier2014) and would also be expected to have a stronger influence on the southern outlet glaciers and lower elevations.
Elevation change derived from DEM differencing: (a) 1997–2004; (b) 2004–2007; and (c) elevation change as a function of elevation for both epochs.
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 km2; 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−1 (Table 1). The area below 1100 m lost elevation (−1.55 ± 0.02 m a−1), while the area above 1100 m gained elevation (0.21 ± 0.02 m a−1). 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 2 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 2 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.
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 2 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.
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.
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−1 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−1) 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, Reference Guðmundsson, Björnsson, Pálsson and Haraldsson2009; Hodgkins and others, Reference Hodgkins, Carr, Pálsson, Guðmundsson and Björnsson2012, Reference Hodgkins, Carr, Pálsson, Guðmundsson and Björnsson2013; Matthews and others, Reference Matthews, Hodgkins, Guðmundsson, Pálsson and Björnsson2015).
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.
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−1 at the termini of the outlet glaciers in the south to small positive values of up to +1.4 m w.e. a−1 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−1 for 1997–2004 and −0.86 ± 0.13 m w.e. a−1 for 2004–07.
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−1. 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.
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.
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.
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, Reference Björnsson, Pálsson, Sigurðsson and Flowers2003) 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 (Reference Björnsson, Pálsson, Sigurðsson and Flowers2003) 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 (Reference Palmer, Shepherd, Björnsson and Pálsson2009) 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−1 (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−1) 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−3 (i.e. Sorge's Law; Huss, Reference Huss2013). Thus, we can calculate the geodetic mass balance for the entire ice cap, of −0.35 ± 0.03 m w.e. a−1. The glaciological mass balance (B n) for the same time period is −1.13 ± 0.17 m w.e. a−1, while modelled B n falls in between −0.86 ± 0.13 m w.e. a−1. 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, Reference Huss2013) 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 (Reference Pálsson2012) 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−1 (note that we use a density assumption 850 ± 60 kg m−3 (Huss, Reference Huss2013) while Pálsson and others (Reference Pálsson2012) used 900 kg m−3). Using the same correction, the adjusted modelled balance is −1.07 ± 0.16 m w.e. a−1. The glaciological mass balance (B n) for the same time period is −1.36 ± 0.20 m w.e. a−1. This epoch shows much stronger agreement between methods, agreeing with Pálsson and others’ (Reference Pálsson2012) 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−1 but the geodetic mass balance is −1.27 ± 0.15 m w.e. a−1. 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−1 but the geodetic mass balance is −0.35 ± 0.03 m w.e. a−1. 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 (Reference Huss2013) 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−1 and the geodetic balance is −0.99 ± 0.11 m w.e. a−1. 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 (Reference Tennant, Menounos, Ainslie, Shea and Jackson2012) 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 (Reference Tennant, Menounos, Ainslie, Shea and Jackson2012) 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 (Reference Tennant, Menounos, Ainslie, Shea and Jackson2012) 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 (Reference Cogley2009) 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−1 (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−1 although only 3% of the differences, all with the GeB measurements more negative than the GlB measurements, were statistically significant. Similarly, Zemp and others (Reference Zemp2013) 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−1 (again with GeB slightly more negative than GlB) and an RMSE of 0.23 m w.e. a−1. Conversely, in this study, the GlB (−1.29 ± 0.19 m w.e. a−1) is more negative than the GeB (−0.99 ± 0.11 m w.e. a−1), 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.
ACKNOWLEDGEMENTS
The in situ mass-balance survey of Langjökull was a joint effort of the Glaciology Group, Institute of Earth Sciences, University of Iceland and the National Power Company (Landsvirkjun). We thank Philippe Crochet and Tómas Jóhannesson from the Icelandic Meteorological Office for providing the gridded climate data and for useful discussions about the climatology of Langjökull. The 2007 lidar data were collected by the UK Natural Environment Research Council Airborne Research and Survey Facility (Grant IPY 07-08). Additional funding was provided by the United States National Science Foundation (Grant No. DGE-1038596), St Catharine's, St John's and Trinity Colleges and the University of Cambridge B.B. Roberts and Scandinavian Studies Funds. We thank Cameron Rye for initial help coding the mass-balance model. In addition, the authors would like to thank the creators and development teams of the freely available software used in this research, including ImageJ (http://rsb.info.nih.gov/ij/), MultiSpec (https://engineering.purdue.edu/~biehl/MultiSpec/), QGIS (http://www.qgis.org/), Plot (http://plot.micw.eu/), and Zotero (http://zotero.org/).
Open access
