Downscaled precipitation applied in modelling of mass balance and the evolution of southeast Vatnajökull, Iceland

Abstract

Simulations of the post-Little Ice Age evolution of three outlet glaciers of southeast Vatnajökull, Iceland – Skálafellsjökull, Heinabergsjökull and Fláajökull – are presented. A coupled shallow-ice-approximation ice-flow and degree-day mass-balance model is applied that is calibrated with a 14 year record of in situ mass-balance measurements. The measured mass balance cannot be realistically represented by constant horizontal and vertical precipitation gradients. High-resolution (1 km) precipitation fields, derived from downscaled orographic atmospheric circulation models of precipitation, are required to capture the spatial variation of the winter mass balance. The observed ice volume around 1890 (15–30% larger than in 2000) can be simulated with 1°C lower temperatures and a 20% reduction in the annual precipitation, relative to the reference climate period, 1980–2000. The sensitivity of each glacier’s annual balance to a change in temperature is −1.51 to −0.97 m w.e. a⁻¹ °C⁻¹ and +0.16 to +0.65 m w.e. a⁻¹ for a 10% increase in precipitation. A steady-state experiment applying a step increase in temperature of 2°C (3°C), and precipitation increase of 10%, results in a >50% (80–90%) decrease in ice volume.

Introduction

The climate in Iceland is influenced by the atmospheric circulation of the North Atlantic and the oceanic boundaries defined by the warm Irminger current and the cold East Greenland current (e.g. Einarsson, 1984; Ólafsson and others, 2007). Iceland is located in the northern part of the North Atlantic storm track, and high amounts of precipitation are delivered to the temperate glaciers in Iceland, which have mass turnover rates on the order of 1.5–3.0 m w.e. a⁻¹ (Björnsson and others, 2013). The mass-balance sensitivity of the glaciers and ice caps is in the range −3.0 to −0.6 m w.e. a⁻¹ °C⁻¹ (Guðmundsson and others, 2011; Pálsson and others, 2012; Jóhannesson and others, 2013), which is among the highest in the world (De Woul and Hock, 2005). Due to the high spatial variability of the precipitation and the relative sparseness of precipitation observations, which are mostly confined to lowland areas, numerical atmospheric models, rather than statistical methods, have been used to simulate the spatial and temporal structure of the precipitation fields (Rögnvaldsson and others, 2004, 2007; Crochet and others, 2007; Jóhannesson and others, 2007 and references therein). Undercatch of rain gauges is a known problem, especially during winter and/or strong winds (e.g. Sigurðsson, 1990; Crochet, 2007), so series of observed precipitation are often quite unreliable. Estimates of precipitation by observations of snow accumulation and runoff records in the vicinity of the area of interest are, in general, a better indicator for precipitation at high elevation and in complex terrain than conventional precipitation measurements (e.g. Rögnvaldsson and others, 2004). As precipitation observations from lowland stations generally cannot be reliably extrapolated into complex or high topography, they introduce large errors in mass-balance modelling (e.g. Stahl and others, 2006; Huss and others, 2008). However, in most mass-balance models precipitation is extrapolated from meteorological stations outside the glaciers (e.g. Andreassen and Oerlemans, 2009; Aðalgeirsdóttir and others, 2011). Mass-balance models are of varying complexity, ranging from simple temperature-index models or positive degree-day (PDD) models (e.g. Laumann and Reeh, 1993; Jóhannesson and others, 1995; Braithwaite and Zhang, 2000) to surface energy-balance models (e.g. Hock, 2005 and references therein). A number of studies have used data from distant meteorological stations to produce precipitation grids for the glaciers (e.g. Auer and others, 2007; Engelhardt and others, 2012; Huss and others, 2012; Marzeion and others, 2012; Radić and others, 2014) or data from atmospheric models providing input for mass-balance modelling (e.g. Rasmussen and others, 2007; Machguth and others, 2009; Paul and Kotlarski, 2010; Andreassen and others, 2012; Jarosch and others, 2012; Van Pelt and others, 2012).

Ice-dynamical models are based on the Stokes equations or their approximations, from zeroth-order shallow-ice approximation (SIA) to higher-order models, which include the longitudinal stress gradients (e.g. Hindmarsh, 2004; Le Meur and others, 2004). A number of comparative studies indicate that numerical models of reduced complexity and full system models give similar length and volume changes for valley glaciers (Leysinger-Vieli and Guðmundsson, 2004; Oerlemans, 2008; Lüthi, 2009). Numerical ice-flow models simulating the larger ice caps in Iceland and their response to changes in mass balance have been developed in recent years (Aðalgeirsdóttir, 2003; Flowers and others, 2005; Marshall and others, 2005). A vertically integrated SIA ice-flow model coupled with a PDD surface mass-balance model (Jóhannesson and others, 1995; Jóhannesson, 1997) has been used to simulate the evolution of ice caps in Iceland (Aðalgeirsdóttir, 2003; Aðalgeirsdóttir and others, 2006; Jóhannesson and others, 2007; Guðmundsson and others, 2009a). Air temperature measured outside the glaciers represents variations in the incoming radiation flux (better than the dampened temperature of the boundary layer above the melting glacier surface), so the ablation of the glaciers can be described with PDD models (e.g. Guðmundsson and others, 2009b). The evolution of the outlet glacier Hoffellsjökull (Fig. 1) during the whole post-Little Ice Age (LIA) period ∼1890–2010 has successfully been simulated with the same coupled model (Aðalgeirsdóttir and others, 2011).

Fig. 1. (a) Iceland in the North Atlantic Ocean and major ocean surface currents: East Greenland Current (EGC), East Iceland Current (EIC), Irminger Current (IC) and North Atlantic Current (NAC). (b) Iceland and the larger ice caps: Vatnajökull (V), Langjökull (L), Hofsjökull (H), Mýrdalsjökull (M), Öræfajökull (Ö), the plateau of Breiðabunga (BB) and the locations of the meteorological stations at Fagurhólsmýri (F) and Hólar in Hornafjörður (HH). The box outlines (c). (c) Eastern Breiðamerkurjökull, Skálafellsjökull, Heinabergsjökull, Fláajökull and Hoffellsjökull. The surface topography is from the 2010 lidar DEM with 100 m contours shown in grey, and ice divides in black. The locations of mass-balance stakes are indicated by black dots; stakes not used for the calibration of the mass-balance model are shown in parentheses. Nunataks are shown in light grey, including Snjófjall (S) and Litlafell (L) in Heinabergsjökull. Proglacial lakes are shown in blue and the road is delineated in black.

In this paper, the PDD mass-balance model coupled with the vertically integrated SIA ice-flow model, solved with the finite-element method, is used to simulate the evolution of three outlet glaciers of southeast Vatnajökull and their sensitivity to climate change. The first part of the paper describes how three different datasets are used to simulate the mass balance in the coupled model runs: (1) precipitation model based on constant vertical and horizontal precipitation gradients, using data from meteorological stations outside the glacier (e.g. Jóhannesson and others, 1995; Aðalgeirsdóttir and others, 2011), (2) constant winter mass-balance grid (manually interpolated from in situ balance measurements) revised with downscaled precipitation model output data from a numerical atmospheric model and (3) downscaled orographic precipitation from a linear model which simulates physical precipitation processes.

In the second part of the paper, we describe a series of model simulations, using the downscaled precipitation data from the linear model, with step changes in precipitation and temperature. These simulations are carried out to assess the sensitivity of the modelled glaciers to perturbations in the climate. Observations of area and volume changes from the LIA maximum around 1890 until 2010 are available for the outlet glaciers of southeast Vatnajökull (Hannesdóttir and others, 2014a), and provide validation data for the model simulations. The LIA maximum glacier geometry was successfully simulated, as well as the evolution of the glaciers during the time period when downscaled precipitation model data are available (1959–2010).

Study Area

The three outlet glaciers of southeast Vatnajökull – Skálafellsjökull, Heinabergsjökull and Fláajökull – are located in one of the warmest and wettest areas of Iceland (e.g. Ólafsson and others, 2007). These non-surging outlets are 300 m thick on average, descend from a 1500 m high plateau towards the lowlands, and their termini are at 30–60 m a.s.l. (Fig. 1; Table 1). The outlet glaciers are 23–25 km long and have an average slope of 3–4°. They range in area from 100 to 170 km² and have a volume of 30–55 km³ (Table 1). The snowline at the end of summer, which is used as a proxy for the equilibrium-line altitude (ELA), has been derived from a series of Moderate Resolution Imaging Spectroradiometer (MODIS) images during the period 2007–11 (Hannesdóttir and others, 2014a), and is in the range 900–1100 m on the three outlets. (Table 1). Proglacial lakes formed in front of the outlet glaciers around the mid-20th century. Skálafellsjökull has a long, narrow accumulation area, and funnels down an icefall, where the elevation drops by 200 m over a distance of 1 km. A small proglacial lake was formed in the 1960s and has remained of similar size (2.2 km²). Heinabergsjökull is divided into three branches by two mountains (Fig. 1), and has a long, flat ablation area. The glacier terminates in a proglacial lake, which was formed in the 1940s, has increased in size since then and is now 9.5 km². Fláajökull has a wide accumulation area and a gently sloping 4 km wide tongue that terminates in a lake that formed during or prior to the 1940s, and has area 9.7 km².

Table 1. Characteristics of the four outlet glaciers in 2010

Data

Basal topography

The three outlet glaciers were surveyed with a radio-echo sounder and differential GPS (DGPS) equipment in the period 2000–05 (Magnússon and others (2012) provide details of the method). Point measurements were carried out in the ablation area, whereas continuous profiles were surveyed in the accumulation area (Fig. 2a). The vertical accuracy of the basal topography is 5–20 m, depending on the location. A small depression below the terminus of Skálafellsjökull reaches 100 m below sea level. Heinabergsjökull has excavated a 10 km long trench, which underlies the majority of its ablation area, reaching maximum depths of 220 m below sea level. The bed of Fláajökull is 200 m below sea level in the ablation area (Fig. 2a). Previous model simulations on neighbouring Hoffellsjökull show that using a filled trench or the current basal topography, had negligible effect on the evolution simulations (Aðalgeirsdóttir and others, 2011). We thus use the measured bed DEM (unfilled trench) in all model runs.

Fig. 2. (a) Basal topography derived from radio-echo sounding measurements (shown with white lines and dots) of Skálafellsjökull, Heinabergsjökull and Fláajökull. The glacier margin in 2000 is shown. Proglacial lakes are shown in light blue, and nunataks in white. (b) Surface topography from the 2010 lidar DEM shown with 100 m contours and glacier outlines at different times from Hannesdóttir and others (2014a). The location of the MODIS-derived snowline in the period 2007–11 is shown in light blue.

Glacier area and volume changes 1890–2010

The evolution of the area and volume of the three outlet glaciers during the period 1890–2010 (Fig. 2b; Table 2) has been derived from a new glacier inventory. DEMs of the ∼1890 glacier surface were created from glacial geomorphological features and historical data from the time of the LIA maximum (Hannesdóttir and others, 2014b), whereas DEMs from 1904, 1945, 1989, 2002 and 2010 were generated from maps, aerial images, DGPS surveys and a lidar survey (Hannesdóttir and others, 2014a). The most accurate DEMs were produced with airborne lidar technology in late August–September of 2010 and 2011 (IMO and IES, 2013). The glaciers started to retreat from their terminal LIA moraines after 1890 and retreated during most of the 20th century. They halted or advanced slightly in the 1960s–80s due to lower temperatures, but the retreat accelerated after ∼2000. During the first decade of the 21st century, the glaciers experienced the highest rate of mass loss in the post-LIA period (Hannesdóttir and others, 2014a).

Table 2. Observed area (A; km²) and volume (V; km³) (from Hannesdóttir and others, 2014a) and simulated (using temperature from Hólar and LT-DP) area and volume in 1989, 2002 and 2010, both from uncoupled and coupled runs (see Fig. 11)

Mass-balance observations

Mass-balance measurements have been carried out on Vatnajökull since 1991, and the majority of the stakes are located on the northern and western part of the ice cap (Björnsson and Pálsson, 2008). Six stakes were on a profile from the terminus up to the ice divide on Breiðamerkurjökull, three stakes were measured on Hoffellsjökull, and a few survey sites were located in the accumulation area of the studied glaciers (Fig. 1; Björnsson and Pálsson, 2008). The variation of mass balance with elevation on southeast Vatnajökull is shown in Figure 3a. Ablation of up to 9mw.e.a⁻¹ is observed during summers on Breiðamerkurjökull and Hoffellsjökull, and negative winter balances at the termini (Björnsson and Pálsson, 2008). Two new stakes were added to the mass-balance survey network in 2009 in the accumulation area of Skálafellsjökull and Fláajökull. Digital mass-balance maps, for every year since 1996, have been manually interpolated (Björnsson and others (2002) provide details of the method), using the in situ mass-balance measurements and the observed mass-balance gradient on southeast Vatnajökull (Fig. 3a). The resulting maps are integrated over the whole glacier area to calculate the mean specific mass balance shown in Figure 3b for the three outlets in this study.

Fig. 3. (a) The annual mass balance (b _n) measured at stakes on Breiðamerkurjökull (circles), Skálafellsjökull and Fláajökull (crosses) and Hoffellsjökull (triangles) in the period 1996–2012. (b) Mean specific winter (b _w, blue), summer (b _s, red) and net balance (b _n, black) for Skálafellsjökull (triangles), Heinabergsjökull (circles) and Fláajökull (crosses) for the period 1996–2013. New mass-balance stakes were added to the network in 2009 (the timing is shown with a vertical line).

Temperature and precipitation data from meteorological stations

Temperature and precipitation records are available from two lowland meteorological stations south of Vatnajökull (Fig. 1; Table 3), at Fagurhólsmýri (16 m a.s.l.; 8 km south of Öræfajökull) and Hólar in Hornafjörður (16 m a.s.l.; 15 km south of Hoffellsjökull). The temperature record at Hólar is available for the period 1884–90 and since 1921 (Fig. 4), while precipitation measurements started in 1931. The temperature record at Fagurhólsmýri goes back to 1898, and precipitation has been measured from 1921 until 2008 (Fig. 4). Temperature records from other stations around Iceland were used to extend the local records back to 1860 and fill in gaps by an iterative expectation maximization algorithm (described in Aðalgeirsdóttir and others, 2011). To extend the precipitation record of Fagurhólsmýri back to 1860, Aðalgeirsdóttir and others (2011) applied a linear regression between the monthly values of temperature (Hólar) and precipitation (Fagurhólsmýri) and tuned it with the available precipitation records.

Table 3. Mean summer (June–August) temperature (T) at Hólar in Hornafjörður and annual precipitation (P) at Fagurhólsmýri averaged over the respective time periods. The temperature and precipitation during the period 1860–90 are estimated by Aðalgeirsdóttir and others (2011)

Fig. 4. (a) Mean summer (June–August) and annual temperature (grey lines) at Hólar in Hornafjörður and a 5 year running average (black lines). (b) Annual precipitation at Fagurhólsmýri. Grey boxes indicate the time period of reconstructed temperature and precipitation (for details of the reconstruction see Aðalgeirsdóttir and others, 2011). The average values of the baseline period 1980–2000 used for reference in the computations are shown with orange lines.

Downscaled precipitation data

The two different numerical precipitation datasets used in this study are based on atmospheric models. First there are data from the RÁV (Reikningar á veðri) project (Rögnvaldsson and others, 2011), including dynamically downscaled precipitation in Iceland available at 3 km resolution since 1995 and 9 km resolution since 1958. The atmospheric data were produced with the non-hydrostatic atmospheric model WRF (Weather Research Forecasting; Skamarock and others, 2008). The model takes into account atmospheric physics and dynamics, with atmospheric moisture and precipitation processes parameterized using the scheme of Thompson and others (2004). The model is forced by atmospheric analyses from the European Centre for Medium-Range Weather Forecasts (ECMWF), and is widely used in atmospheric research, including in Iceland, where it has, for example, been used to reproduce precipitation on Mýrdalsjökull, southern Iceland (Ágústsson and others, 2013).

Secondly, a linear theory (LT) model of orographic precipitation (Smith, 2003), which includes airflow dynamics, condensed water advection and downslope evaporation, has been used to create a 1 km resolution precipitation dataset for Iceland for the period 1958–2006 (Crochet and others, 2007). The model is driven by a coarse-resolution 40 year reanalysis dataset from ECMWF until 2001 and from 2002, using available ECMWF analysis. The simulated precipitation is in good agreement with precipitation observations accumulated over various timescales (from wind-loss corrected rain-gauge data and glacier mass-balance measurements), both in terms of magnitude and distribution. The results indicate that the model captures the main physical processes governing orographic generation of precipitation in the mountains of Iceland (Crochet and others, 2007). The model allows simulations of the distribution of snow accumulation on glaciers in more detail than has been possible in previous glacier mass-balance studies in Iceland (Crochet and others, 2007). A new precipitation dataset for the period 2007–10 has been made by combining rain-gauge data and monthly climatic precipitation maps derived from the LT model (Crochet and others, 2007; Jóhannesson and others, 2007) through anomaly interpolation (Crochet, 2013). The rain-gauge network used in construction of the daily precipitation maps is composed utilizing data from manual and automatic stations. The data have been quality-controlled and corrected to account for measurement errors (e.g. due to wind loss, wetting and evaporation (Crochet, 2013)). For the sake of simplicity this dataset will also be referred to as the LT model downscaled orographic precipitation (LT-DP).

Models

The mass-balance model

Our mass-balance model is a PDD (temperature index) model that has been developed for temperate glaciers in Iceland and the Nordic countries (Jóhannesson and others, 1995; Jóhannesson 1997) and has previously been used in several glaciological studies in Iceland (Aðalgeirsdóttir and others, 2006, 2011; Jóhannesson and others, 2007; Guðmundsson and others, 2009a). Glacier accumulation and ablation are computed from (1) daily precipitation and temperature values from meteorological stations outside the glacier or (2) distributed fields of daily precipitation from the LT model and temperature values from meteorological stations outside the glacier. Melting of snow and ice is computed from the number of PDDs, using different degree-day factors (amount of melting per PDD) for snow and ice and a constant snow/rain threshold (1°C) (Jóhannesson, 1997). Refreezing of meltwater and wetting of the snow by melt- or rainwater are computed as a specified fraction of the remaining snow. This is a small mass-balance component for the temperate Icelandic glaciers where the winter cold wave in the snowpack is eliminated by percolating meltwater every summer even at the highest altitudes (e.g. Björnsson and Pálsson, 2008). Refreezing at depth in older firn therefore does not occur. The main effect of refreezing and retention for Icelandic glaciers is a delay in runoff from the glacier in the early part of the melt season that has only a small effect on the annual mass balance calculated here (Jóhannesson and others (1995) provide more details).

The in situ mass-balance data from southeast Vatnajökull were used to calibrate the mass-balance model using the LT-DP. The precipitation model extends to 2010, and mass-balance measurements started in 1996, resulting in a 14 year calibration period. The mass-balance model calculates daily mass-balance values, but the model output is defined for summer (1 May–30 September) and winter (1 October–30 April) for every balance year. The degree-day scaling factors for snow (ddf_s) and ice (ddf_i), are determined by minimizing the root mean square (RMS) between measured and modelled mass balance (e.g. Jóhannesson and others, 1995). The ddf_s is derived from seven points in the accumulation area, and the ddf_i from six points in the ablation area (Fig. 1). The best fit (lowest RMS value) was found when using ddf_s =3.7 mm w.e. °C⁻¹ d⁻¹ and ddf_i = 5.5 mm w.e. °C⁻¹ d⁻¹ (Table 4). These values are in the range of previously determined degree-day factors from studies on Langjökull and Vatnajökull (Guðmundsson and others, 2003; Flowers and others, 2007; Jóhannesson and others, 2007; Aðalgeirsdóttir and others, 2011). The modelled and measured winter and summer mass balance are plotted in Figure 5. The mass-balance model explains 87% of the variance of the winter balance and 92% of the variance of the summer balance in the period 1996–2010. However, the correlation is worse when considered separately for the accumulation and ablation areas.

Table 4. Degree-day factors for snow and ice used in mass-balance modelling studies on Vatnajökull, determined by minimizing RMS between observed and modelled mass balance

Fig. 5. Measured vs modelled (using LT-DP and ddf_i = 5.5 mm w.e. C⁻¹ d⁻¹ and ddf_s = 3.7 mm w.e. °C⁻¹ d⁻¹) winter (a) and summer (b) balance of southeast Vatnajökull in the period 1996–2010. The linear relationship is shown and the R ² value indicates the linear regression coefficient.

The ice-flow model

Our ice-flow model is based on the SIA and vertically integrated continuity equation, neglecting longitudinal stress gradients, bed isostatic adjustments and seasonal variations in sliding (e.g. Aðalgeirsdóttir, 2003; Aðalgeirsdóttir and others, 2006). Glen’s flow law with an exponent n = 3 is assumed as a constitutive equation for ice (e.g. Cuffey and Paterson, 2010). In terms of numerical implementation, this version of the model uses the finite-element method with triangular elements, and has previously been applied to simulate the evolution of Hoffellsjökull (Aðalgeirsdóttir and others, 2011). The average slope of the studied outlet glaciers is in the range 3–4°, and the aspect ratio (typical thickness to characteristic length for the ice body) is on the order of 10⁻³, which is acceptable for application of the SIA method (e.g. Le Meur and others, 2004).

Several values for the rate factor (A) in Glen’s flow law are applied to select one that best simulates the observed glacier geometry. In previous studies (Aðalgeirsdóttir and others, 2006; Guðmundsson and others, 2009a) the applied rate factor was 6.8 × 10⁻¹⁵ s⁻¹ kPa⁻³, and in the more recent study of Aðalgeirsdóttir and others (2011) the rate factor A = 4.6 × 10⁻¹⁵ s⁻¹ kPa⁻³ was determined from a series of model runs for Hoffellsjökull to best fit the observed velocity field. Based on the results of several ice-flow model studies, the recommended value for the rate factor for temperate ice is A = 2.4 × 10⁻¹⁵ s⁻¹ kPa⁻³ (Cuffey and Paterson, 2010). Model runs including explicit basal sliding were made, applying a Weertman-type sliding law (Paterson, 1994), where the sliding velocity (V _s) is assumed proportional to the basal shear stress, τ _b, to a power of m (V _s = Cτ _b ^m ). C is the sliding parameter, and the exponent m is assumed equal to 3 in our model runs. No data are available to determine the relative contributions of deformation and sliding velocities to the glacier flow. The chosen horizontal resolution of the input files and the triangular elements of the model mesh was 200 m. The ice divides are kept at a fixed location in the model and no flow is allowed across them. Thus no ice-divide migration due to mass-balance changes is permitted. The ice-flow model is run both coupled and uncoupled with the mass-balance model using LT-DP in the time-dependent simulations. When coupled, the changing glacier surface elevation is updated each year. This affects the calculated mass balance, whereas in uncoupled runs the mass balance is computed on the same topography throughout.

Mass-Balance Simulations with Various Precipitation Data

In this section we describe the model results derived using various precipitation data for the mass-balance model: (1) constant precipitation gradients, (2) an average winter mass-balance grid based on in situ measurements revised with precipitation data from the WRF model, and (3) LT-DP as input. Constant temperature and precipitation forcing is applied in all the simulations, using the averages of the baseline period 1980–2000. The model simulations all start with a smoothed measured 2000 glacier surface DEM (obtained by running the model for 2 years), and a rate factor of A = 4.6 × 10⁻¹⁵ s⁻¹ kPa⁻³.

Constant precipitation gradients and new stake data

Simulations with the mass-balance/ice-flow model were carried out using constant vertical and horizontal precipitation gradients: 0.5 per 100 m in the vertical, 0.0005 and −0.0008 per 100 m in the horizontal east and north directions, respectively. Temperature at Hólar in Hornafjörður and precipitation at Fagurhólsmýri were used to force the model, with the same degree-day factors for snow and ice as in previous studies on south Vatnajökull (Aðalgeirsdóttir and others, 2006). These simulations revealed a surface lowering of >90 m in the accumulation area, and a volume reduction of 20–30% between the start and end of the simulation (Fig. 6a; Table 5). The results indicate that the mass balance is underestimated in the accumulation area of the glaciers by the constant precipitation gradients.

Table 5. Difference (%) between volume and area at the start and end of 300 year simulations using the mass-balance model with constant precipitation gradient (CPG), the winter mass balance (b _w) revised with downscaled precipitation data from the WRF model and the mass-balance model using LT-DP, with rate factors between A = 4.6 × 10⁻¹⁵ s⁻¹ kPa⁻³ and A = 2.4 × 10⁻¹⁵ s⁻¹ kPa⁻³, and including basal sliding (C = 10 × 10⁻¹⁵ m a⁻¹ Pa⁻³). Constant temperature and precipitation forcing is applied using the average of the baseline period 1980–2000, and starting with the smoothed 2000 glacier surface DEM (Figs 6 and 8)

Fig. 6. The elevation difference between the initial measured surface from 2000 and the final simulated glacier surface of a 300 year model run using constant temperature and precipitation forcing of the baseline period 1980–2000. (a) The mass-balance model using a constant precipitation gradient (CPG). (b) Interpolated mean winter mass-balance grid for the period 1996–2006 revised with the downscaled winter precipitation from the WRF model. (c, d) Mass-balance model using downscaled precipitation from the LT model (LT-DP). In the simulations in (a–c) a rate factor of A = 4.6 × 10–15 s⁻¹ kPa⁻³ is used. In the flow model and the simulation shown in (d) a rate factor of A = 2.4 × 10–15 s⁻¹ kPa⁻³ is applied.

Experience from fieldwork led to the suspicion that winter accumulation on the southern slopes of Breiðabunga is higher than on the ice divide (as measured at BB0 and B19; Fig. 7a). Two new stakes (Sk01 and Fl01; Fig. 7b) were therefore added to the mass-balance network in the fall of 2009, which confirmed that snow accumulation is greater south of the plateau (Fig. 7b). It is clear that the stakes on the ice divides do not adequately represent the winter mass-balance variation of the studied outlets (Fig. 7a).

Fig. 7. (a, b) Manually interpolated average winter balance (b _w) maps for (a) 1996–2009, based on in situ measurements, and (b) 2010–13, when stakes Skf01 and Fl01 have been added to the mass-balance network. (c) Average winter precipitation (av.wP) map for 1996–2006 simulated with the WRF model. (d) Average winter precipitation map for 1996–2004 simulated with the LT model.

Winter mass-balance maps from in situ measurements revised with downscaled precipitation data from the WRF model

In order to explore the possible deficit of the winter mass balance in the accumulation area of the three glaciers, as found when using the constant precipitation gradients in the mass-balance model, two independent modelled precipitation datasets were analysed: the downscaled winter precipitation fields derived from the WRF model and from the LT model (Fig. 7c and d). Both datasets indicate that the winter precipitation distribution is more complex in this area than can be accounted for by constant precipitation gradients. The greatest snow accumulation is to the south of the Breiðabunga plateau according to both model results (Fig. 7c and d).

The manually interpolated winter mass-balance maps of the period 1996–2006 were revised with the downscaled precipitation data from the WRF model. Summer melt was deduced from the degree-day model, using the average temperature from Hólar in Hornafjörður in the period 1980–2000. Forcing the ice-flow model with the interpolated mass-balance maps improves the simulated glacier geometry. Less surface lowering occurs in the accumulation area, and the glaciers reach equilibrium without losing the same ice volume as when forcing it with the constant precipitation gradients (Fig. 6a and b; Table 5). The outlets advance beyond the 2000 margin (Fig. 6b), resulting in a 5–14% increase in area compared to the area in 2000 (Table 5).

LT-DP and tuning of the ice-flow rate factor

Figure 6c shows the elevation difference between the start and end of the simulation. The surface lowering in the accumulation area is similar to when using the revised mass-balance maps with the WRF data, but less thickening in the ablation area (cf. Fig. 6b and c). The volume loss is on the order of 2–12%, and the simulated areal extent is slightly larger than the observed area (Fig. 6c; Table 5). The glaciers undergo less surface lowering in the accumulation area than in earlier simulations (Fig. 6a–c), and the modelled volume and area are closer to the observations, the difference being within 3% (Table 5). In view of the good agreement between observed and simulated glacier geometry, a series of simulations with various rate factors and adding explicit sliding were carried out (Table 5). The best fit (both volume and surface profiles) is obtained with a rate factor of A=2.4 × 10⁻¹⁵ s⁻¹ kPa⁻³, representing stiffer ice (Fig. 6d). Adding sliding (C = 10 × 10⁻¹⁵ m a⁻¹ Pa⁻³) gives the same results as having softer ice (A = 4.6 × 10⁻¹⁵ s⁻¹ kPa⁻³) (Table 5).

Summary of simulations using various precipitation data for the mass-balance model

A further comparison of the simulations described above is presented in Figure 8, where longitudinal surface profiles from the ice divide to the terminus are shown, and compared with a number of observed surface profiles. The substantial lowering of the glacier surface using constant precipitation gradients in the mass-balance model is evident on the modelled surface profiles on all three outlets. The simulations using the revised winter mass-balance map with precipitation from the WRF model indicate an overestimation of the winter mass balance, as all three outlets advance ∼1 km beyond the observed margin, and thickening in the ablation zone of up to 100 m (Fig. 6b). Simulations using the LT-DP show the best fit with the observed surface profiles (Fig. 8).

Fig. 8. Longitudinal profiles (locations shown in Fig. 2b) showing the glacier thickness and extent at different times according to observations (filled-in mauves, blues and green areas) from Hannesdóttir and others (2014a) and the various simulations (lines) using constant temperature and precipitation forcing (average of the baseline period 1980–2000), except for the dark-blue surface profile (−1°C and −20% precipitation). Two different rate factors are applied (A = 4.6 × 10–15 s⁻¹ kPa⁻³ and A = 2.4 × 10–15 s⁻¹ kPa⁻³). Three different precipitation distributions are applied: CPG = constant precipitation gradient; winter mass balance (b _w) map revised with downscaled precipitation from the WRF model; and LT-DP. The white-shaded horizontal bands indicate the elevation range of the snowline, derived from the MODIS images.

Figure 9a and b show the modelled winter mass balance and the manually interpolated winter mass-balance map, based on the in situ measurements (including the new stakes). There is good agreement between the two maps, both showing the highest winter mass-balance values south of the ice divide. The modelled and measured specific winter balance is given for each outlet in the figure, showing similar values. The mass-balance model was not able to reproduce the abnormally low summer balance of 2010 (Fig. 3b), which was affected by the tephra from the Eyjafjallajökull eruption deposited on the ice cap (Guðmundsson and others, 2012) and caused increased absorption of shortwave radiation and glacier melt (e.g. Björnsson and others, 2013). For further validation of the mass-balance model using the LT-DP, the end-of-summer snowline (which is used as a proxy for the ELA, derived from a set of MODIS images from 2007 to 2010, from Hannesdóttir and others, 2014a) is plotted on the modelled net mass-balance maps in Figure 9c–f. The location of the MODIS snowline agrees well with the modelled ELA in all four years.

Fig. 9. (a, b) Winter mass balance (b _w) in 2010: (a) modelled using LT-DP and (b) manually interpolated b _w (based on in situ measurements) in m w.e. a⁻¹. The specific winter balance is shown for each glacier: Skálafellsjökull (S), Heinabergsjökull (H) and Fláajökull (F). (c–f) Modelled net balance (b _n) in m w.e. a⁻¹ in 2007 (c), 2008 (d), 2009 (e) and 2010 (f) and the ELA (b _n = 0) drawn in pink. The MODIS end-of-summer snowline (a proxy for the ELA) is shown with black dots (from Hannesdóttir and others, 2014a).

Simulations with LT-DP

Step changes in temperature and precipitation

Steady-state experiments were performed to analyse the sensitivity of the model to step changes in temperature and precipitation relative to the baseline period 1980–2000 (Table 6). In reality, glaciers are constantly responding to climate variations and are never in a steady state, but such simulations provide insight into the sensitivity of glaciers to changes in climate. All simulations start with a geometry that has been run into a steady state. Each model simulation is forced with constant mass-balance grids produced by the mass-balance model. The applied step changes in temperature are in the range of the observed temperature difference between the late 19th century and 2010 observed at Hólar in Hornafjörður (Table 3), and possible future climate scenarios. The modelled glaciers respond to step changes in precipitation and temperature by increasing or decreasing their volume and extent (Table 6). The mass-balance sensitivity to a 1°C temperature decrease is equivalent to a 40% increase in annual precipitation (Table 3), and similarly the sensitivity to a 1°C increase in temperature is close to a 40% decrease in annual precipitation (Table 6).

Table 6. Mass-balance sensitivity of the three outlet glaciers to given deviations from the climate of the baseline period 1980–2000, and the per cent difference between the volume, and area, at the start and end of steady-state simulations, using step changes in temperature and precipitation, relative to the average of the baseline period. The bottom row shows the results of using the average precipitation and temperature of the period 2000–10

Simulations were made to assess the sensitivity of the glaciers to possible future temperature and precipitation changes, according to predictions for the 21st century in Iceland (e.g. Jóhannesson and others, 2007; Nawri and Björnsson, 2010). Temperature is predicted to increase at a rate close to 0.3°C decade⁻¹ until 2050, and 0.2°C decade⁻¹ until 2100, with superimposed decadal variations determined by natural climate variability (Nawri and Björnsson, 2010). Changes in precipitation are projected to be moderate, and the increase in annual precipitation from 1961–90 to 2071–2100 is likely to be in the range 0–10% in most regions (Jóhannesson and others, 2007; Nawri and Björnsson, 2010). A 2°C warming and a 10% increase in annual precipitation would eventually lead to >50% volume loss, 8–10 km terminus retreat and up to 35% area loss (Fig. 10b; Table 6). A temperature increase of 3°C and a 10% increase in annual precipitation could lead to 80–90% volume loss (Table 6; Fig. 10c). The glaciers undergo similar volume loss if the annual precipitation is kept the same as in 1980–2000 (Table 6), as a higher percentage of precipitation will fall as rain at higher temperatures.

Fig. 10. (a) The simulated LIA glacier geometry obtained with a constant 1°C cooling and 20% decrease in annual precipitation, relative to the baseline period 1980–2000, and the reconstructed LIA extent in dark blue (Hannesdóttir and others, 2014b). (b, c) The glacier geometry at the end of simulations forced with 10% step increase in annual precipitation and step increase in temperature (b) 2°C and (c) 3°C, relative to the baseline period. The red glacier margin is the result of simulations without any change in annual precipitation, only increase in the temperature, and the lidar margin of 2010 is shown in blue.

Simulating the ∼1890 volume

Model runs were carried out to simulate the reconstructed ∼1890 glacier geometry and volume of the three outlets (as described in Hannesdóttir and others, 2014b). Temperature measurements, from the late 19th century from Hólar in Hornafjörður (Table 3) indicate 1°C lower temperatures relative to the baseline period. Three different precipitation variants (90%, 85% and 80% of the average annual precipitation of 1980–2000) accompanied by the 1°C lower temperature were applied (Fig. 10). The simulation that best resembles the reconstructed ∼1890 glacier geometry was obtained using a 20% reduction of annual precipitation, which is in accordance with the precipitation estimates over the period 1860–90 (from Aðalgeirsdóttir and others, 2011; Table 3).

The simulated ∼1890 volume of the three outlet glaciers is within the error estimates of the reconstructed (observed) volume derived from geomorphological field evidence (Table 2). The simulated ∼1890 areal extent is slightly smaller than the observed area, except for Fláajökull. The disagreement between the observed and simulated areas can be explained by the small tongues on the west side of Fláajökull, reaching beyond the mapped 1890 glacier margin in the simulations (Fig. 10a). The ELA of the simulated LIA glaciers is around 800 m on Skálafellsjökull, 750 m on Heinabergsjökull and 800 m on Fláajökull, in good agreement with the estimated LIA ELA, which is based on the elevation of the uppermost lateral moraines (Hannesdóttir and others, 2014b). This is a known method, where the lateral moraines are used as a proxy for the ELA, since lateral moraines are only deposited below the ELA, where glacier ice is emerging (e.g. Benn and Evans, 2010).

Simulating the 1959–2010 evolution

The time-dependent simulations for the period 1959–2010 were initialized with a steady-state glacier geometry, and run both with coupled and uncoupled mass-balance ice-flow models (Fig. 11). The model simulates the observed volume changes well (Fig. 11). The simulations show a slight mass gain in the 1990s, which has also been observed as slight advances of the termini and gain in the glacier surface elevation in the accumulation area (Hannesdóttir and others, 2014a). The observed volume decrease from 2002 to 2010 is not entirely captured by the model, and this is especially noticeable for Heinabergsjökull. The difference between the observed and simulated volumes in 2010 is 9% of the observed 2010 volume. This is larger than the 1.5% and 3.2% differences for Skálafellsjökull and Fláajökull, respectively.

Fig. 11. The volume evolution of the time-dependent simulations, from coupled (dotted lines) and uncoupled (solid lines) runs using the temperature from Hólar in Hornafjörður and LT-DP during the period 1959–2010. The observed volumes of 1890, 1904, 1945, 1989, 2002 and 2010 are shown with open circles (from Hannesdóttir and others, 2014a). The simulated volume of 1890 is shown with x, assuming that temperatures were 1°C lower than during the baseline period, 1980–2000, and the annual precipitation 20% less (see Fig. 10a for glacier extent).

There is a 1–2% difference in volume between the coupled and uncoupled simulations in the model runs; all glaciers experience greater volume loss in the coupled runs when the mass balance is computed on a lowering surface (Fig. 11; Table 2). This emphasizes the importance of using coupled models for future simulations, when the glaciers will undergo increased surface lowering and retreat, as a result of a warmer climate.

Discussion

The importance of using downscaled precipitation in the simulations

Located on the southeast coast of Iceland, Skálafellsjökull, Heinabergsjökull and Fláajökull receive high amounts of precipitation in southerly and easterly winds associated with the frequent passage of atmospheric fronts and lows along the coast. From the model simulations it is evident that the precipitation distribution is not adequately described by constant horizontal and vertical precipitation gradients. Orographically enhanced precipitation and snowdrift complicates the snow distribution in the accumulation area of the three outlet glaciers. The two new mass-balance stakes (Skf01 and Fl01) confirm the modelled winter mass-balance pattern of the LT-DP and WRF simulations. Good agreement between the modelled and measured winter balance during the only overlapping year (2010) on the three outlets (Fig. 9a and b) and between the MODIS-derived and modelled ELA (Fig. 9c–f) indicates that the 1 km model resolution is adequate to simulate the precipitation distribution in the complex topography. The 3 km resolution of the downscaled precipitation from the RÁV data is not sufficient for the outlets of southeast Vatnajökull. However, it is adequate to simulate the atmospheric flow and the spatial structure of the precipitation field on the relatively flat Mýrdalsjökull ice cap (Ágústsson and others, 2013). This may be attributed to the more complex topography along the southeast coast than in the accumulation area of Mýrdalsjökull.

Constant precipitation gradients may be sufficient when simulating one outlet glacier (as has successfully been done for Hoffellsjökull (Aðalgeirsdóttir and others, 2011)) with more mass-balance data to constrain the model than the three outlets. Modelling a group of outlets separated by steep mountains that influence the orographically enhanced precipitation, however, is more complex. Also, a wealth of data are available for Hoffellsjökull, but not the outlet glaciers for this study. In addition to the mass-balance measurements, data include glacier velocity fields from GPS surveys and satellite data, and energy-balance measurements from an automatic weather station, which were used in the calibration of the mass-balance model (Aðalgeirsdóttir and others, 2011).

The average temperature at Hólar in Hornafjörður in the period 2000–10 was 10.5°C, compared to 9.6°C during the baseline period 1980–2000 (Table 3). Simulations using the average temperature and precipitation of the period 2000–10 result in volume loss of 7–12% (of the 2000 volume) (Table 6), which is less than the volume loss due to a 1°C warming relative to the baseline period (23–34%). Due to this difference, it can be concluded that considerably more precipitation fell in the mountains and on the outlet glaciers of southeast Vatnajökull than was observed at the lowland meteorological stations during the same period. Recent experiments with numerical weather models in this area indicate that changes in precipitation in mountainous terrain are not necessarily correlated with changes in precipitation in lowlands (Ágústsson, 2014).

Sensitivity of the outlet glaciers to climate perturbations

The importance of precipitation for the mass balance of these outlet glaciers is highlighted when comparing the results of simulations forced with and without a change in precipitation. A 1°C cooling and unchanged precipitation, compared to a 1°C cooling and 20% less annual precipitation, relative to the baseline period 1980–2000, results in a 30–50% increase in ice volume and a 10–18% increase, respectively (Table 6). To counterbalance a volume increase due to a 1°C cooling, the precipitation would need to decrease by ∼40% (Table 6), which is likely outside the range of natural variability. This is similar to previous results of mass-balance modelling studies where increases in precipitation of 20–50% and 30–40% were required to balance increased ablation due to a 1°C temperature rise in a variety of climate regimes (Oerlemans and others, 1998; Braithwaite and Raper, 2002). The three outlets respond differently to warming and cooling scenarios (Table 6). Considerable deviations are observed between the modelled and observed surface profiles for Heinabergsjökull in most of the simulations (Fig. 8), whereas there is better agreement between modelled and observed surface profiles for the other two outlets. Heinabergsjökull is more sensitive to changes in temperature and precipitation. For example, the increase or decrease in volume of Heinabergsjökull in response to changes in precipitation and/or temperature is nearly twice that of Skálafellsjökull (Table 6). The hypsometry (area distribution with elevation) and basal topography of these outlets differ. Heinabergsjökull terminates in an overdeepened basin, with a large part of the bed below sea level (Fig. 2a), and its terminus is at low elevation. The proglacial lake of Heinabergsjökull probably affects the mass balance and enhances the ablation by melting and calving, as has been observed and calculated for Hoffellsjökull (Aðalgeirsdóttir and others, 2011). This additional mass loss is not taken into account in the flow model, as no calving mechanism is included. This can explain the deviation between the modelled and measured ice volumes in the first decade of the 21st century (Fig. 11). The lake area has not increased during the 8 year period. However, the surface topography of the terminus has changed during this period, indicating that the glacier tongue is now floating. A possible change in mean glacier density may also explain why the model does not capture the 2002–10 volume loss, as a constant density is assumed in the simulations.

Conclusion

The known geometry and volume evolution of three outlet glaciers of southeast Vatnajökull has been simulated using a coupled numerical mass-balance/ice-flow model. Forcing the mass-balance model with constant horizontal and vertical precipitation gradients does not successfully describe the spatial variance in the winter mass balance in this area. The results indicate that the 1 km resolution of the LT model is adequate to simulate the precipitation fields on the glaciers, and that downscaled precipitation data are necessary for mass-balance modelling of glaciers in complex mountain terrain. The outlet volume around 1890 is successfully simulated by forcing the model with temperatures 1°C lower than average for the reference period 1980–2000, which is in line with data from nearby meteorological stations outside the glaciers and a decrease in annual precipitation of 20%. The calculated sensitivity to 1°C cooling is similar to a 40% increase in annual precipitation, or in the range 1.0–1.5mw.e. a⁻¹. Warming of 2–3°C (relative to the baseline period 1980–2000) by 2100, which some forecasts predict, will result in a >50–80% decrease in ice volume of the southeast outlet glaciers. For time-dependent future simulation, downscaled precipitation fields are necessary.