The importance of accurate glacier albedo for estimates of surface mass balance on Vatnajökull : evaluating the surface energy budget in a regional climate model with automatic weather station observations

A simulation of the surface climate of Vatnajökull ice cap, Iceland, carried out with the regional climate model HIRHAM5 for the period 1980–2014, is used to estimate the evolution of the glacier surface mass balance (SMB). This simulation uses a new snow albedo parameterization that allows albedo to exponentially decay with time and is surface temperature dependent. The albedo scheme utilizes a new background map of the ice albedo created from observed MODIS data. The simulation is evaluated against observed daily values of weather parameters from five automatic weather stations (AWSs) from the period 2001–2014, as well as in situ SMB measurements from the period 1995– 2014. The model agrees well with observations at the AWS sites, albeit with a general underestimation of the net radiation. This is due to an underestimation of the incoming radiation and a general overestimation of the albedo. The average modelled albedo is overestimated in the ablation zone, which we attribute to an overestimation of the thickness of the snow layer and not taking the surface darkening from dirt and volcanic ash deposition during dust storms and volcanic eruptions into account. A comparison with the specific summer, winter, and net mass balance for the whole of Vatnajökull (1995–2014) shows a good overall fit during the summer, with a small mass balance underestimation of 0.04 m w.e. on average, whereas the winter mass balance is overestimated by on average 0.5 m w.e. due to too large precipitation at the highest areas of the ice cap. A simple correction of the accumulation at the highest points of the glacier reduces this to 0.15 m w.e. Here, we use HIRHAM5 to simulate the evolution of the SMB of Vatnajökull for the period 1981–2014 and show that the model provides a reasonable representation of the SMB for this period. However, a major source of uncertainty in the representation of the SMB is the representation of the albedo, and processes currently not accounted for in RCMs, such as dust storms, are an important source of uncertainty in estimates of snow melt rate.


Introduction
Worldwide, glaciers and ice caps are losing mass at increasing rates as a response to climate change (e.g.Vaughan et al., 2013).Major changes in the dimensions of glaciers are expected to affect the sea level and climate throughout the world, and it is therefore important to describe and understand the glacier climate.Glacier retreat and mass loss at significantly increasing rates are also observed for Icelandic glaciers (Björnsson et al., 2013), which could potentially contribute to the rise in sea level by 1 cm (Björnsson and Pálsson, 2008;Björnsson et al., 2013).The runoff from Vatnajökull ice cap is economically important to hy-Published by Copernicus Publications on behalf of the European Geosciences Union.
dropower production in Iceland and the present and future mass balance is thus of keen interest.Numerical highresolution regional climate models (RCMs), such as MAR (Gallée and Schayes, 1994), RACMO2 (Meijgaard et al., 2008), or HIRHAM5 (Christensen et al., 2006), are valuable tools for estimating the meteorological parameters and mass balance variability at the surface of glaciers.However, to carry out reliable future projections, or reconstruct the past climate, it is important to evaluate how well models simulate the present climate Evaluation of RCMs is important, not only because it reveals possible biases in the model but also because it could yield recommendations for model improvements.Much work has gone into evaluating RCMs over Greenland (e.g.Box and Rinke, 2003;Noël et al., 2015;Rae et al., 2012;Langen et al., 2017;Fettweis et al., 2017) and Antarctica (e.g.Lenaerts and Van Den Broeke, 2012;Agosta et al., 2015), but less effort has gone into evaluating them over Iceland (e.g.Ágústsson et al., 2013;Nawri, 2014).
However, a long-term meteorological monitoring programme has been conducted on Icelandic glaciers since the 1991-1992 glaciological year (e.g.Björnsson et al., 1998).Therefore, Icelandic glaciers are excellent candidates for evaluating modelled meteorological and SMB components.Compared to Greenland, observations are recorded in a relatively small area, offering a good opportunity to evaluate the spatial and temporal variability of the HIRHAM5 model on a regional scale.As albedo in Iceland is significantly different from that of Greenland or Antarctica, e.g.due to frequent dust storms and occasional volcanic eruptions, model evaluations over Iceland provides important insight into the effect of albedo changes on the glacier energy balance.
Due to the large spatial and temporal variation in albedo of Icelandic glaciers (spanning from less than 0.1 for dirty ice in the ablation zone to 0.9-0.95 for new snow), and the large sensitivity of melt to variations in albedo, it is crucial to have correct estimates of the albedo when modelling the surface mass balance.However, accurate modelling of the albedo can be challenging.For example, volcanic eruptions and dust storms can significantly lower the glacier albedo, and thus increase the amount of melt (e.g.Conway et al., 1996;Gascoin et al., 2017;Wittmann et al., 2017), but are difficult to include in albedo models.Accurate simulations of the ice albedo is also problematic, as for some glaciers it varies with elevation (e.g.Knap et al., 1999) but not for others (e.g.Greuell et al., 1997).In addition, the ice albedo may decrease with time (e.g.Reijmer et al., 1999), increase with time (e.g.Oerlemans and Knap, 1998), or remain constant (e.g.Greuell et al., 1997) depending on the glacier.
Here we present a 1981-2014 SMB data set of Vatnajökull ice cap modelled by HIRHAM5 at 5.5 km resolution.HIRHAM5 is a state-of-the-art, high-resolution RCM that has been well validated over Greenland (e.g.Box and Rinke, 2003;Lucas-Picher et al., 2012;Rae et al., 2012;Langen et al., 2017).In this study, HIRHAM5 incorporates an up-dated albedo scheme, using a background MODIS ice albedo field, in the aim of capturing the effect of dust and tephra on ice albedo in the ablation zone.This method of determining the ice albedo has previously been used by, for example, van Angelen et al. (2012).Model simulation results are compared to observations from automatic weather stations (AWSs) and in situ mass balance observations, in an effort to improve the performance of the model.The possible physical reasons for any model biases are discussed, and recommendations for corrections are made where possible.

HIRHAM5
In this study we employed the RCM HIRHAM5 (Christensen et al., 2006), which was developed at the Danish Meteorological Institute.It is a hydrostatic RCM which combines the dynamical core of the HIRLAM7 numerical forecasting model (Eerola, 2006) and physics schemes from the ECHAM5 general circulation model (Roeckner et al., 2003).Model simulations have been successfully validated over Greenland using AWS and ice core data (e.g.Box and Rinke, 2003;Stendel et al., 2008;Lucas-Picher et al., 2012;Langen et al., 2015;Rae et al., 2012;Langen et al., 2017).
While the original HIRHAM5, as described in Christensen et al. (2006), used unchanged ECHAM physics, an updated model version, which includes a dynamic surface scheme that explicitly calculates the surface mass budget on the surface of glaciers and ice sheets, is used in this study.This new scheme takes melting of snow and bare ice into account and resolves the retention and refreezing of liquid water in the snow pack (Langen et al., 2015(Langen et al., , 2017)).In addition, the fivelayer surface scheme in ECHAM has been expanded to 25 layers.

New albedo parametrization
The updated model also features a more sophisticated snow albedo scheme (Nielsen-Englyst, 2015) than that used in the original HIRHAM5; whereas the previous scheme was purely temperature dependent, the new scheme depends both on the age of the snow and the surface temperature.The scheme is similar to that used in Oerlemans and Knap (1998), which assumes that the albedo decays exponentially as it ages, but in this study an additional temperature component is applied.If there is snow on the surface, the change in the snow albedo from one time step to the next depends on whether the surface is in a dry (< 271 K) or wet regime (≥ 271 K).In the dry regime, the surface temperature is too low for any melting to occur, while in the wet regime the temperature in the surface layer is high enough for the surface to be melting.The snow albedo changes over a time step, δt, as where α mx is the minimum snow albedo value that can be reached from ageing of the snow and τ x is a timescale which determines how fast the albedo reaches its minimum value.These two variables take on different values depending on whether the snow is in the dry (d) or wet (w) regime.
Observations from the AC and ELA stations were used to determine α mx and τ x .The optimal variables were found by minimizing the weighted mean RMSE between the modelled and measured albedo by varying the values of α mx and τ x .The best-fit values were found to be α md = 0.65, α mw =0.41, τ md = 5 days, and τ mw = 10 days.
Albedo is only refreshed to the maximum value if snowfall constitutes more than 95 % of the total precipitation.A partial refreshment is possible as the albedo is only reset to the maximum allowed value if the amount of snowfall on that day (S 0 ) is higher than 0.03 m w.e.This threshold was chosen to provide the best fit with the AWS observations.The rate of refreshment b is given by where S f is the amount of snowfall during the model time step in m w.e. and S 0 is the critical amount of snowfall in m w.e. per model time step needed to completely refresh the albedo.Using this rate, the albedo is then refreshed using where α max is the maximum albedo for freshly fallen snow, set equal to 0.85 as this provides the best average fit with the observations.In the case of shallow snow cover, the surface albedo will be affected by the albedo of the underlying ice.A smooth transition between the snow and bare ice albedo is therefore implemented, and the final albedo is thus expressed as where d is the snow depth, and d s is a characteristic scale for snow depth.Following Oerlemans and Knap (1998), the characteristic scale is set to 3.2 cm snow depth.If no snow is present, the albedo is set to the bare ice albedo.The bare ice albedo is determined from a background ice albedo map which was created using MODIS observations from the period 2001-2012.How this map was created is described in Sect.3. The extent to which this bare ice MODIS albedo map improves the simulations will be estimated by comparing the results with those from a model simulation using a constant ice albedo in Sect.4.8.

Experimental design
In this study, HIRHAM5 is run at a resolution of 0.05 • (equivalent to ∼ 5.5 km) on a rotated pole grid for the period 1980-2014.The model uses 31 irregularly spaced vertical atmospheric levels from the surface to 10 hPa with a model time step of 90 s in the dynamical scheme.The model is configured for a domain containing all of Greenland and Iceland.The model is forced at the lateral and lower boundaries by the ECMWF ERA-Interim reanalysis data set (Dee et al., 2011), which uses observations from satellites, weather balloons, and ground stations to create a comprehensive reanalysis of the atmosphere.The model is forced by temperature, wind, relative humidity, and surface pressure at the lateral boundary, and sea surface temperature and sea ice fraction at the lower boundary at 6 h intervals.
The new snow/ice surface scheme discussed above is run offline in this study, meaning that the subsurface scheme is run separately from the atmospheric code.This is done by forcing the subsurface scheme every 6 h by radiative and turbulent surface fluxes, as well as snow, rain, evaporation, and sublimation data from a HIRHAM5 experiment (Mottram et al., 2016) with a previous version of the albedo and refreezing schemes (e.g.Langen et al., 2017).While a full, highresolution HIRHAM5 run is computationally very expensive, the offline model offers a fast and flexible option to test new model implementations and allows for a quick and thorough spin-up of the subsurface.The offline model was initialized with values from a previous offline model run with a different albedo scheme and then a model spin-up was performed by integrating the model for 150 years repeating the forcing from 1980.The largest adjustments occurred during the first 75 years of the spin-up, after which the variation was much smaller than the interannual variability.At the end of the run, the solar radiation, surface mass balance, runoff, snow depth, and refreezing had all converged, as had the temperature, liquid and snow content in all 25 subsurface layers.The final state of the spin-up was then used as the initial condition for the 1980-2014 model simulation.The reported values of albedo, upward longwave and shortwave radiation, and surface mass balance in the following are all from the offline run.
A disadvantage of this method is that it neglects feedbacks between the atmospheric circulation and the surface conditions like the albedo and temperature.However, since the surface temperature of Vatnajökull is typically near the melting point during the summer, both in reality and in the model, changes in the albedo should not have a large effect on upward longwave radiation and the turbulent fluxes.Thus, while the updated surface scheme is important for the mass balance components, the error due to the neglected feedbacks is likely small in the model calculations.

Model uncertainty
Due to nonlinearities in the HIRHAM5's model dynamics and physics, it has an implicit uncertainty due to internal model variability originating from nonlinear processes (e.g.Giorgi and Bi, 2000;de Elía et al., 2002).This variability is caused by numerical sensitivity, uncertainty in the boundary and initial conditions, and errors due to model parametrizations (e.g.Box and Rinke, 2003), including, for example, the albedo parameterization, the vertical gradients in the boundary layer, or cloud radiative effects.In addition, using a constant value of z 0 for both snow and bare ice could lead to large errors in the turbulent fluxes (e.g.Brock et al., 2000).

Observational data
The primary observational data set used in this study was collected by AWSs at selected locations on Vatnajökull.Since 1994, 1-13 stations have been operated on the ice cap during the summer months (e.g.Oerlemans et al., 1999;Guðmundsson et al., 2006).The temperature, relative humidity, wind speed, and wind direction at 2 m above the surface have been measured during the entire period (1992-present), while the radiation components have been measured since 1996.For this study, data from five AWSs were considered -three on Brúarjökull (B) and two on Tungnaárjökull (T) (see Fig. 1).The uncertainties of the AWS observations vary depending on the sensor.The temperature and humidity sensors have an accuracy of 0.2 K and 2 % for temperature and humidity, respectively, while the accuracy of the wind speed is 0.2 m s −1 (Guðmundsson et al., 2009).The radiative fluxes were measured using either Kipp and Zonen CM14, CNR1 or CNR4 sensors that have a maximum manufacturer-reported uncertainty of ±10 % for daily totals (e.g.Kipp and Zonen, 2002).However, the uncertainty has independently been evaluated to be lower (3-5 %) when used in an ice sheet environment (van den Broeke et al., 2004;Guðmundsson et al., 2009).The turbulent fluxes, combining sensible and latent heat fluxes, and surface pressure were not measured at the stations, but were estimated using the methods described in Sect.3.1.
In addition to AWS data, in situ mass balance measurements were used to evaluate the simulated surface mass balance (SMB) at several sites on Vatnajökull.Conventional in situ mass balance measurements have been carried out every glaciological year since 1991-1992, with 60 stations measured each year on average.The measurement sites are shown in Fig. 1.The uncertainty of the mass balance measurements has been estimated to be ±0.3 m w.e.
The SMB measurements are conducted at the beginning and end of the accumulation season in order to measure both the winter and summer balance.The winter balance is measured in the beginning of the melt season by drilling down to the previous summer layer and weighting the snow column.The summer surface is used as the reference level even if some snow accumulation had occurred by the time the summer balance measurements were conducted.The snow thickness on top of the summer surface at the time of the autumn survey has been measured since 1995.This is needed when comparing with the simulation of snow accumulation.
Observations of the broadband albedo in the shortwave spectrum (0.3-5.0 µm) from the MODerate Resolution Imaging Spectroradiometer (MODIS) were used to create a background map of the ice albedo at all glacier grid points in HIRHAM5, which was used in the implemented HIRHAM5 albedo scheme.MODIS product MCD43A3 v006 was used for the background map (Schaaf, 2015).The MODIS estimates of the albedo on Vatnajökull are in good agreement with AWS data (Gascoin et al., 2017).The MODIS data were extracted in geographical coordinates (long-lat) at a resolution of 0.005 • , i.e. close to the original MODIS resolution of 500 m.This was done using the MODIS reprojection tool with the bilinear interpolation method.These MODIS data in latitude-longitude coordinates were then resampled to match the rotated HIRHAM5 long-lat grid coordinates by bilinear interpolation using MATLAB's interpn function (MATLAB, 2015).
In order to determine the bare ice albedo at each grid point, daily MODIS data over Iceland from the period 2001-2012 were used.Years with volcanic eruptions were discarded, as the volcanic ash lowered the albedo values far below the average.The minimum autumn albedo value was then determined in each grid point using values from July-September and that value used to create a bare ice albedo map of the glaciers.The final albedo map had ice albedo values in the range 0.03-0.3for Vatnajökull.The spectral properties of ice in the ablation zone are controlled by tephra layers in the ice, which are exposed as the glacier melts (Larsen et al., 1996).Additional tephra or dust deposition will therefore only have a small effect on the spectral properties of the ice, as the ice surface is already covered in dark bands.In addition, field observations suggest that the new particles are generally washed off from year to year.Applying one background map for the entire period should therefore provide the same results as applying a map created for each year.In addition, it allows us to run the model for years where no MODIS observations are available or where the amount of observations over the ice cap are sparse due to, for example, clouds.

AWS point models
The turbulent energy fluxes were calculated from AWS measurements using a one-level eddy flux model (Björnsson, 1972;Guðmundsson et al., 2009) which uses Monin-Obukhov similarity theory (Monin and Obukhov, 1954) and implements different roughness lengths for the vertical profiles of wind, temperature, and water vapour (Andreas, 1987).The model is described in detail in Guðmundsson et al. (2009).Uncertainties of this model for example pertain to the aerodynamic roughness length for momentum z 0 .The majority of z 0 values recorded over melting glacier surfaces vary over 2 orders of magnitude (between 1 and 10 mm), but over fresh snow or smooth ice surfaces the roughness length is generally around 0.1 mm (Brock et al., 2006).An order of magnitude increase in z 0 can more than double the estimated turbulent fluxes (Brock et al., 2000), so the chosen roughness length parametrization can greatly affect the performance of the model.Generally, a constant value of z 0 is prescribed for snow and/or ice surfaces (Brock et al., 2006), which is an oversimplification as the roughness may vary significantly over the ablation season (e.g.Grainger and Lister, 1966).
However, since measurements of the evolution of z 0 over the entire measurement period are not available, a constant roughness length of 1 mm was chosen in the calculation of the non-radiative fluxes.Sensitivity tests were conducted to estimate how large an error this choice of roughness length could lead to at the used AWS sites.A roughness length of 0.1 mm would decrease the calculated turbulent fluxes by 16-22 %, while using a roughness length of 10 mm would increase the calculated fluxes by 10-19 %, depending on the station.Since the contribution of the turbulent fluxes to the total energy balance is generally low, this translates into an increase or a decrease in the total energy balance at the stations by a maximum of 7 %.
The surface air pressure at the station is also needed to calculate the turbulent fluxes, but it is not measured at the AWS sites.Instead it is estimated at the relevant elevation h using synoptic observations from meteorological stations operated by the Icelandic Met Office and the following relationship: where P (h 0 ) and T (h 0 ) are the air pressure and air temperature, respectively, observed at an elevation h 0 (e.g Wallace et al., 2006).This method has previously been applied successfully at various locations on Vatnajökull and Langjökull (e.g.Guðmundsson et al., 2006Guðmundsson et al., , 2009)).Comparisons between station values and model values are made by bilinearly interpolating the model output to the measurement position using the four closest model grid points and using only glacier-surface type grid cells.

Validation method
In order to remove the effect of seasonally varying magnitudes of the energy balance components, the percentage errors listed in Tables 2-4 are calculated as the root mean square error (RMSE) divided by the observations.
HIRHAM5 uses an elevation model over Iceland which has been interpolated onto the 5.5 km model grid.Since errors in the elevation of the glacier surface can introduce significant biases in temperature and pressure which are not caused by physical model errors (Box and Rinke, 2003), any elevation bias in the model has to be taken into account before validating the results.The elevation bias was calculated as the difference between the model elevation and GPS observations at each site (Table 1).
The temperature was corrected for the elevation bias in order to compare the model results to the AWS measurements at AWS locations.This was done using a constant lapse rate of 6.5 K km −1 , which resulted in temperature corrections on the order of 0.1-0.3K. Pressure is corrected using Eq. ( 5) decreasing the bias down to 0.1 to 0.5 hPa.Thus, although the HIRHAM5 elevation is consistently overestimated, the resulting differences are not large enough to introduce significant biases in temperature and surface pressure.

Meteorological variables
As the sensible and latent heat fluxes are computed using the surface pressure p sl , air temperature at 2 m, T 2 m , relative humidity r 2 m , and wind speed u, these model variables were evaluated at all five stations at the measurement height.How well these variables are simulated should indicate the model's ability to simulate the turbulent fluxes.
The comparison of modelled and observed mean daily values during the summer months from the period 2001-2014 is shown in Fig. 2 and Table 2.The surface pressure, p sl , which was not observed at the stations but estimated using Eq. ( 5), is generally forecast with only a small error.At each station there is a high positive correlation (r > 0.9) between modelled and estimated pressure (Eq.5), for the entire time series and for each individual year.
The model also captures the 2 m temperatures, T 2 m , satisfactorily.The largest deviation from the observations is found at the B AB station, which underestimates the temperature by 0.8 K on average.The temperature is also underestimated at the four other stations, but by at most 0.6 K.The model simulates the variation in temperature well; for example, it captures the temperature dampening over a melting glacier surface.This is expressed in the high correlation values for all five stations (r ∼ 0.9).
The measured relative humidity, r 2 m , at all five stations is generally high, with only 1-3 % of the data points at each station falling below 70 %, and the minimum daily value between 42 and 58 %.The model simulates a lower mean humidity than the measured at all five stations, with 8-20 % of the points at each stations having values lower than 70 % and minimum daily values between 18 and 30 %.Since the exchange coefficient for moisture is a function of the atmospheric temperature profile, the underestimation of the relative humidity could be due to a too low temperature gradient between the atmosphere and the surface.This is consistent with the underestimation found in the 2 m temperature.The correlation of between 0.68 and 0.7 indicates that the model simulates the humidity fluctuations satisfactorily.
The lowest wind speed level in HIRHAM5 is at 10 m and the AWS wind speeds are measured at between 2 and 4 m, depending on the year, the HIRHAM5 wind speed is extrapolated to the measurement height using a logarithmic profile with a roughness length of 1 mm.At all five locations, HIRHAM5 simulates winds that are too weak on average.This could be due to the uncertainty arising from the interpolation of the model winds from second-lowest level (30 m) to the lowest level (10 m) under stable conditions, as the wind speed can change significantly over the 20 m interval.

Longwave radiation
As shown above, HIRHAM5 underestimates the temperature at all five stations, with the largest underestimation at the B AB station.As a result, a similar underestimation of incoming longwave radiation is obtained at all five stations, with the largest difference occurring at the B AB station 3. The average percentage difference is approximately 8 % for all five locations (see Table 3), and falls well within the 10 % uncertainty of the AWS observations.However, Fig. 3a also shows that 25-30 % of the simulated days have errors larger than 10 %.
The incoming LW radiation is mainly emitted from clouds and atmospheric greenhouse gases, and therefore a source of the underestimation could be that the model underrates cloud formation and/or simulates clouds that are too optically thin in the LW region of the spectrum.An underestimation of the temperature in the atmosphere could also be causing the underestimation.
Figure 3b shows the comparison of the modelled and measured outgoing LW radiation.There is a small overestimation at the T AC station, and a small underestimation of the other four stations, but in general the model reproduces the daily values well (r ∼ 0.76).The average percentage deviation between the modelled and measured values is only around 3 %, combined with between 0.5 and 2 % of the HIRHAM5 data points having deviations larger than 10 %.
Due to an underestimation of the incoming LW radiation, and only small negative or positive biases in the outgoing LW, the net LW (incoming-outgoing) radiation has a mean negative bias at all AWS locations (−7.9 W m −2 ).

Shortwave radiation and albedo
Figure 4 and Table 3 show the comparisons of the modelled and measured components of the shortwave (SW) radiation as well as the surface albedo.On average, the incoming SW radiation is underestimated at all five stations.This underestimation is also present in the means at all five stations for most years, except in 2002,2004,2005, and 2014 at the B AB station.This suggests that there are errors in either the modelling of the clouds, e.g.due to an overestimation of the cloud fraction, the amount of cloud formation, or the optical thickness of the clouds in the shortwave region, and/or because of errors in the clear-sky fluxes.
The albedo comparison is shown in Fig. 4b.The modelled albedo at the two AB stations has the largest deviation from the observations; this is partly due to the modelled snow cover, which either does not completely disappear or disappears later in the year than the AWS data show.At the B AB station, the ice layer is generally exposed in the model (except in 2001 and 2011-2013), although the snow cover always persists longer than in reality.One exception occurs in 2001, where the modelled albedo never drops down to the ice value, whereas observations show albedo values as low as 0.03.This one year therefore highly contributes to the aver-age overestimation of the albedo.This very low albedo value could be due to a layer of dust or tephra beneath the station, so it may not represent the ice albedo.However, very low ice albedo values down to 0.05 are not uncommon in the ablation zone of Vatnajökull (e.g.Gascoin et al., 2017).Comparisons with the mass balance measurements (discussed in Sect.4.6.1)show that the winter balance is overestimated during approximately half of the measured years, which contributes to delay the albedo drop in the model.
At the T AB station, a too-thick modelled snow cover in winter is also the cause of some of the discrepancy.Comparisons with mass balance measurements (Sect.4.6.1)show that the winter balance is always overestimated at this station.An overestimation of the snow thickness at the beginning of summer, combined with an underestimation in the radiation and turbulent fluxes, leads to persistent snow cover at the end of summer.As a result, the ice surface is never exposed in the model during any of the modelled years, and the albedo never drops much below 0.4 (the minimum snow albedo), even though the AWS data shows that the ice surface was exposed during all but two years, i.e. 2008 and 2010.During these two years, the simulated albedo fits well with observations.
Another issue which affects both stations is that the MODIS albedo at these points is not as low as the measured albedo.The MODIS ice albedo at these stations is 0.10 (B AB ) and 0.16 (T AB ), whereas the observations show the albedo can drop as low as 0.01 at both stations.The albedo drops below the MODIS value every year at the B AB , and during 2001-2005 and 2011 at the T AB stations.This is presumably due to the heterogeneity of the albedo in the ablation zone, which means that a low in situ albedo value at a point cannot be captured at the current HIRHAM5 resolution.At the ELA station, the mean albedo value is underestimated (Table 3).Close to the equilibrium line, the albedo is highly variable both temporally and spatially; for example, there is a large difference in albedo depending on whether the previous year's summer surface was exposed or not.In general, the model overestimates the albedo during years where the summer surface was exposed, and underestimates the albedo during years where it was not.In addition, the winter mass balance at this station is always underestimated (Sect.4.6.1),meaning that the thickness of snow layer in spring is underestimated and the effect of the underlying ice layer will therefore be overestimated, leading to the underestimation in albedo.
The smallest difference between modelled and observed albedo is found at the two AC stations.The B AC station generally provides the best fit with the observations, while the model tends to underestimate the albedo at the T AC station.An exception to this is found in 2010 and 2011, where the albedo was overestimated by the model at both stations due to ash deposition from the Eyjafjallajökull and Grímsvötn eruptions (e.g.Gudmundsson et al., 2012).
A general reason for the model overestimating the albedo is that it does not take the albedo changes due to dust storms or volcanic dust deposition into account.For instance, the very low albedo values obtained at the T AC station (Fig. 4b) are due to tephra deposition on the glacier during the 2010 eruption of Eyjafjallajökull (e.g.Gudmundsson et al., 2012;Gascoin et al., 2017).Even though dust events do not cause as large changes in albedo as a volcanic eruption, they can still significantly lower the albedo (e.g.Painter et al., 2007;Wittmann et al., 2017) .As previously mentioned, the albedo in HIRHAM5 often reaches its yearly minimum value later in the summer than the observed.Such discrepancy could be explained by dust events, advancing or delaying the drop in surface albedo.Wittmann et al. (2017) investigated 10 dust events which occurred at the B ELA station in 2012, and found a lowering in the albedo during all events and showed that the dust storms have a significant effect on the resulting energy balance.
The error in the outgoing shortwave radiation is caused by errors in the albedo and the incoming SW.At the B AB station, the incoming radiation is slightly underestimated but the albedo is overestimated; hence, the outgoing SW is overestimated.The values at the four other stations are all underestimated, due to larger underestimations of the incoming SW radiation and lower albedo errors.As both the incoming and outgoing SW radiation are underestimated at most stations, the net SW shows a negative bias of ∼ −6 to −12 W m −2 at the AC and ELA stations, and of −22 and −28 W m −2 at the two AB stations.The resulting average model error at all five stations is −15.5 W m −2 .

Turbulent fluxes
As HIRHAM5 underestimates meteorological variables at all stations, similar underestimation is obtained for the turbulent fluxes (Table 3 and Fig. 5).The two AC stations have the largest differences and also the lowest correlation (0.45 and 0.49) between the AWS estimate and the HIRHAM simulation.The other three stations also have significantly lower values in the HIRHAM5 model than in the AWS model, but with higher correlation coefficients (0.69-0.73).
It is important to bear in mind that this comparison is a model-model comparison, so while the eddy flux model may give a good estimate of the turbulent fluxes, model errors still affect the results, e.g.due to the use of a constant roughness length.

Total energy balance
After the simulated components of the energy balance were evaluated against AWS observations, the total energy balance was estimated (see Table 3).The energy balance (E) is found using where LW net is the net LW radiation, SW net is the net SW radiation, and H s+l are the turbulent fluxes.Overall, the melt energy is underestimated, owing to all elements of the energy balance generally being underestimated.This is in large part due to the underestimation of the modelled incoming radiation.We attribute this to an error in the modelling of the clouds, but since both the incoming SW and LW radiation are underestimated, inaccurate cloud representation cannot be the the only source of the error.Errors in the interaction of clouds and radiation, e.g.error in the optical thickness of the clouds, or in the clear-sky fluxes, could partly explain these discrepancies.The underestimation of the incoming LW radiation could also be due to errors in the vertical atmospheric temperature gradient.Since the simulated outgoing LW radiation generally only has a small negative bias, the deviation in net LW radiation is governed by the incoming radiation.Errors in the simulated albedo mean that both the in-and outgoing SW radiation greatly contribute to the deviation in net SW radiation.These errors can be partly attributed to ash and dust deposition during volcanic eruptions and dust storms, which are not taken into account in HIRHAM5.In addition, errors in the simulated albedo also stem from snow cover that disappears too slowly compared to AWS records in the ablation zone.As a result, modelled albedo drops too slowly compared to the measured albedo.The underestimation of the net SW and LW radiation and the turbulent fluxes leads to underestimated melt energy.This contributes to the overestimation of the modelled snow thickness.
In order to estimate how much the different components contribute to the energy difference on a year-to-year basis, the mean difference between modelled and observed energy components during each summer (April-October) is shown for each station (Fig. 6).
At the B AC station, the contribution of the long-and shortwave radiation and turbulent fluxes to the energy difference is consistent for the entire period, with the error of each component being almost equal, varying between −25 and 0 W m −2 .At the T AC station, the error due to the three components is also of the same order of magnitude, except in 2010 and 2011 where the error in the net SW radiation is much larger than that in the other components.This is due to a large drop in the albedo as a result of the Eyjafjallajökull (2010) and Grímsvötn (2011) eruptions.The mean difference between observations and the simulations of the SW radiation for noneruption years is −3 W m −2 , whereas the radiation difference in 2010 is −106 W m −2 .Assuming the larger deviation from the mean in 2010 is only due to the volcanic eruption, the increase in available energy due to the eruption is 103 W m −2 .If it is further assumed that the surface was always at melting point, the increase in melt due to the 2010 Eyjafjallajökull eruption over the 128-day measuring period would be ∼ 3.1 m w.e. at this station.
At the ELA site, the contribution from the modelled turbulent fluxes to the energy balance deviation generally varies between ±10 W m −2 , except in 2013 where the bias is around −25 W m −2 .Modelled longwave radiation is consistently underestimated by −10 W m −2 .The deviation in the shortwave radiation is more variable, as expected from the results of the albedo comparison.Depending on whether bare ice was exposed or not, the albedo is generally either overor underestimated.For example, at B ELA , the ice surface was reached in, for example, 2007 and 2012, resulting in an overestimation of the albedo.In, for example, 2002 and 2009, however, the albedo was high the entire summer as no ice was exposed, resulting in an underestimation of the predicted albedo.
At the T AB station, both the net LW radiation and the turbulent fluxes agree well with observations for the entire period.The net SW radiation, however, is always underestimated, especially in the period 2001-2003 and 2011.These years, the measured albedo at the station goes below 0.1, while the HIRHAM5 albedo stays around 0.4.As previously discussed, this albedo bias, and hence underestimated net SW radiation, occurs because of an overestimation of the snow cover at the station due to an overestimation of the winter accumulation and possibly also the proximity of the equilibrium line.An underestimation of the incoming SW radiation, which we attribute to an error in cloud cover amount of clearsky fluxes, also contributes to this error.
At the B AB station, the longwave radiation bias is relatively constant with values close to 0 W m −2 for much of the measurement period.The absolute deviation due to the turbulent fluxes is less than 10 W m −2 for most of the period, although with slightly larger deviations from the period 2007-2010.The SW radiation is always underestimated at this station, mostly due to the previously discussed overestimation of the albedo.

At AWS sites
Scatter plots of measured and HIRHAM5 simulated SMB are shown in Fig. 7 and the average deviations are shown in Table 4.
The winter mass balance comparison allows to evaluate of the winter precipitation in HIRHAM5.The simulated mass balance at the B ELA and B AC are always underestimated, while the T AC stations is underestimated during all years but one (2012).The simulated value at the T AB station is overestimated over the whole period.The modelled mass balance at the B AB station has an almost equal amount of years which are over-and underestimated.Apparently the model either carries too much precipitation when the clouds reach the glacier, resulting in too much precipitation at the ice sheet margin, or more melting occurs at the ablation area stations during the winter months than the model estimates.
The summer SMB results are in good agreement with the results of the energy balance calculations.The summer SMB is generally overestimated, although it is underestimated occasionally at all stations except T AB .The ELA station has the largest amount of underestimated points, which is consistent with the findings from the energy balance calculations.Besides the errors introduced due to the underestimation of the energy balance, possible over-or underestimations of the modelled summer accumulation contribute to these errors as well.
Due to the difference in the summer and winter balance, the net balance at the B AC , T AC , and B ELA stations is generally underestimated in HIRHAM5, while the balance at the two AB stations is generally overestimated.This is due to a general overestimation of the winter balance in the ablation area, due to either an underestimation of the winter melt or an overestimation of precipitation, as discussed above.

At all measurement sites
SMB is also measured at 25-120 non-AWS sites, depending on the year (Fig. 1).In order to estimate how well the model represents the SMB at non-AWS sites, the data from all the sites between 1995 and 2014 were compared with the HIRHAM5 simulation (Fig. 8; Table 4).
The winter balance at all measured points is slightly overestimated by HIRHAM5 on average.However, this is mostly due to a large difference between measured and simulated SMB at the ice-covered, high-elevation, central volcano Öraefajökull (the white dots in Fig. 8).Only one site has been measured on this glacier for a few years only (Guðmundsson, 2000), in a spot that always receives a large amount of precipitation.However, since HIRHAM5 consistently overestimates the accumulation by 100-200 %, this one point has a large effect on the mean error.This is a well-known issue with hydrostatic models like HIRHAM5, as they char-acteristically overestimate the precipitation on the upslope and peaks in complex terrain.The reason for this is that the precipitation is calculated as a diagnostic variable -i.e. it is not governed by an equation that is a derivative of time, meaning that when the required conditions for precipitation are met in the local atmosphere, the precipitation appears instantaneously on the surface.Thus, the scheme does not allow horizontal advection of snow and rain by atmospheric winds, which is a key process in complex terrain, as it can force the precipitation downslope (e.g.Forbes et al., 2011).Without this effect, precipitation is generally overestimated at high peaks like Öraefajökull.Removing this location from the comparison, the total difference drops to one-third the difference with respect to the AWS sites only (−0.09 m w.e.).The reason the difference is smaller than for the AWS sites only is that more sites close to the edge of the ice cap are included.The winter balance at the measurement points at the outer parts of the icecap generally is overestimated in the model, and therefore these points partly offset the underestimation in the middle of the ice cap.
On average, the summer ablation is underestimated, which is consistent with the findings from the AWS stations that there is an average underestimation of the energy available for melt.The mean error and RMSE is only slightly larger than at the AWS sites.
The mean net balance is overestimated by approximately the same amount as the summer balance, partly due to the low mean deviation in the winter SMB.Due to the large deviation at Öraefajökull in the winter SMB, the Öraefajökull points clearly have the largest bias.If these points are excluded, a RMSE closer to that for the AWS locations is found (1.1 m).

Reconstructing the SMB of Vatnajökull
Spatial maps of the (uncorrected) average winter, summer, and net SMB from the 1980-1981 glaciological year until 2013-2014 are shown in Fig. 9.The approximate location of the average ELA is marked in the figure.The model captures the position of the ELA fairly well, but at, for example, Brúarjökull, where the average ELA is at 1200 m, the position of the average ELA is at a too high elevation.The average deviation between observation and model over the observation period at each measurement location is also shown in Fig. 9, in order to give an indication of the average error of the model at different parts of the ice cap.The winter balance (Fig. 9e) is generally overestimated at low elevations and underestimated at high elevations, except for at Öraefajökull, where there is a large overestimation of the winter balance, as discussed in the previous section.As can be seen in Fig. 9e, there is generally a low SMB bias at high elevations and a high SMB bias at low elevations during the summer.This is consistent with the comparisons with AWS stations, as we found that the bias in the energy available for melt was smaller at high elevation than at low elevation (see   This was partly due to a larger albedo bias for stations in the ablation zone than for stations in the accumulation zone. In addition to the spatial maps, the winter, summer, and net mass balances of Vatnajökull were calculated for the entire simulation period, and the results were compared with an estimate of the specific balance from 1995 to 2014, created by interpolation of the mass balance measurements (e.g.Pálsson et al., 2015); see Fig. 10.The model prediction of the mean specific summer mass balance generally fits well with the interpolated observations, with an overall difference of only 0.06 m w.e.The largest deviations are obtained in 1995, where ablation is overestimated in the simulation, and in 1997, 2005, and 2010-2012, where ablation is underestimated, most likely due to ash depositions on the glacier following the 1996 Gjálp eruption, the 2004 and 2011 Grímsvötn eruptions or the 2010 Eyjafjallajökull eruption, which are not taken into account in the model.
Excluding the years where the albedo was affected by volcanic eruptions, the average difference becomes smaller but the model also predicts slightly too much ablation, as the difference becomes −0.02 m w.e.
There is a shift in the summer and annual mass balance calculated by the model and the in situ MB measurements around 1996, with a generally more negative mass balance after 1996 than before.This is consistent with the increase of the annual mean temperature of Iceland in the mid-1990s, which resulted in a mean annual temperature ∼ 1 K higher in the decade after than the decade prior to 1995.This is likely linked with atmospheric and ocean circulation changes around Iceland, as there was a rapid increase in ocean temperatures off the southern coast in 1996 (Björnsson et al., 2013).
The specific winter mass balance is overestimated in HIRHAM5 for the entire measurement period with an average of 0.54 m w.e.Due to this difference, and only the small negative mean difference in summer mass balance, the annual mass balance of Vatnajökull is overestimated every year with an average difference of 0.50 m w.e.
However, this is mostly due to the large overestimation of the winter accumulation on Öraefajökull; comparison with the mass balance measurements showed that the model overestimated the winter accumulation by 100-200 % compared with the observations.In an attempt to estimate how much this error affects the results, a simple correction was added to the Öraefajökull points by reducing the simulated winter SMB by 50 %.The correction was added to four model grid points around Öraefajökull, due to the high (> 10 m yr −1 ) annual specific mass balance in these points (see Fig. 9a).The resulting modelled winter and annual specific balance are shown in Fig. 11.The winter balance is still overestimated, but the difference between modelled and interpolated values has been reduced to only 0.1 m w.e.In addition, the average difference between the HIRHAM5 and interpolated annual SMB drops to only 0.08 m w.e.

Comparison with constant ice albedo simulation
In order to quantify the changes in the model performance resulting from the new albedo scheme used in this study, which utilizes an albedo map based on MODIS data (Gascoin et al., 2017), the results are compared to those of a previous run using a constant ice albedo of 0.3.The average difference in albedo and mass balance over the period 2001-2014 in each grid point are shown in Fig. 12, as well as the position of the AWS stations.
There is little to no difference between the two runs in the accumulation zone, due to the year-round snow cover.In the ablation zone, however, using the MODIS ice albedo map has a large effect on the simulated albedo.The largest differences are found on the southern outlet glacier Skeiðarárjökull, which is unfortunately a glacier where no mass balance or AWS measurements have been conducted.The B AB and B ELA stations are located in areas that are affected by the ice albedo, either because ice is exposed (B AB ) or because the underlying surface contributes to the albedo (B ELA ).The T AB station is located in the ablation area, but the ice surface is never exposed in the model due to an overestimation of the winter accumulation.The albedo estimate at this station was therefore not improved by using the MODIS albedo.
When the model is run with the constant ice albedo of 0.3, the amount of ablation will be lower and thus the specific summer balance will be higher.Compared to the simulation using the MODIS map (Fig. 11), the constant ice albedo simulation results in an increase in the specific summer SMB by an average of 0.37 m w.e., or 18 %, per year for the period 1995-2014.The increase in the summer SMB ranges from 14 cm (in 2014) to 85 cm (in 2001) and the percentage increase varies between 8 % (in 2011) and 39 % (in 1995).As the winter balance is not dependant on the ice albedo, there are no changes in the specific winter SMB between the two simulations.

Conclusions
The comparison of a HIRHAM5 simulation with data from five AWSs on Vatnajökull ice cap allows us to evaluate the model performance.By comparing observations from April to October with model output, it was found that the model simulates the surface energy balance components and surface mass balance well, albeit with general underestimations of the energy balance components.Even though the energy balance was generally underestimated, the model simulated the near-surface temperature well.The reason for this is that the comparisons only use observations from the summer months, where the glacier surface is generally at the melting point, and thus the energy is used for melting and not for raising the temperature of the surface.
The modelled incoming radiation is underestimated on average in both the shortwave and longwave spectrum, which we suggest is due to biases in the modelling of the cloud cover combined with errors in the optical thickness in the short-or longwave spectrum, or errors in the clear-sky fluxes.
Whereas the modelled outgoing LW radiation component is within the uncertainty of the LW observations at the five stations, which is consistent with the ability of the model to capture surface temperatures, there was a larger difference between the modelled and measured outgoing SW radiation.This is partly due to the underestimation of the incoming SW radiation and partly due to inaccuracies in the simulated albedo.The albedo was simulated using an iterative, temperature-based albedo scheme (Nielsen-Englyst, 2015) with a bare ice albedo determined from MODIS data (Gascoin et al., 2017).The simulated albedo was generally overestimated during the summer and did not reach the lowest yearly value as early in the year as the measured albedo, particularly in the ablation zone.This was attributed to an overestimation of the snow cover in the ablation zone, an overestimation in the MODIS ice albedo compared with AWS observations, and the fact that the model does not account for the effect of volcanic dust deposition during eruptions and dust events on the albedo.A possible means of capturing dust storms or eruptions into the model is to implement a stochas- -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7  overestimated by 0.06 m w.e. on average -i.e.there is generally too little ablation in the summer, with too much ablation in 1995, and too little ablation in years with, or following, volcanic eruptions.The winter balance is overestimated by 0.5 m w.e., mostly due to a large overestimation at the high elevation glacier Öraefajökull.This overestimation of accumulation at high elevation is characteristic for hydrostatic RCMs (Forbes et al., 2011).If the overestimation at these points is corrected, we estimate that the simulated winter balance would fit well with the observations, as the overestimation of the balance would drop to around 0.1 m w.e.That the model catches the changes in the specific mass balance well over the mass balance measurement period, and also captures the shift in mass balance in the mid-1990s, gives us confidence that the model estimates the specific mass balance of Vatnajökull well over the entire simulated period from 1980 to 2014.HIRHAM5 is therefore a useful tool to expand the time series of the specific SMB beyond the measurement years.However, as ERA-Interim reanalysis data only go back to 1979, the model would need to be forced at the lateral, for example, by output of a general circulation model.However, using other reanalysis data probably leads to different errors; this needs further investigation.The model could also be a useful tool to estimate the future evolution of the SMB of the ice cap, but this would also require a different forcing at the lateral boundary like general circulation model output.This would most likely introduce larger biases than the ones found using ERA-Interim, and the magnitude of these biases would need to be estimated and corrected before using the model for future projections.
Data availability.HIRHAM5 output is freely accessible from http: //prudence.dmi.dk/data/temp/RUM/HIRHAM/GL2/, as are MODIS data from https://modis.gsfc.nasa.gov/data/.Measurements from automatic weather stations and from in situ mass balance surveys are partially owned by the National Power Company of Iceland and are therefore not publicly available at this time.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.(a) The average location of the AWS sites.Only the labelled sites were used in this study.(b) The average location of the mass balance sites from 1995 to 2014.The coloured lines connect mass balance sites along a transect.Not all mass balance sites were measured every year.
AWS data from the period 2001-2014 for three Brúarjökull stations and two Tungnaárjökull stations are considered, as well as SMB point measurements from 1995 to 2014.All stations were operated during the summer months, but since 2006 the lowest Brúarjökull station has been operated year round.Comparisons are made between daily averages from the HIRHAM5 model and the in situ observations collected at the AWSs.HIRHAM5 daily means are calculated from 6hourly outputs, while the AWS daily means are calculated from observations at 10 min intervals.L. S. Schmidt et al.: Evaluating the surface energy balance in HIRHAM5 over Vatnajökull, Iceland

Figure 2 .
Figure 2. Scatter plots of the measured (a) surface pressure, (b), air temperature at 2 m, (c) relative humidity at 2 m, and (d) wind speed at 2 m, by stations on Bruarjökull (red) and Tungnaárjökull (blue) versus the same components simulated by HIRHAM5 at the same locations.

Figure 3 .
Figure 3. Scatter plots of the measured longwave radiation components, LW↓ and LW↑, by stations on Brúarjökull (red) and Tungnaárjökull (blue) versus the LW radiation components simulated by HIRHAM5 at the same locations.The dashed line corresponds to ±10 %, i.e. the manufacturer-reported uncertainty of the AWS measurements.

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Scatter plots of the measured shortwave radiation components, (a) SW↓, (b) albedo, and (c) SW↑, by stations on Bruarjökull (red) and Tungnaárjökull (blue) versus the shortwave radiation components simulated by HIRHAM5 at the same locations.The dashed line corresponds to the uncertainty of the measured AWS components.

Figure 7 .
Figure 7.Comparison of the winter, summer, and net mass balance from the period 1995-2014 between the mass balance measurements at the five AWS sites and the HIRHAM5 simulation.

Figure 8 .
Figure 8.Comparison of SMB measurements from Vatnajökull ice cap from 1995 to 2014 and HIRHAM5 simulated values.Different colours represent different outlet glaciers; see Fig. 1b.The white dots are from a point on Öraefajökull.

Figure 9 .
Figure 9.The average (a) winter, (b) summer, and (c) net SMB simulated by HIRHAM5 from the 1980-1981 glaciological year to 2013-2014.The contour lines marks the approximate location of the ELA, which generally lies between approximately 1100 and 1300 m elevation.Panels (d)-(f) show the average deviation between model and observations over the observation period (1992-2014) for each measurement location for the (d) winter, (e) summer, and (f) whole glaciological year.

Figure 10 .Figure 11 .
Figure 10.Average summer (red lines), winter (blue lines), and net (green lines) specific surface mass balance for the whole of Vatnajökull.The solid lines are the mass balance of Vatnajökull based on mass balance measurements and manual interpolation, while the dashed lines are the mass balance as simulated by HIRHAM5.

Figure 12 .
Figure 12.Difference in (a) mean albedo and (b) mean SMB in m w.e. for the period 2001-2014 between two runs with HIRHAM5, one using a MODIS bare ice albedo map and the other with a constant ice albedo.The locations of the AWS stations used in this study are shown with black circles.

Table 1 .
Average measured elevation and average bias of the interpolated HIRHAM5 elevation at each station for 2001-2014.

Table 2 .
Comparison of the surface pressure p sl , air temperature at 2 m, T 2 m , relative humidity r 2 m , and wind speed u, from HIRHAM5 simulations and AWS measurements during the summer months (April-October) for the period 2001-2014.The HIRHAM5 bias (HIRHAM5-AWS), the root-mean-square error (RMSE), the percentage error, and the correlation (r) are shown.

Table 3 .
Comparison of incoming and outgoing long-and shortwave radiation, albedo (α), turbulent fluxes (H s+l ), and total energy (E) from HIRHAM5 simulations and AWS measurements during summer months (April-October) from the period 2001-2014.The HIRHAM5 bias (HIRHAM5-AWS), the root-mean-square error (RMSE), the percentage error, and the correlation (r) are shown.

Table 4 .
Comparison of HIRHAM5 and mass balance measurements, both at AWS sites and for all measuring sites on Vatnajökull.