Modelling the 20th and 21st century evolution of Hoffellsjökull glacier, SE-Vatnajökull, Iceland

. The Little Ice Age maximum extent of glaciers in Iceland was reached about 1890 AD and most glaciers in the country have retreated during the 20th century. A model for the surface mass balance and the ﬂow of glaciers is used to re-construct the 20th century retreat history of Hoffellsjökull, a south-ﬂowing outlet glacier of the ice cap Vatnajökull, which is located close to the southeastern coast of Iceland. The bedrock topography was surveyed with radio-echo soundings in 2001. A wealth of data are available to force and constrain the model, e.g. surface elevation maps from ∼ 1890, 1936, 1946, 1989, 2001, 2008 and 2010, mass balance observations conducted in 1936–1938 and after 2001, energy balance measurements after 2001, and glacier surface velocity derived by kinematic and differential GPS surveys and


Introduction
Iceland (103 000 km 2 ) lies in the North Atlantic Ocean, just south of the Arctic Circle.Due to the warm Irminger Current, the island enjoys a relatively mild and wet oceanic climate and a small seasonal variation in temperature.The average winter temperatures are around 0 • C near the southern coast, where the average temperature of the warmest month is only 11 • C and the mean annual temperature is about 5 • C (Einarsson, 1984).At present about 11 % of the country is covered by glaciers (Björnsson and Pálsson, 2008).The Icelandic ice caps are temperate, characterized by high annual mass turnover (1.5-3.0 m water equivalent (w.e.)) and are highly dynamic.They are sensitive to climate variations and have responded rapidly to changes in temperature and precipitation during historical times (Björnsson, 1979;Björnsson et al., 2003;Björnsson and Pálsson, 2008;Guðmundsson et al., 2011).The recorded volume and area changes of these glaciers are therefore good indicators of climate change.
Historical records of the glacier extent reach back to the settlement of Iceland in the late 9th century.During the settlement, glaciers were smaller than at present.They started to advance in the 13th century at the onset of the Little Ice Age that lasted until late 19th century when most glaciers in Iceland reached their maximum extent.In the 20th century, the climate was significantly warmer than during the Little Ice Age, with higher temperatures in the period 1930-1945 and again at the end of the century.Icelandic glaciers have retreated during most of the 20th century with the exception of a standstill or an advance in   (Sigurðsson and Jóhannesson, 1995;Björnsson and Pálsson, 2008).
In this study, a two-dimensional Shallow Ice Approximation (SIA) ice-flow model coupled with a positive degreeday (PDD) mass-balance model (Aðalgeirsdóttir et al., 2004(Aðalgeirsdóttir et al., , 2006) ) is used to simulate the evolution of Hoffellsjökull outlet glacier of the Vatnajökull ice cap.The model is run from the Little Ice Age maximum extent, at the end of the 19th century, throughout the 20th century and then used to predict the future response of the glacier to an ensemble of climate change scenarios.Temperature measurements from the meteorological station Hólar in Hornafjörður (HH) and precipitation measurements from the station Fagurhólsmýri (F), (15 and 85 km away from Hoffellsjökull, respectively; locations shown in Fig. 1a), are used to force the coupled model.A sensitivity study of various model parameters and model assumptions is carried out and an ensemble of climate change scenarios (Jóhannesson et al., 2011) is used to assess the relative importance of natural climate variability and a deterministic anthropogenic warming trend for the evolution of the glacier during the 21st century.

Study area
Hoffellsjökull is a south-flowing outlet glacier of Vatnajökull, the largest ice cap in Iceland (Fig. 1).The accumulation zone is a part of the eastern sector of the ice cap between two mountain ranges, the subglacial Breiðabunga and Goðahnúkar, at 1350-1450 m elevation.From the accumulation area, the ice flows in two branches, west and east of the central nunatak Nýju Núpar.The branches meet below the nunatak and the ice is funnelled through a 2 km wide ice fall at 600-700 m elevation where the elevation drops by 300 m over a 4 km distance.Below the ice fall, the glacier spreads out forming two branches, west and east of a prominent ridge.Currently the eastern branch terminates in a lagoon at ∼40 m a.s.l. and the western branch in a lagoon at ∼20 m a.s.l.(local names are Hoffellsjökull and Svínafellsjökull for the eastern and western branch, respectively).
The first glaciological expeditions to this area were conducted in [1936][1937][1938], when a group of Swedish and Icelandic glaciologists measured ice flow, surface mass balance and surface topography.They also carried out a detailed analysis to understand the relative roles of accumulation and melting in the total mass balance of the glacier and to establish a relationship between the climate and the advance and retreat of the glaciers.An up to 8 m thick winter snow layer was measured in the accumulation area (∼4 m w.e.).Ice melt of up 10 m w.e. was measured in the lowest part of the ablation zone in summer, and 2 m w.e. was melted during winter.Taking into account ∼2 m of annual rainfall, the runoff from this part of the glacier was estimated as ∼14 m w.e. per year; a surprisingly high value (Ahlmann, 1939;Ahlmann and Thorarinsson, 1943;Thorarinsson, 1939Thorarinsson, , 1943)).
The glacier was revisited in 2001 when the surface and bed topographies were mapped.In 2002, automatic weather stations were placed at 2 locations on the glacier and a still ongoing surface mass balance survey was initiated.In addi- tion, the glacier mass balance has been measured since 1996 at four locations close to the ice divide of Hoffellsjökull to the west and north.The mass balance of the nearby southern outlet, Breiðamerkurjökull, has been measured since 1996 (Fig. 1b).

Geometry
Hoffellsjökull was surveyed with a radio-echo sounder and Differential Global Positioning System (DGPS) equipment in 2001.Continuous profiles, approximately 1 km apart, were measured in the accumulation zone and a dense network of point measurements were carried out in the ablation zone.Digital Elevation Models (DEMs) of the surface and bedrock were created from these data (Fig. 2; Björnsson and Pálsson, 2004).The estimated errors are at most 1-5 m (bias less than 1 m) for the surface map and depending on the location, 5-20 m for the bedrock map.Glacier extent and maps of the ice surface of Hoffellsjökull are available from the years 1904 (lowest parts, geodetic survey in summer), 1936 (geodetic survey in summer), 1946 (aerial photographs in autumn), 1989 (aerial photographs in autumn), 2001 (DEM from DGPS survey in spring), 2008 (DEM from SPOT5 HRS images in autumn, Korona et al., 2009) and 2010 (airborne LiDAR in autumn).
The Cryosphere, 5, 961-975, 2011 www.the-cryosphere.net/5/961/2011/The most accurate glacier DEM is created from the airborne LiDAR survey conducted in the autumn of 2010 (5 × 5 m pixel resolution, with an accuracy of <20 cm in elevation and <0.5 m in position, Fig. 3).It is used as a reference map for co-registering the SPOT5 HRS-DEM (pixel resolution of 40 × 40 m), using the ice free areas surrounding Hoffellsjökull and the correlation method described by Guðmundsson et al. (2011).This comparison revealed a horizontal shift of the HRS-DEM by 15 m and 5 m towards east and north, respectively, and a vertical offset of 2.3 m.Gaps in the HRS DEM, due to low contrast of the SPOT5 stereo image pairs at some sections in the accumulation area of Hoffellsjökull, were filled by smoothly adjusting the LiDAR DEM to the HRS DEM.Available elevation values within the accumulation area were used to calculate the offset and tilting between the 2008 and 2010 DEMs.After correcting the HRS DEM with the LiDAR DEM, we estimate the residual vertical bias error of the HRS DEM to be <0.5 m; as obtained in a similar study by Guðmundsson et al. (2011).
The DEMs from before 2001 were created from digitized contour lines of available older maps (Fig. 2).The ice surface elevation of 1890 (Fig. 2b, Little Ice Age maximum, LIA max ) was reconstructed from the more recent surface DEMs assuming that although the surface has lowered significantly, the spacing of the elevation contours in the accumulation zone only changes slightly.The location of the LIA max terminal moraines and the location and elevation of the side moraines on both sides of the glacier (Björnsson and Pálsson, 2004;Jónsson, 2004) as measured by DGPS were used in the reconstruction of the 1890 ice surface elevation of the lower part of the glacier as well as the ice surface map from 1904 which was used as a constraint for the curvature of the elevation contours.A conservative vertical error estimate for the reconstructed 1890 DEM is 15-20 m, 10-15 m for the 1936 DEM, 5-10 m for the 1946 DEM and 5 m for the 1989 DEM.
Additional information on glacier extent and glacier variations may be extracted from written historical descriptions of local farmers from the 18th and 19th century.They describe the land use (i.e.grassing of cattle and sheep, use of the forest growing in the valley etc.), the advance of the glacier during the LIA, and the frequency and size of jökulhlaups from glacier dammed lakes.From these descriptions it is known that the glacier terminus advanced at least 5-7 km during the LIA (Jónsson, 2004).
www.the-cryosphere.net/5/961/2011/The Cryosphere, 5, 961-975, 2011 Much of the current ablation area of Hoffellsjökull lies in an over deepened trench.Presently, the trench reaches down to ∼300 m below sea level under the eastern branch of the glacier (Fig. 2).This trench is the upper part of a valley carved into the bedrock by glaciers during the past glaciations and that extends towards the Iceland shelf margin.Before the LIA advance, the trench was filled with loose sediments, similar to the sediment plains in large areas of SE-Iceland.The historical documents describe that the glacier advanced during the 18th and 19th centuries over a vegetated plain (at ∼40-80 m a.s.l.), where now there is an 11 km 2 trench excavated by the glacier.The radio-echo sounding measurements show that the volume of the trench is about 1.6 km 3 (Björnsson and Pálsson, 2004).This volume consisted of loose, fine grained sediments that are easily transported by fluvial processes at the glacier sole.The advance of the glacier over sediment plains started ∼400 yr ago.Sediment removal rates similar to that described by Björnsson (1996) at the nearby outlet Breiðamerkurjökull would not be sufficient to remove all the sediments.Jökulhlaups, however, are highly effective in removing and transporting sediments.Jökulhlaups from the many marginal lakes of Hoffellsjökull started ∼1840 (Björnsson, 1976;Hjulström, 1954;Arnborg, 1955a, b), with increasing frequency and larger volumes in the period 1920-1960 as the glacier lowered and retreated.We suggest that most of the loose sediments were already removed by the mid 20th century.Longitudinal and cross sections in Fig. 4 show the thickness and shape of the glacier at the different times.The maximum thickness of the glacier of ∼560 m is reached about 5 km upstream from the terminus.
A widespread break-up of the over-deepened eastern branch of the terminus area happened in 2010 as is shown in the high resolution LiDAR DEM in Fig. 3.This is likely to be the beginning of a formation of a terminus lake with a calving glacier front, that will gradually grow as the glacier retreats.

Mass-and energy balance
The surface mass balance was measured at 10 locations (fully covering the glacier) during the glaciological years 1935/1936-1937/1938 by the Swedish-Icelandic expedition (Fig. 5; Ahlmann and Thorarinsson, 1943), at two locations on the glacier (one close to the terminus and the other in the central part of the accumulation zone) since the glaciological year 2001/2002, at two locations close to the western and northern margins of the accumulation zone, and at a number of sites on nearby outlet glaciers of SE-Vatnajökull since 1996 (location of sites shown in Fig. 1b, c).In spring (April-May) cores are drilled through the winter snow layer and the density is measured.Stakes or wires are left in the core holes (in the ablation zone, ∼10 m long wires are drilled into the ice with a steam drill).The summer balance is measured from the extension of the stakes or wires in the spring and late autumn (September-October; Björnsson et al., 1998Björnsson et al., , 2003)).
The mass balance and climate conditions during the 1936-1938 Swedish-Icelandic expedition were similar to the first decade of the 21st century (Figs. 5 and 8).Furthermore, the relationship between mass balance and elevation at Hoffellsjökull during the period 1936-1938 is similar as observed at the nearby outlet glacier Breiðamerkurjökull (see Fig. 1 for location) since 2001 (Fig. 5).Digital mass balance maps of Hoffelssjökull for each year after 2001-2002 have been generated using the observed mass balance and the observed elevation dependency of mass balance at stakes on SE-Vatnajökull.The resulting mass balance maps are integrated over the whole glacier area to calculate the mean specific mass balance given in Table 1.
The mass balance measurements have been supplemented by glacio-meteorological observations at both mass balance sites on the glacier; station Hof at ∼1140 m a.s.l., slightly above the current ELA, and station HoSp at ∼100 m a.s.l.within the ablation zone (Fig. 1).The weather parameters observed on the glacier have been used to calculate the full energy balance at the two AWSs with the methodology described by Guðmundsson et al. (2009a).Both the mass balance and the energy balance data have been used to calibrate and validate the PDD mass balance model applied here (see Sect. 5).

Surface velocity
The ice surface velocity was measured during the 1936-1938 Swedish-Icelandic expedition at several stakes along the transverse profile DD' shown in Fig. 4 (Thorarinsson, 1943).There is a large variability in the rate of movement from one year to another (Fig. 6a).The horizontal ice velocity has been observed since the glaciological year 2001/2002.The positon of mass balance stakes is measured with kinematic or differential GPS land survey instruments in all visits to the sites.The summer and sometimes annual or winter velocities are calculated from measured positions in spring, midsummer and autumn (Fig. 6b).Late summer and a sparse annual velocity maps have been deduced from correlation of 2.5 m resolution SPOT5 HRG images (Fig. 7 and Table 2); with an accuracy of about half the pixel size (Berthier et al., 2005).The magnitude of the SPOT5 velocity in profile DD' is similar to the observations made in the 1930s (Fig. 6a).A good spatial coverage was obtained by using two SPOT images from late summer 2002, but reliable signals were only obtained for limited areas when an attempt was made to determine the annual velocity fields from the two late summer SPOT5 images from 2002 and 2003.

Temperature and precipitation records from meteorological stations
In the present study, the mass balance of Hoffellsjökull is calculated by a degree-day model that uses daily mean temperature at Hólar in Hornafjörður and precipitation at Fagurhólsmýri as input parameters (the location of the stations is shown in Fig. 1).For this study, monthly mean data were available or estimated for the study period.To obtain daily records, a representative one-year data set of daily mean temperature and daily accumulated precipitation was constructed with random sampling of each day of the year from the period 1981-2000 which has daily data series The Cryosphere, 5, 961-975, 2011 www.the-cryosphere.net/5/961/2011/available.This was done to maintain the statistical relationship between the daily temperature and precipitation, which is important for realistic mass balance modeling.Daily temperature and precipitation time series were then constructed by adding the difference between the observed monthly mean value and the corresponding monthly mean for the period 1981-2000 to the data set of representative daily values.
The monthly records at the stations Hólar in Hornafjörður and Fagurhólsmýri are, however, not available for the whole study period and therefore it is necessary to use records from other stations to fill in the data gaps and extend the available series.Historical temperature and precipitation records from a number of stations around Iceland are available (Fig. 1 ).An iterative expectation-maximization (EM) algorithm (Dempster et al., 1977) was used to construct a continuous temperature record for Hólar in Hornafjörður extending back to 1860 using the temperature records from these stations.Precipitation measurements at Fagurhólsmýri were started in 1924.To create a record that goes back to 1860 to match the temperature record, we used a linear regression between the monthly values (separate regression for each month) of the temperature at Hólar in Hornafjörður and precipitation at Fagurhólsmýri, tuned with available precipitation observations (Fig. 8).The correlation between precipitation and temperature in the period 1931-2009 is 0.52 when all data are used and 0.5 when every second year is used for reconstruction and the others for testing.During summer months (May-August) there is almost no correlation with temperature.The winter balance is however not expected to be correlated with the precipitation of the summer months.Testing the method with data from the winter months only gives a higher correlation of 0.58 (using different data for optimization and testing).This regression method reconstructs the long term variability in the precipitation, but the resulting annual variability is smaller than the observed.

Future climate scenarios
In total 13 future climate scenarios from the Nordic Project Climate and Energy Systems (CES) were used to force the ice sheet model into the future.These scenarios are based on three dynamical downscalings of global Atmosphere-Ocean General Circulation Model (AOGCM) climate change simulations using the A1B emission scenario (Nakicenovic et al., 2000) performed with Regional Climate Models (RCMs) (ECHAM5-r3/DMI-HIRHAM5, HadCM3/MetNo-HIRHAM and ECHAM5-r3/SMHI-RCAO) and a data set of 10 global AOGCM simulations, also based on the A1B emission scenario, submitted by various institutions to the IPCC for its fourth assessment report (IPCC, 2007).These 10 GCMs were chosen from a larger IPCC data set of 22 GCMs based on their surface air temperature (SAT) performance compared with the ERA-40 reanalysis in the period 1958-1998 in an area in the N-Atlantic encompassing Iceland and the surrounding ocean (Nawri and Björnsson, 2010).The temperature and precipitation simulated by the models was averaged over Iceland to provide single time-series intended to represent the climate development for the whole country.Based on the downscaled RCM model output the temperature change of the GCM-based scenarios was increased over all seasons and the entire period by 25 % in the interior of Iceland, where the large ice caps are located (Nawri and Björnsson, 2010).
Before year 2010, the glacier model is forced with daily mean records constructed from the monthly mean observed temperature and precipitation as previously explained.Possible natural variations in the climate are important for nearfuture projections as the magnitude of the expected anthropogenic change has then not exceeded the random variability of the climate.Therefore, many different climate scenarios were derived for the CES project (Jóhannesson et al., 2011).Expected values of temperature and precipitation in 2010 were estimated by statistical autoregressive (AR) modelling of the past records, thereby taking into account the warming that has been observed in recent years as well as the inertia of the climate system so that the very high temperatures of the last few years have only a moderate effect on the derived expected values.These expected values are intended to be representative for the deterministic part of the recent variation in the climate when short-term climate variations have been removed by the statistical analysis.Scenarios of the monthly mean temperature and precipitation were calculated from year 2010 to the end of the 21st century by fitting a least squares line to monthly values simulated by the RCM www.the-cryosphere.net/5/961/2011/The Cryosphere, 5, 961-975, 2011 or GCM from 2010 onwards and shifting the simulated timeseries so that the 2010 value of the least squares line matched the expected 2010 value based on the AR modelling of the past climate.
The trend analysis of the future climate eliminates the direct use of a past baseline period in the derivation of the scenarios and provides a consistent match with the recent climate development.Furthermore, the statistical matching of the past climate observations with the trend lines of the future climate provides an implicit bias correction (Jóhannesson et al., 2011).Figure 9 shows the simulated annual mean temperature and precipitation averaged over Iceland for the 13 scenarios used in this study, compared with the average 2000-2009 climate.All the scenarios indicate a slight increase in the precipitation during this century (∼10 %) and a warming of 1.1-1.5 • C near the middle of the 21st century and ∼2-3 • C by the end of the century, relative to the 2000-2009 average.When compared to the reference period 1981-2000, when the net balance of most of the Icelandic ice caps was close to zero (Guðmundsson et al., 2011;Aðalgeirsdóttir et al., 2006;Guðmundsson et al., 2009b), the temperature increase is 2.0-2.4 • C near the middle of the 21st century and ∼3-4 • C at the end of the century.

Mass balance model
The mass-balance of Hoffellsjökull is simulated with a positive degree-day model (PDD), that has been applied on several Icelandic glaciers, both independently (Jóhannesson, 1997;Jóhannesson et al., 1995Jóhannesson et al., , 2006) ) and coupled with an ice flow model for Langjökull, Hofsjökull and S-Vatnajökull ice caps (Aðalgeirsdóttir et al., 2006;Guðmundsson et al., 2009b; see locations of the ice caps in Fig. 1).Constant linear vertical and horizontal precipitation gradients (0.5 per 100 m in vertical, 0.005 and −0.008 per 1000 m in horizontal east and north direction, respectively) are used in the model to take into account observed precipitation variations, assuming a constant snow/rain threshold of 1 • C. In a previous study for Vatnajökull ice cap (Aðalgeirsdóttir et al., 2006) the model parameters; temperature lapse rate and the separate degree-day scaling factors for snow (ddf s ) and ice (ddf i ) were calibrated until the modeled mass balance matches the available mass balance observations on the whole S-Vatnajökull (up to 23 stakes, Fig. 1) from the mass balance years 1991/1992 to 2004/2005, using temperature at Hólar in Hornafjörður and precipitation at Fagurhólsmýri as an input.With a temperature gradient of 0.56 • C per 100 m and the degree-day factors ddf s = 4.45 mm • w.e.C −1 d −1 and ddf i = 5.30 mm • w.e.C −1 d −1 , the model explains 92 % and 95 % of the variance of the winter and summer balance at S-Vatnajökull, respectively (Jóhannesson et al., 2007).
Another model calibration of the degree-day factors for Hoffellsjökull was done in this study, using the two mass balance measurements on the outlet glacier itself and the stake at 1260 m a.s.l.near the ice divide of Hoffellsjökull that have observations since 2001.The modelled melt was also compared to the water equivalent of the daily melting energy calculated at the two AWSs on the glacier (locations in Fig. 1c) for validation of the degree-day model.By assuming the same temperature gradient as in the first calibration, this second calibration determines the degree day factors to be 4.0 ± 0.5 mm • w.e.C −1 d −1 and 5.3 ± 0.7 mm • w.e.C −1 d −1 for snow and ice, respectively (the given errors are one standard deviation (σ ) of the degree-day factors optimized separately for each year).These two model calibrations indicate that the degree day factors are robustly determined and the observed surface mass balance is well reproduced.

Ice flow model
The dynamic ice flow model is similar to that used by Aðalgeirsdóttir et al. (2006) for Hofsjökull and S-Vatnajökull and Guðmundsson et al. (2009b) for Hofsjökull and Langjökull, except that the numerical implementation in the model applied here uses a staggered finite element method on a triangular grid rather than a finite difference method to solve the continuity equation (Sigurðsson, 1992), allowing for variable grid size in the model domain.This model is based on the vertically integrated continuity equation and the shallow ice approximation (SIA), neglecting longitudinal stress gradients and surge dynamics, and excludes bed isostatic adjustments and seasonal variations in sliding (e.g.Aðalgeirsdóttir, 2003;Aðalgeirsdóttir et al., 2006).The geometry of Hoffellsjökull with a mean thickness of approximately 300 m and length of 20 km (Fig. 4) justifies the application of a SIA ice flow model.Although longitudinal stress gradients ignored by SIA models affect the flow of glaciers in areas of complicated or steep bed geometry such as in the area around the ice fall near the centre of the glacier, a comparative study has shown that ice volume variations and the retreat and advance rates of a SIA-based model are approximately the same as those computed with a full system model (Leysinger-Vieli and Guðmundsson, 2004).Since the focus of this study is not to model the ice flow in the area of the ice fall, but rather to study ice volume variations, large-scale geometry changes and the advance and retreat of Hoffellsjökull in response to changes in the climate forcing, it is appropriate to apply a SIA model.
Basal sliding was not explicitly included in the model used by Aðalgeirsdóttir et al. (2006) and Guðmundsson et al. (2009b), but implicitly included in the calibration of the power-law constitutive relationship (Glen's flow law, Glen, 1955;Cuffey and Paterson, 2010), which relates strain rate to deviatoric stresses.Aðalgeirsdóttir et al. (2006) and Guðmundsson et al. (2009b) used a series of model runs with varying flow parameters to select the parameterization The Cryosphere, 5, 961-975, 2011 www.the-cryosphere.net/5/961/2011/that best simulated the measured glacier geometry.The rate factor (A) calibrated in this manner is on the order of 6 × 10 −15 s −1 kPa −3 ; a value that has been recommended for temperate ice (Paterson, 1994), but characterizes softer ice than the average of a series of model calibrations that is now recommended (Cuffey and Paterson, 2010).Here, two approaches are tested: Method I (MI) that uses the same flow law parameter, A, as in the two studies discussed above, implicitly including basal sliding, and Method II (MII) that includes a Weertman type sliding law (Paterson, 1994) where the sliding velocity is assumed to be proportional to a power of the basal shear stress, τ b , (V slid = C • τ m b ); C is the sliding parameter and the exponent m = 3 in our calculations.When including explicit basal sliding in the model, a flow parameter that characterizes stiffer ice is required (e.g.Aðalgeirsdóttir, 2000).The available data do not allow the determination of the relative contributions of deformation and sliding velocities for the glacier.Many combinations of A and C will produce a satisfactory fit to observations.A series of model runs were carried out to obtain the flow and sliding parameters for Hoffellsjökull that resulted in a good simulation of the observed 20th century evolution of the glacier geometry.The obtained values for the rate factor and the sliding parameter are A = 4.6 × 10 −15 s −1 kPa −3 and C = 10 × 10 −15 m a −1 Pa −3 , respectively.
The ice divide is kept at a fixed location in the model computations presented here and no flow is allowed across that boundary.This is not an entirely realistic boundary condition as some ice flow may in fact occur across topographic ice divides and it can be expected that the ice divides of the eastern outlets of the Vatnajökull ice cap can shift as a consequence of the response of the ice cap to mass balance variations.As a first approximation, no ice flow across the topo-graphic ice divide is assumed (consistent with the SIA formulation of the ice flow model) and fixed boundaries of the main ice flow basins may be assumed to be reasonable since the location of the ice divide is to a large degree controlled by the basal topography.However, the assumption of fixed ice divides may be expected to become increasingly inaccurate when simulated changes in the geometry of the glacier become relatively large compared with the original size of the glacier.

Steady state experiments
Steady state experiments (model runs with constant input parameters) were carried out to (i) test the performance of the model, (ii) investigate the stability of the ice geometry and (iii) derive an appropriate initial LIA max geometry used to quantify the sensitivity of Hoffellsjökull to changes in individual climatic as well as physical parameters.The forcing and response of the glacier are in reality always varying and a steady state does not exist in nature but it is nevertheless useful to model steady states to understand the sensitivity of the model.Applying the observed ice surface geometry as initial state for transient runs may yield unrealistic ice geometry changes due to model approximations and deficiencies.
The LIA max volume of Hoffellsjökull was simulated by running a spin-up of the coupled mass balance -ice flow model (Table 3) assuming the average climate of the coldest recorded period in Iceland from 1860-1890 to remain constant (baseline period 1 in Fig. 8).The average temperature at Hólar in Hornafjörður and precipitation at Fagurhólsmýri was reconstructed to be about 1.0 • C and 0.37 m w.e. a −1 www.the-cryosphere.net/5/961/2011/The Cryosphere, 5, 961-975, 2011 Table 3. Volume (V ) and area (A) of Hoffellsjökull.
V o and A o are observed, V s MI is simulated with method MI (A = 6.8 × 10 −15 s −1 kPa −3 ) implicitly including basal sliding, V s MII is simulated with method MII (A = 4.6 × 10 −15 s −1 kPa −3 and C = 10 × 10 −15 m a −1 Pa −3 ).All volumes correspond to autumn, where the 2001 volume, measured in the spring, has been corrected with the summer balance in Table 1.Here, the error of the corrected 2001 map is estimated to be within 2 m.The errors shown are from the surface maps, calculated using the standard error formula (this error is independent between each observed year).In addition to this error, there is an uncertainty of <3 km 3 (same for all years), due to errors in the bedrock map.Those two volume errors need to be displayed and interpreted independently due to their different nature.Table 4. Mean annual temperature (T ) at Hólar in Hornafjörður and precipitation (P ) at Fagurhólsmýri, averaged over the years from t 1 to t 2 .V o is the volume at the year t 2 , from Table 3. V s is the volume of a stable spin-up glacier corresponding to the average climate over the years t 1 -t 2 .3 and the red line in Fig. 10a).This result indicates that in 2001 the glacier was still responding to the colder and/or wetter climate before 1981-2000 and is adjusting to the previous changes in the climate forcing.

Year
Model runs with the trench under the terminus area filled (horizontal dotted line in Fig. 4) resulted in the same steady state LIA max volume, indicating that the size of the trench does not have a large impact on the model simulations.All the model runs were therefore done with the measured 2001 shape of the trench.
The three steady state runs in Fig. 10a show that the glacier model responds to both temperature and precipitation variations.Assuming a warming of T = 1 • C and no precipitation change ( P = 0 m w.e. a −1 ) when shifting the climate forcing from baseline period 1 to 2 results in a 25 % smaller ice volume compared with a simulation when precipitation is also increased (red and green lines in Fig. 10a).This • C −1 d −1 ).In (A) and (C) A = 6.8 × 10 −15 s −1 kPa −3 , C = 0 (no sliding, method MI); in (A) and (B) ddf s = 4.0 mm w.e.
• C −1 d −1 ; in (B) and (C) the reference climate of 1981-2000 is kept fixed.Note: in all the model runs, the initial ice surface is the steady spin-up state corresponding to LIA max and the red curve is the same in (A), (B) and (C) (same model run).model result indicates an importance of precipitation forcing for maintaining the glaciers in this area.
The sensitivity test for the model parameters (Fig. 10b) shows that tuning the deformation and sliding parameters, A and C, does not have a unique solution as mentioned above.Many combinations of A and C give approximately the same steady state volume and model behavior.The two selected combinations, methods MI and MII, indicate that implicitly taking basal sliding into account and assuming softer ice (red line) gives a similar result as stiffer ice with explicit basal sliding (green line).
Figure 10c shows the sensitivity of the steady state ice volume to the adopted values of the degree-day coefficients for The Cryosphere, 5, 961-975, 2011 www.the-cryosphere.net/5/961/2011/snow and ice, ddf s and ddf i .The figure shows that a change in ddf s has a greater effect than a change in ddf i .The likely reason for this larger effect could be that the degree-day coefficient for snow affects larger area of the glacier, the accumulation area, and increasing the melt of snow raises the equilibrium line altitude and lowers the albedo, a positive feedback which increases the melt.

Reconstruction of the 20th century evolution
The observed 1890 glacier geometry was used as an initial configuration for a simulation through the 20th century.The year 1895 was chosen as the starting year of the coupled mass balance -ice flow model run, as it marks the beginning of warming after the coldest part of the LIA (see temperature time-series in Fig. 8).The ∼20 % volume reduction of the ice cap from 1895 to 2010 was successfully simulated by the coupled model (Fig. 11a) and the model simulates the measured volume changes during the 20th century well (Table 3).
Only a small difference in the volume evolution is found when using methods MI and MII and a reasonable agreement is obtained between the observed and modeled surface velocity field in both cases (Fig. 7b and c).The modeled velocity in the year 2002 compares well with the spatial pattern of the SPOT derived summer surface velocity (Table 2 and comparison on profile DD' in Fig. 6a), indicating that the model captures the large-scale flow pattern of the glacier.The 20th century glacier runoff anomalies, relative to a reference runoff computed at the start of the model simulations, from the area that was ice-covered at LIA max (234 km 2 ) was computed from the mass balance and precipitation fields as simulated by the model (Fig. 11b).The negative values indicate less runoff from the area than the reference value in 1895.The fastest retreat rates and highest runoff rates were obtained for the warm years between 1925-1960and 2000-2005 (Figs. 8 and 11b) (Figs. 8 and 11b).

Simulation of future response to an ensemble of climate scenarios
The coupled mass balance -ice flow model, which has been calibrated with available observations, was used to simulate the future response of the glacier during the 21st century, using the temperature and precipitation scenarios shown in Fig. 9.For comparison, a model simulation with constant climate corresponding to the last decade (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) was also carried out (∼2 • C and ∼1 • C warmer than during the baseline periods 1 and 2, respectively).According to this constant climate model run (red dashed curve in Fig. 11a), the ice volume will be reduced by 30 % with respect to the 2010 volume (Table 3) at the end of the 21st century.If the climate warms as suggested by most of the climate change scenarios, the glacier will have almost disappeared at the end of the 21st century.
The runoff anomalies relative to the 1895 runoff reference, corresponding to the same model runs are shown in Fig. 11b.The runoff from the LIA max ice-covered area is projected to increase for the next 30-40 yr and be similar or slightly larger than the runoff during the warm years between 1925 and 1960 when the ice-covered area was considerably larger.The runoff is simulated to start decreasing after year 2050 as a consequence of the reduction in the ice-covered area and approach a level near the end of this century that is similar to the runoff during cooler periods of the 20th century when the glacier surface mass balance was close to zero.
The modelled surface mass balance and the glacier extent at 50 yr intervals throughout the simulation are shown in Fig. 12 (the DMI-HIRHAM ECHAM5 scenario is used after 2010).The accumulation area shrinks to zero during the first half of the 21st century and the glacier thins and is nearly vanished at the end of the simulation.The average mass balance -elevation relationship for 10 yr periods at the same 50 yr intervals is shown in panel F in Fig. 12.Comparison with observed mass balances (Fig. 5) indicates that the degree-day mass balance model simulates the current mass balance distribution of the glacier reasonably well.

Discussion and conclusion
The volume change of the Hoffellsjökull outlet glacier during the last decade can be estimated from observations using two independent methods, (i) integrating the surface mass balance (SMB) in Table 1 over the surface area in Table 3 (assuming density of ice 900 kg m −3 to get ice equivalent volume change), and (ii) integrating the elevation difference between the DEMs over the glaciated area (dDEMs).Volume changes from these two methods for two time periods are shown in Table 5.During these periods, volume loss occurs not only by surface melt, but also by melting and calving into the frontal lake.Annual volume loss in the lake can be roughly estimated (e.g.Björnsson et al., 2001) from the incoming radiation of ∼100 W m −2 (the average sort wave radiation measured at the lower automatic weather station on Hoffellsjökull) and a lake area evolving from negligible size in 2001 to ∼2 km 2 in 2010, to be 0.005 km 3 a −1 (ice equivalent volume) during the first years, increasing to 0.02 km 3 a −1 in 2010.This estimated volume loss in the lake is added to the integrated SMB shown in Table 5.The simulated volume changes with methods MI and MII are −1.8 km 3 for the period 2001-2008 and −2.4 km 3 for the period 2001-2010.Both the simulated and the integrated SMB volume change estimates are lower than the estimates derived from the DEMs but the estimates are not inconsistent when uncertainties are taken into account.The simulations do not take into account the volume loss due to calving into the frontal lake.The simulated volume changes would change slightly with a flow model that includes volume loss from calving, but an inclusion of a calving model is beyond the scope of this study.
www.the-cryosphere.net/5/961/2011/The Cryosphere, 5, 961-975, 2011  For a given constant reference climate the same glacier volume can be reconstructed by using both method MI (rate factor A = 6.8 × 10 −15 s −1 kPa −3 , implicitly including the basal sliding) and MII (rate factor, A = 4.6 × 10 −15 s −1 kPa −3 , together with a Weertman type sliding law with a constant sliding parameter, C = 10 × 10 −15 m a −1 Pa −3 ) (Fig. 10b).The same is true for a simulation of the time evolution from 1895 to 2010 (blue and red solid curves in Fig. 11a; Table 3).Using the flow parameter corresponding to stiff ice without adding the sliding leads to a steady state with about ∼10 % larger ice volume (Fig. 9b).The tuning of the deformation and sliding parameters, A and C, is non-unique and without better constraints for the viscosity of the ice and basal condition it is not possible to determine the relative contributions of internal deformation and basal sliding to the observed surface velocity (Fig. 7; Table 2) and select unique values for the flow and sliding parameters.
Model studies of the recovery of depressions formed in the surface of Vatnajökull during a subglacial eruption and subsequent outburst flood in 1996 have been carried out to infer the possible range of the flow parameter, A (Aðalgeirsdóttir et al., 2000;Jarosch and Guðmundsson, 2007).The results indicate that about 1/3 lower value for A than has been recommended for temperate ice by Paterson (1994) gives the best fit to the surface velocity measurements.These results are in line with other model studies and are reflected in an update for the recommended flow parameter presented by Cuffey and Paterson (2010).Information about the relative importance of basal sliding may be inferred from seasonal variations in the ice surface velocity.The SPOT velocity maps show 30-50 % faster flow within the ablation zone during late summer than the annual average (Figs.6, 7, Table 2).The velocity measurements at stake Hof (∼1120 m a.s.l.) shown in Fig. 6b show that the early summer velocity is The Cryosphere, 5, 961-975, 2011 www.the-cryosphere.net/5/961/2011/The sensitivity of the model results to the applied degreeday factors was tested by varying the values in the steady state experiments (Fig. 10c).The degree-day factors determined by using mass balance data from the two stakes on Hoffellsjökull and data from the two automatic weather stations are similar to the degree-day factors determined by optimization of 23 mass balance stakes on S-Vatnajökull.The sensitivity tests show that varying the degree-day factors by the standard deviation obtained from an optimization for each year separately does not have a major impact on the steady state volume (Fig. 10c).
The advantage of the degree-day melt model is that it uses only temperature and precipitation away from the glacier as an input (the only available long term observation).Detailed comparison of the degree-day melt model and a more ad-vanced energy balance model have been carried out on both Langjökull (Guðmundsson et al., 2009a) and Vatnajökull ice cap (Guðmundsson et al., 2003).They found the degree-day parameters to be relatively constant with time within the ablation season and with elevation, varying mainly with the surface conditions, that is whether the melting surface is snow or ice.This was found in spite of the fact that the relative importance of the different energy balance components varies substantially during the ablation season.The energy balance at the glacier surface may be altered in a warmer climate (e.g.Björnsson et al., 2005;Guðmundsson et al., 2009a), e.g.due to increased long-wave radiation, stronger glacier winds, earlier start of melting so that the highest incident radiation will fall on a surface with an already reduced albedo, and ablation extending into the autumn resulting in melt driven by periods of high temperatures and strong winds.The results of Guðmundsson et al. (2009a) indicate, however, that given no extreme changes in surface albedo values, the degree-day model may provide reasonable predictions of increased ablation in response to regional temperature changes within ±3 • C from a given reference period for which the degreeday parameters are optimized.In our study the temperature changes are within this range and therefore it is concluded that the optimized degree-day parameters are appropriate for the simulations with the coupled surface mass balance -ice flow model that are presented here.
Runoff from the area covered by the LIA max glacier for the whole modeled period 1895-2100 is shown in Fig. 11b.The runoff accounts for both melting of snow and ice and the precipitation in the form of rain in the area.The results show large interannual variability as well as an increase in runoff with temperature in the first half of the 21st century while the volume of the glacier is declining.It should be noted that this model output cannot be validated against observations as no river discharge measurements are available from this area.Towards the end of the simulation, the runoff is reduced because the glacier has nearly disappeared.Even though there is a large spread in the simulated results depending on the applied scenario, the general trend of an initial increase, followed by a decrease in runoff during the second half of the 21st century is similar for most of the scenarios.It is found that until ∼2030, the large interannual variations in the climate lead to interannual runoff variations of a similar magnitude as the average runoff increase with respect to the period 1981-2000.Around the middle of the century, most of the climate change scenarios indicate that the simulated runoff changes due to the warming of the long-term climate have become larger than the magnitude of the interannual variations.
The simulated volume and runoff changes of the Hoffellsjökull outlet glacier are similar to the results for the larger Hofsjökull and Langjökull ice caps that have been forced with the same climate scenarios (Jóhannesson et al., 2011).The future of the Icelandic glaciers and ice caps is dependent on their elevation range (e.g.Guðmundsson et al., www.the-cryosphere.net/5/961/2011/The Cryosphere, 5, 961-975, 2011Cryosphere, 5, 961-975, 2009b) ) and the future climate.If the current climate will persist into the future the glaciers are predicted to lose about 20-30 % of their volume, depending on their elevation range.
In response to the predicted climatic changes, with a substantial increase in average temperatures, most of the glaciers will almost vanish in 100-200 yr.

Fig. 1 .
Fig. 1. (A) Iceland and the largest ice caps, Vatnajökull (Va), Hofsjökull (Ho) and Langjökull (La).Locations of the weather stations, used to reconstruct the temperature and precipitation records, are shown with letters: Reykjavík (R), Fagurhólmýri (F), Haell (H), Stykkishólmur (S), Teigarhorn (T), Vestamannaeyjar (V), Akureyri (A) and Hólar in Hornafjörður (HH).(B) The surface topography of Vatnajökull ice cap.Dots show the sites of mass balance and velocity measurements, blue dots show the location of the mass balance stakes at Breiðamerkurjökull (Bre).The red box indicates the position of the frame to the right.(C) Hoffellsjökull surface topography in 2001.The ice divide and the model domain are indicated with the red curve enclosing a glaciated area of ∼212 km 2 in 2001.Black triangles show the locations of automatic weather stations on the glacier.N is the location of the central nunatak Nýju Núpar.

Fig. 2 .
Fig. 2. (A) Measured bedrock topography of Hoffellsjökull (2001).Blue colours indicate elevation below sea level.(B-E) Surface topography at different times, showing retreat and thinning during the 20th century.The radio-echo sounding (RES) line network and sites of RES point survey is shown in (E).The red line is the 2001 ice divide for Hoffellsjökull (212 km 2 ).

Fig. 3 .
Fig. 3.A topographic relief shading showing the August 2010 Li-DAR DEM of the terminus and the lower part of Hoffellsjökull.The LIA terminus moraines can be seen in front of the terminus as well as the extensive break-up of the eastern branch of the terminus into the lake in front of the glacier that occurred in 2010.

Fig. 4 .
Fig. 4. Two longitudinal profiles and two cross sections showing the thickness and the location of the terminus at different times.The map shows the location of the sections and the terminus position at the different times.

Fig. 5 .
Fig. 5. Winter, summer and annual balance measured at Hoffellsjökull and Breiðamerkurjökull.Measurements from the 1930s expedition (at the lowest stakes only annual balance was measured and therefore winter balance set to zero) are shown with black symbols.The red dots show the average and standard deviation of the mass balance measurements in the period 2001-2010 at the two stakes on Hoffellsjökull and the two stakes close to the ice divide.The blue dots are similar measurements from Breiðamerkurjökull that have been used to determine the elevation dependency of mass balance when computing the mean specific mass balance in Table 1 (see Fig. 1 for location of the stakes).

Fig. 7 .
Fig. 7. (A) Late summer ice surface velocity determined by crosscorrelation of 2.5 m resolution SPOT5 HRG satellite images acquired on 27 August and 22 September 2002.The numbers 1-3 point out locations where annual surface velocity signal could be deduced by correlating two SPOT5 HRG satellite images from late September 2002 and 2003 (see Table 2).(B-C) Depth-averaged modeled ice velocity corresponding to the simulated 2002 ice surface geometry, using A = 6.8 × 10 −15 s −1 kPa −3 implicitly including the basal sliding (in B), and using A = 4.6 × 10 −15 s −1 kPa −3 and C = 10 × 10 −15 m a −1 Pa −3 (in (C)).(D) The contribution of the modeled basal sliding to the velocity shown in (C).

Fig. 8 .
Fig. 8. Mean annual temperature (upper panel) at Hólar in Hornafjörður and precipitation (lower panel) at Fagurhólsmýri.The gray areas indicate the periods when the time-series were reconstructed.The horizontal lines indicate the average climate for the two baseline periods 1860-1890 and 1981-2000.The dots show the temperature and precipitation in 1895, the first year of the model simulations.

Fig. 9 .
Fig. 9. Simulated mean annual temperature (upper panel) and precipitation (lower panel) development averaged over whole of Iceland; scenarios from the Nordic CES Project.The average climate of 2000-2009 is plotted for comparison (dashed line).The DMI HIRHAM ECHAM5 climate scenario is shown in bold because the results from this run are shown in bold in Fig. 11 and the modelled surface mass balance and glacier extent from this run are shown in Fig. 12.
, during the period 1860-1890 than in 1981-2000 (baseline period 2 in Fig. 8, Table 4) when the surface balance of most of the Icelandic ice caps was near zero, as mentioned earlier.The model was then forced from the simulated LIA max configuration with a step change in temperature and precipitation corresponding to the difference between the periods 1860-1890 and 1981-2000 towards a new steady state.The steady state volume corresponding to the period 1981-2000 was found to be ∼10 % smaller than the observed volume in 2001 (Table

Fig. 10 .
Fig. 10.(A) Steady state volume of the model glacier corresponding to the 1860-1890 baseline climate ( T = 0 • C & P = 0 m a −1 ; Table 4 and Fig. 6) and sensitivity of the model to step changes in T and P (the 1981-2000 baseline climate corresponds to T = 1 • C & P = 0.37 m a −1 ).The dashed straight line shows the 2001 measured volume (57.0 km 3 ).(B) Sensitivity to the rate factor A (in 10 −15 s −1 kPa −3 ) and the sliding parameter C (in 10 −15 m a −1 Pa −3 ).(C) Sensitivity to shift by one σ of ddf s and ddf i (in mm w.e.• C −1 d −1 ).In (A) and (C) A = 6.8 × 10 −15 s −1 kPa −3 , C = 0 (no sliding, method MI); in (A) and (B) ddf s = 4.0 mm w.e.• C −1 d −1 and ddf i = 5.3 mm w.e.• C −1 d −1 ; in (B) and (C) the reference climate of 1981-2000 is kept fixed.Note: in all the model runs, the initial ice surface is the steady spin-up state corresponding to LIA max and the red curve is the same in (A), (B) and (C) (same model run).

Fig. 11 .
Fig. 11.(A) Simulated evolution of Hoffellsjökull ice volume during the 20th and 21st centuries, initiated with the observed LIA max glacier geometry.The green dots are volume estimates from DEMs.The volume change 1895-2010 is simulated by using the T and P records in Fig. 6 and the two methods MI (blue curve) and MII (red curve) for implicit and Weertman type sliding, respectively.The future evolution is computed with MI by using the climate scenarios in Fig. 8 and by maintaining the average climate of 2000-2009 (dashed red curve).(B) Simulated specific runoff anomalies (precipitation and glacier melt) relative to the specific runoff at the start of the simulation (2.5 m a −1 in 1895) from the area covered by the LIA max glacier (∼234 km 2 ).Negative values indicate smaller runoff than at the start of the simulation.

Fig. 12 .
Fig. 12. (A-E)The modeled surface mass balance at 50 yr intervals computed with method MI and the DMI HIRHAM ECHAM5 climate scenario (shown with bold line in Fig.9).The average modeled surface mass balance for 10 yr periods corresponding to the first decade of each 50 yr interval is plotted as a function of elevation in (F) (location of the profile is shown in (A).

Table 5 .
Ice equivalent volume changes for two periods(2001- 2008 and 2001-2010)estimated from observed surface mass balance and volume loss in the frontal lake, from the difference in the measured DEMs and simulated with the coupled surface mass balance-ice flow model.