Sensitivity of Glacier Runoff to Winter Snow Thickness Investigated for Vatnajökull Ice Cap , Iceland , Using Numerical Models and Observations

Several simulations of the surface climate and energy balance of Vatnajökull ice cap, Iceland, are used to estimate the glacier runoff for the period 1980–2015 and the sensitivity of runoff to the spring conditions (e.g., snow thickness). The simulations are calculated using the snow pack scheme from the regional climate model HIRHAM5, forced with incoming mass and energy fluxes from the numerical weather prediction model HARMONIE-AROME. The modeled runoff is compared to available observations from two outlet glaciers to assess the quality of the simulations. To test the sensitivity of the runoff to spring conditions, simulations are repeated for the spring conditions of each of the years 1980–2015, followed by the weather of all summers in the same period. We find that for the whole ice cap, the variability in runoff as a function of varying spring conditions was on average 31% of the variability due to changing summer weather. However, some outlet glaciers are very sensitive to the amount of snow in the spring, as e.g., the variation in runoff from Brúarjökull due to changing spring conditions was on average 50% of the variability due to varying summer weather.


Introduction
Glacier mass loss is a key contributor to sea-level change (e.g., [1]) as well as changes in river volume and seasonality (e.g., [2,3]).Worldwide, melt water from glaciers e.g., affect sediment transport, the availability of freshwater, and hydropower production (e.g., [4]).As the climate warms, the amount of runoff from glaciated areas will initially increase until a peak runoff is reached [5], following which the amount of runoff will decrease (e.g., [3,5,6]).On seasonal timescales, glacier mass loss affects the river flow by storing winter precipitation and releasing the most water during the summer months.A warming climate will change the timing of this seasonality, e.g., due to an earlier onset of glacier melt (e.g., [2]).
Like many other glaciers and ice caps worldwide (e.g., [7]), Icelandic glaciers are currently losing mass and retreating because of a warming climate [8].Icelandic glaciers currently store the water equivalent of ∼20 years of precipitation [9], which is released over a range of timescales [2].Due to the large discharge from Icelandic glaciers, glacial runoff is an important component of Icelandic hydrology, and many of the main Icelandic rivers have a significant glacial component.Since hydropower stations produce over 80% of the electricity in Iceland, knowledge of the amount of runoff from glaciers is not only important due to the effect on sea level, but also for harnessing of hydropower.
The surface energy balance of a glacier, and therefore the runoff, has a large sensitivity to changes in albedo.Icelandic glaciers have a large spatial and temporal variation in albedo, spanning from less than 0.1 for dirty ice in the ablation zone to 0.9-0.95 for new snow.The low ice albedo is due to layers of volcanic ash, which dominate the spectral properties of the ice (e.g., [10,11]).How much glacier ice is exposed during the summer and the timing of when the winter snow melts away is therefore expected to have a large effect on the runoff over the ablation season.Due to climate change, the amount of solid precipitation is expected to change (e.g., [6,7,12]) and it is, therefore, important to investigate the sensitivity of ablation to snow thickness changes.
The ablation of glaciers can be estimated using either physical or empirical methods (e.g., [13]).While physical models, i.e., energy balance models, can provide direct estimates of the energy balance, a lack of observations and the many input parameters can make them impractical.Empirical models, i.e., degree-day models, describe the statistical relations between melt and temperature or other weather parameters.These models can be more practical but need to be tuned for every location and provide less information about the surface climate.Neither of these models can estimate the accumulation, and therefore a physically based atmospheric model is needed to simulate the mass balance.
Global Climate Models (GCMs) and Regional Climate Models (RCMs) combine atmospheric and surface physics.They are characterized by the explicit solution of the hydrodynamic equations of the atmosphere, which is of vital importance to the mass balance of glaciers as it e.g., determines the transfer of moisture from the source regions (e.g., [13]).RCMs are forced by a GCM or reanalysis at the boundaries of their domain, and by downscaling the GCM they provide added value to the simulations of the cryosphere.Since they are running on a smaller domain, they can run at higher resolution and include more complex models of the cryosphere.Numerical high-resolution RCMs, such as RACMO [14], WRF [15], MAR [16], and HIRHAM5 [17], are therefore valuable tools to estimate mass balance at the surface of glaciers, and much work has gone into evaluating these models over Greenland (e.g., [18,19]) and Antarctica (e.g., [20,21]).To advance and evaluate models used for regional downscaling through international cooperation, the Coordinated Regional Downscaling Experiment (CORDEX) was started in 2009 [22].Through this project, a collection of historical and future simulations from many different RCMs were made publicly available.
Another type of models used for regional downscaling is Numerical Weather Prediction (NWP) models.These models use physical models of the atmosphere and oceans to predict the weather based on current weather conditions.RCMs are generally derived by NWPs, but the two categories of models serve different functions: NWPs predict the weather while RCMS project climate scenarios.The major difference between the two types of models is therefore the scale of the simulations.NWPs generally operate over smaller domains and much smaller time scales (days as opposed to decades for RCMs).For this reason, NWPs are often run at a higher resolution than RCMs, even if currently a big challenge is the execution of climate simulations at 1 km resolution (convection-permitting simulations), for example in the frame of CORDEX-CORE.
However, projects to use NWPs as a tool for climate reanalysis have been established.Over the Iceland domain, the NWP HARMONIE-AROME has been used for the ICRA reanalysis project [23] to estimate the climate from 1980 to 2016.The project has e.g., been focusing on adding an improved surface scheme for Icelandic glaciers, to achieve more accurate mass balance estimates.Overall, the reanalysis showed a good correlation with available Automatic Weather Station (AWS) observations over Iceland [23].
Although RCMs and NWPs are important tools for climate and weather predictions, they often have systematic biases in temperature and/or precipitation (e.g., [24,25]), which need to be corrected to get the most accurate estimates of changes in mass balance and runoff.Since runoff is a nonlinear function of precipitation, model bias correction is especially important for hydrological modeling (e.g., [25,26]).Numerous bias correction methods have been developed to remove these model errors (e.g., [27,28]), for example by adjusting the modeled mean of the corrected variables to match observations more closely.Often a time-invariant method is used, where the error correction found over the evaluation period is assumed to be applicable to any time period.However, this method works best for short-term modeling, as the bias correction cannot be expected to be constant over longer timescales (e.g., [28]).
In this study, the RCM HIRHAM5 [17] and the NWP HARMONIE-AROME [29] are used to simulate the runoff from Vatnajökull and to quantify the impact of the snow thickness in spring on the runoff the following summer.The precipitation biases in both HIRHAM5 and HARMONIE-AROME are evaluated, and a simple bias correction is applied to the model precipitation.The impact of the spring conditions is then tested by setting up a model experiment which uses the simulated spring conditions at the end of the 31st of March for all years in the period 1980-2015 as a starting point for model runs forced by summer weather of all years in the same time period.Distributions of the likely runoff from Vatnajökull are created based on those 36 years, and we investigate to what extent spring snow thickness affects summer melt.While several studies have investigated the effect on increasing accumulation on the runoff (e.g., [30][31][32]), this is the first attempt to quantify the effect of spring snow cover on glacier melt.

Study Site
Vatnajökull ice cap is the largest temperate ice cap in Europe, with an area of 8100 km 2 .It is located close to the southeastern coast of Iceland, and covers active volcanic areas where eruptions are frequent (e.g., [33]).The surface elevation spans from 0 to 2110 m above sea level (a.s.l.).Most of its precipitation falls as snow, as the average temperatures on the ice cap is close to or below freezing.The summer balance is normally negative for the entire ice cap, but the central portions might turn slightly positive during repeated cold spells.At the highest areas of the glacier, days with melting only occur 10-20 days of the year, while at lower levels the ablation season generally lasts 3-4 months [9].Vatnajökull is affected by frequent dust storms and tephra from volcanic eruptions.Tephra layers in the glaciers ice dominate its spectral properties [10], but the effect of tephra on the accumulation zone is generally small.In addition to surface melting, geothermal activity at the bed provides an on average small contribution to the total melt [9].
Surface and bedrock maps of Vatnajökull have been produced using interpolated profiles of radio-echo sounding and precision altimetry [9].These have been used in conjunction to delineate the locations of glacial water catchments which feed glacial rivers.

Energy and Mass Balance
AWSs have been operated on Vatnajökull since 1994, with 1-13 stations measuring on the ice cap during the summer months (e.g., [34,35]).The 2 m temperature, 2 m relative humidity, 4 m wind speed, and the 4 m wind direction have been measured for the entire period, while the radiation components have been measured since 1996.Observations from five AWSs, three on Brúarjökull and two on Tungnaárjökull that have been operated at approximately the same location since 2001 are used to evaluate the simulated incoming radiative fluxes in HIRHAM5 and HARMONIE.In addition, estimations of the albedo at the stations in the ablation zone are used to evaluate how well the model captures the timing of the albedo drop from snow to ice values.The average locations of the stations are shown in Figure 1.The observations are conducted at 10 min intervals, but daily averages are used for the comparisons.To evaluate the modeled snow thickness in spring, the modeled winter balance was compared to in situ winter mass balance measurement at several sites on Vatnajökull.In situ mass balance measurements have been carried out twice a year every glaciological year since 1991-1992, with an average of 60 measured locations every year.The uncertainty in the measurements has been estimated to be ±0.3 m water equivalent (w.eq.) [8,36].There are no precipitation observations for Vatnajökull, but it is generally assumed that the winter mass balance is equal to the winter precipitation (e.g., [10]).Almost all winter precipitation on Vatnajökull is snowfall and there is only little melt in the summer besides in the lowest part of the ablation zone.The winter balance is therefore considered a good measure of the winter precipitation and will therefore in this study be used to evaluate the simulated precipitation in the two models.

River Runoff
Substantial hydrological measurements have been conducted in glacial rivers in Iceland mainly due to their importance for hydropower production.Therefore, time series spanning over decades are available for several glacial rivers.In this study, the runoff into two rivers, Jökulsá á Dal and Skaftá, is used to evaluate the model runoff; both have a large glacial component.The outlines of the glacial water catchments (GWCs) that contribute to these rivers are shown in Figure 1.
About 80% of the runoff from Brúarjökull, a northern Vatnajökull outlet, flows into Hálslón reservoir, which lies on the uppermost former path of the Jökulsá á Dal river.Hálslón reservoir water level has been monitored since the impoundment of the dam in autumn 2007 [37].The runoff to the reservoir is estimated based on a relationship between reservoir water level and storage in a time series of daily averaged inflow.Before Hálslón reservoir was created, the runoff from Brúarjökull was e.g., measured further down the Jökulsá á Dal river at Hjarðarhagi (VHM110) [37].To expand the time series for Brúarjökull, the Hjarðarhagi observations are used from 1980 to 2007, after which the runoff was routed into Hálslón reservoir from 2007 to present.The size of the Hjarðarhagi catchment area is ∼3300 km 2 [38] and the size of the Hálslón catchment area is about 1800 km 2 , where approximately 1200 km 2 is covered with Brúarjökull.The glacial component of the catchment spans from ∼750 to 1600 m a.s.l.Runoff from non-glaciated areas has not been estimated for the Jökulsá á Dal sites, which adds further uncertainty to the results.About 80% of the catchment area of Hálslón is covered by Brúarjökull while the rest is non-glaciated land with a seasonal snow pack.It can very tricky to distinguish between glacier melt and seasonal snow pack melt from only an inflow time series.However, the range of snow pack storage has been observed since 1956, which can give a rough estimation of the contribution.Since 1956, the average seasonal spring snow pack has been estimated to store ∼0.19 ± 0.08 km 3 of water [37].Since 2007, the yearly runoff into Hálslón has been observed to be between 3.0 and 4.5 km 3 , and we therefore expect the runoff from non-glaciated areas in spring to contribute less than 10% to the yearly observed runoff.However, since the snow pack generally melts between April and middle of July, it can have a large contribution during the spring months.In 2013 there was an unusually high amount of seasonal snow in the catchment corresponding to approximately 0.25 km 3 of stored water.In the runoff time series, 0.58 km 3 water was observed between April and June, meaning the off glacier seasonal snow contributed about 40% to the spring runoff in 2013.
In addition, fall precipitation can also have a large contribution on the runoff time series.This was e.g., observed in autumns of 2015 and 2016.
Skaftá river receives runoff from the Skaftá GWC, which includes the outlet glacier Skaftárjökull and part of the outlet Tungnaárjökull.The closest hydrometric station to the glacier is at Sveinstindur (V299), 25 km down-river from the glacier margin, and has been operated since 1986 [39].The Sveinstindur station has a catchment area of 714 km 2 [40], where approximately 400 km 2 is glacierized.The glacial area spans from approximately 700-1600 m a.s.l.The runoff data comes with a quality flag for each daily datapoint.Only data that are of full quality have been used for this evaluation, while data that have needed corrections due to ice interference on the runoff or other problems have been discarded.To isolate the glacial runoff in the river runoff record, the runoff originating from ice-free areas is estimated using the hydrological model WaSiM [41].WaSiM is a physically based distributed model, which offer methods of varying complexity to simulate different elements of the hydrological cycle depending on available observations.It has in recent years been used for hydrological modeling of Iceland and includes interaction between glaciers and ice-free areas such as e.g., interflow.A detailed description of the setup for Iceland is given in Jónsdóttir [42] and Einarsson [43].The WaSiM model was run with four equally good parameterizations by the Icelandic Meteorological Office [39] and the four outputs were used to assess the uncertainty of the model estimation.The difference between the four outputs is small overall, as the deviation from the mean of the four model outputs is on average ±0.6%.However, it can deviate by up to ±44.0% on specific days.For the comparison presented here, a mean of all four WaSiM runs was used as an estimate for the non-glacial runoff.
There is a substantial uncertainty in the WaSiM simulations, which leads to a negative runoff from the glacier during years with a low runoff during early spring or late autumn.Discarding these negative values means that the total accumulated runoff might be overestimated.However, the days with negative values are removed in both simulated and observed time series, and therefore this should not be an issue for the comparison.
Jökulhlaups from the Skaftár cauldrons (e.g., [44,45]), two subglacial lakes at geothermal areas that collect water and release in jökulhlaup events, provide another complication in the runoff time series.The two cauldrons are located just north of the Skaftá GWC and drain every 2-3 years [44].Since jökulhlaups are not simulated in the model, runoff on days when jökulhlaups are in progress are removed from the time series.The cauldrons do not contribute to the runoff into Skaftá when no jökullhlaup event is in progress, as the water gathers in the subglacial lakes.The cauldrons are in a different water catchment, and not a part of the Skaftá GWC, and therefore jökullhlaups do not affect the water balance of this catchment except on event days.

HIRHAM5 and HARMONIE-AROME
In this study, we apply the RCM HIRHAM5 [17] and the NWP HARMONIE-AROME [29].HIRHAM5 is a hydrostatic RCM which combines the dynamical core of the HIRLAM7 numerical forecasting model [46] and physics schemes from the ECHAM5 general circulation model [47].A detailed description of the model configuration can be found in Christensen et al. [17].It is run over a domain containing Greenland and Iceland, with a horizontal resolution of 0.05 • × 0.05 • on a rotated-pole grid.This corresponds to ∼5.5 km resolution.The atmosphere has 31 vertical levels, from the surface to 10 hPa.The time step is set equal to 90 s.The model is forced by ERA-Interim reanalysis data [48], which is a global atmospheric reanalysis by ECMWF spanning back to 1979, at the lateral boundaries at 6 h intervals.ERA-Interim sea surface temperatures and sea ice fraction are imposed in ocean points.The model output is at 6 h intervals.HIRHAM5 model simulations have been successfully validated over Greenland (e.g., [18,19,[49][50][51][52]) and Iceland [10] using AWS and ice core data.While evaluating the incoming energy and mass balance fluxes over Vatnajökull in HIRHAM5, Schmidt et al. [10] found a generally good agreement between observations and simulations in the incoming energy fluxes.However, one important drawback with the simulation from HIRHAM5 was an overestimation of the winter accumulation in the ablation zone of several outlet glaciers, e.g., Brúarjökull, Dyngjujökull and Tungnaárjökull, as well as in areas with high orographic forcing.This was partly attributed to the model using hydrostatic physics and not allowing horizontal advection of snow.
HARMONIE-AROME is a non-hydrostatic convection-permitting NWP model and is developed in a cooperation between HIRLAM, ALADIN and Meteo-France [29].The forecast model was based on the AROME-France model (e.g., [53]), but now differs from the original model in various aspects [29].Details on the model configuration are described in Bengtsson et al. [29].In autumn 2015, the Icelandic Meteorological Office started a reanalysis project for Iceland (ICRA) using the HARMONIE-AROME model [23].Details on the reanalysis setup over Iceland can be found in Nawri et al. [23].ICRA currently spans from 1 September 1979, until 31 December 2016.It is run over a domain containing only Iceland at a horizontal resolution of 0.025 • × 0.025 • , corresponding to ∼2.5 km.The atmosphere has 65 vertical levels, from the surface to 10 hPa.The time step is set equal to 60 s.Like HIRHAM5, HARMONIE-AROME is forced by ERA-Interim reanalysis data at the lateral boundaries at 6 h intervals.The output is given at 1 h intervals.Since HARMONIE-AROME is non-hydrostatic and calculates precipitation diagnostically, it provides a better representation of the accumulation in areas with high orographic forcing than HIRHAM5 (more details in Section 4.1).

Snow Pack Scheme
While the original HIRHAM5 [17] used unchanged ECHAM physics, an updated model version is used in this study which expands the 5-layer surface scheme in ECHAM to 25 layers.In addition, it includes a dynamic surface scheme that explicitly calculates the surface mass budget, takes melting of snow and bare ice into account, and resolves the retention and refreezing of liquid water in the snow pack [19,51].
This updated surface scheme in HIRHAM5 has previously been adopted to and evaluated for Vatnajökull by Schmidt et al. [10].An updated snow-albedo scheme was tested in this study, depending both on the age of the snow and the surface temperature, and the snow aging parameters were tuned to better fit with observations from AWSs operated on Vatnajökull.The ice albedo, emerging as snow melts away in the ablation zone, was improved by using a map of the ice albedo based on MODIS observations [10,11].HIRHAM5 does not have a scheme for routing of runoff.The only delay taken into account is a runoff timescale based on the surface slope of the glacier.
The snow pack scheme was created for HIRHAM5, and any new implementations to the scheme are first implemented into the offline version of the model [19].The offline version of HIRHAM5 only contains the HIRHAM5 glacier model, which is run offline forced by 6-h surface energy (incoming shortwave (SW↓) and longwave (LW↓) radiation and turbulent fluxes) and mass fluxes (snow, rain, evaporation, and sublimation) from a previous HIRHAM5 experiment.While a fully coupled, high-resolution HIRHAM5 run is computationally very expensive, this offline model offers a quicker alternative to test new model implementations.In this study, we are using this offline model, which we refer to as the HIRHAM5 snow pack scheme.This model can be forced by incoming energy and mass components from another climate model, although it will no longer be a fully coupled system and missing feedback must be taken into considerations.
The ICRA HARMONIE-AROME runs have a less advanced snow pack scheme, which e.g., does not track the thickness of the snow layers, making the HARMONIE-AROME snow pack scheme less suited for this study.In addition, since the HIRHAM5 snow pack model is computationally cheaper than running the online HARMONIE-AROME model, it is more applicable for the many model runs required in this study.

Experiment Design
To use the validated offline snow pack model from HIRHAM5 but also make use of the improved incoming mass flux and higher resolution of the HARMONIE-AROME meteorological forcing, a combination of the two models is used in this study.This is implemented by running the snow pack scheme from HIRHAM5, in the same way as described in Schmidt et al. [10], but forcing the snow pack model with the meteorological forcing from HARMONIE-AROME.The snow pack scheme is forced every 6 h by surface energy (incoming radiation and turbulent fluxes) and mass fluxes (snow, rain, and evaporation) from HARMONIE-AROME.To initialize the snow pack model, a 100-year spin-up was performed by repeating the forcing from 1980.At the end of the spin-up, all variables had converged in all model layers.The final state of the spin-up was then used for a 1980-2015 simulation, forced by 6 h HARMONIE-AROME input.A schematic diagram showing the coupling is shown in Figure 2. To estimate the effect of the spring conditions on the summer runoff, sensitivity runs are conducted by repeating the modeled spring conditions of each year from 1980 to 2015, followed by the weather of all summers in the same period.Firstly, the HIRHAM5 snow pack scheme is run from 1980 to 2015 over all glacier points using HARMONIE-AROME incoming mass and energy fluxes, saving the state of the snow pack (snow thickness, temperature, albedo etc.) every year at the end of 31 March.Snow ablation generally starts between late April and May, so this date was chosen to include the entire ablation season.The state of the surface on 31 March of each year is then used to initialize a model run using the atmospheric forcing (i.e., incoming energy and mass fluxes) for each of the 36 model years starting on 1 April.By the end of the model run, 36 × 36 = 1296 summers have been simulated.In addition, to create a case with only snow melt, the model is also initialized with a thick enough snow layer on the surface, so the underlying ice is never exposed.This is done to estimate how much the melt is affected by the change in albedo when the ice is exposed when compared to only snow melt.

Evaluation of Model Precipitation
As mentioned above, the energy and mass balance in HIRHAM5 has previously been evaluated over Vatnajökull by Schmidt et al. [10], and an ice albedo map using MODIS data [54] was implemented in the model which improved the albedo simulations.Offline model simulations with HIRHAM5 showed that it overestimates the winter balance in the ablation zone as well as in areas with high orographic forcing.HARMONIE, on the other hand, appears to simulate the distribution of precipitation over the glacier more satisfactorily, overall.Scatter plots of the measured winter balance and model results, which have been interpolated onto the measurement locations using bilinear interpolation, are shown in Figure 3a,b and Table 1.The HARMONIE results are clearly closer to the one-to-one line than HIRHAM5, which is also evident from the lower root mean square error (RMSE) (0.6 m w.eq.for HARMONIE vs. 1.2 m w.eq.for HIRHAM5).The larger outliers, which correspond to the measuring site on Öraefajökull, which is the highest area within Vatnajökull, are also closer to the observations in HARMONIE than in HIRHAM5.Only a few observations are available for this location, as the summit has only been surveyed in the period 1993-1998 [55].However, comparisons with geodetic mass balance from 1982 to 2010 [56] show that the simulated mass balance overestimation in both models occurs for the entire investigated period.
To assess the spatial distribution of the errors, mean difference maps were created by comparing the model simulations to interpolated mass balance maps based on measurements (Figure 3c,d).The distribution of the differences is similar in both models, with the steep, high-elevation areas in the south and south-east generally being too wet in the models and the middle of the ice cap generally being too dry.However, in areas with high orographic forcing, HARMONIE simulates the observations better than HIRHAM5.The difference on Öraefajökull is smaller in HARMONIE than in HIRHAM5, and the high-elevation area around Bárðarbunga is significantly better represented in HARMONIE.The precipitation in the center of the ice cap is also less underestimated in HARMONIE.
The other fluxes needed to force the snow pack model, i.e., the incoming radiative and turbulent fluxes, have also been evaluated in both models using observations from five AWSs from 2001 to 2014.The mean differences are shown in Table 1.It is important to note that the turbulent fluxes are not compared to observations but to a one-level eddy flux model [57] which uses the AWS observations of temperature, humidity, and wind speed at 2 m to estimate the turbulent fluxes.The incoming solar radiation is better represented in HARMONIE, with both a smaller model difference and RMSE, as well as a higher correlation.The average difference is slightly smaller in HIRHAM5 for both the incoming longwave radiation and the turbulent fluxes, but the RMSE and correlation are better for both in HARMONIE.Overall, the incoming radiative flux seems to be better simulated by HARMONIE than HIRHAM5, at least at the five AWS locations.
Based on these comparisons, HARMONIE simulates the incoming mass and energy fluxes for the current Icelandic climate conditions better than HIRHAM5.The reason is the more complex model physics in HARMONIE than in HIRHAM5.Due to the higher-resolution and non-hydrostatic physics used in HARMONIE, the model can resolve non-hydrostatic features in the atmosphere not captured by the lower-resolution, hydrostatic HIRHAM5 model.The higher resolution also means that the topography is better resolves, which is especially important for simulations over the steep slopes of Öraefajökull.Another important distinction between the models is the treatment of precipitation.In HIRHAM5, precipitation is treated as a diagnostic variable.When the right conditions for precipitation are met, it immediately appears on the surface, and thus no horizontal advection of rain and snow occurs.In HARMONIE-AROME, the precipitation is treated as a prognostic variable and it therefore allows horizontal advection of snow and rain.This is a very important process in areas with high orographic forcing such as e.g., Öraefajökull.
Since HARMONIE-AROME has a higher resolution and a more complex representation of the atmosphere, which lead to a more accurate simulations of the meteorological variables, incoming mass and energy balance fluxes from HARMONIE were chosen to force the HIRHAM5 snow pack model in order to get the best representation of the incoming mass and energy fluxes for realistic runoff simulations.Using NWP models to force a snow pack scheme has previously been successfully done in studies by e.g., Bellaire and Jamieson [58] and Vionnet et al. [59].

Accumulation Scaling
Although the accumulation of precipitation in HARMONIE is an improvement over the HIRHAM5 results, there are still substantial biases in the winter balance compared to observations.Accurate estimation of the winter accumulation is especially important in the ablation zone, as too much modeled winter snow will delay the timing of the exposure of the underlying low albedo ice surface and thus cause underestimation of the amount of runoff.To get the most accurate representation of the snow thickness, a scaling is applied to each model year.The scaling is based on the average winter balance difference map from 1992 to 2015, as shown in Figure 3c.The same absolute scaling is used for every model year and is applied to each snow fall event.The scaling matrix is shown in Figure 4a.A new spin-up of the subsurface was performed by integrating the model for 100 years repeating the corrected forcing from 1980.Only the incoming precipitation was corrected, all other incoming fluxes were unchanged HARMONIE forcing.A new model run for 1980-2015 was then conducted.A comparison of the winter mass balance observations with the winter mass balance simulated with the corrected forcing is shown in Figure 4b.The largest difference from the uncorrected model comparison is that the Öraefajökull points are now much closer to the one-to-one line and the RMSE has now decreased to 0.37 m w.eq.This is almost within the error of the observations which is estimated to be 0.3 m w.eq.Before the correction, the too high winter accumulation in the ablation zone resulted in a delay of the timing of the albedo drop from snow to ice in the model, and in some instances the ice surface was not exposed at all (not shown).

Comparison to Runoff Measurements
Model comparisons with runoff observations from Skaftá GWC into Skaftá river and from Brúarjökull into Jökulsá á Dal river at Hjarðarhagi  and Hálslón reservoir (2007-2015) are shown in Figure 5 and Table 2.In Skaftá, the average deviation in total summer runoff over the time period is −0.07 km 3 , which is within the uncertainty of the compared time series.The correlation between model and data is also high at 0.78 for the whole period.
For the runoff from Brúarjökull, the deviation in total runoff between observations and model is at −0.6 km 3 .This higher deviation is expected, both because no correction has been added to account for non-glacierized areas and because the runoff is estimated from the water level in the Hálslón reservoir for part of the time series, where the uncertainty is expected to be larger than for the conventional gauging from Skaftá or Hjarðarhagi.The correlation between model and data is still high at 0.89 for the whole period, which suggests that the variability is well captured.For this study, the focus is generally on the runoff over the entire melt season.The snow pack model used here does not include any hydrological routing of runoff, and therefore a daily comparison of has larger deviations than a seasonal comparison.

Snow Thickness on 1 April
The conditions of the ice cap in spring will affect the summer runoff.Although the temperature in the winter snow is of importance, by far the most important variable is the thickness of accumulated snow on the ice in the ablation zone.This is because the snow thickness in the ablation zone determines the timing of the exposure of the underlying low albedo ice.To investigate snow thickness variability in 1980-2015, time series and histograms of the modeled average snow thickness in the ablation zone on April 1st for the whole ice cap as well as for selected outlet glaciers are shown in Figure 6.For the whole glacier, the modeled snow thickness (for ablation zone points below 1300 m) is found to have an average value of 1.2 m w.eq., with a 1.3 m w.eq.difference between the minimum and maximum value.However, since there might be a large variation in snow thickness for the different outlet glaciers, we focus on three outlet glaciers: the north-facing Brúarjökull with an average of 1.4 m w.eq., the south-west facing Skaftá GWC with an average of 1.2 m w.eq.and the south-facing Síðujökull with an average snow thickness of 1.1 m w.eq.The variation over the time series is 1.5-2.1 m w.eq.depending on the outlet.
All four histograms have a two-peaked distribution, which is especially apparent in the histogram for Brúarjökull.This is due to a shift in snow thickness in the mid-1990s.Around 1995, there was a rapid increase in ocean temperatures off the southern coast of Iceland which was likely due to atmospheric and ocean circulation changes around Iceland [60].This resulted in an increase of the mean annual temperature of ∼1 • C after 2000 compared to the mid-1990s (e.g., [8,61]), and observations of the annual surface mass balance show a clear shift from positive to negative values in this period.From the time series of the snow thickness, it is clear that the model catches this shift.

Sensitivity Runs
The HIRHAM5 snow pack model was run to create 36 different simulated spring conditions, i.e., accumulated winter snow, each followed by 36 different summer conditions to assess the sensitivity of the runoff to the thickness of the spring snow.Maps of the results from the sensitivity runs are shown in Figure 7, and time series and histograms of the results for the whole icecap and for selected outlet glaciers are shown in Figure 8.Only the ablation zone is shown in all figures as the different spring conditions are mostly going to affect the runoff from the ablation zone.The equilibrium line altitude (ELA), and therefore the ablation zone, varies both with location and between years, but has not been observed higher than 1300 m for the years under consideration.We have therefore defined the ablation zone as everything below 1300 m.For evaluation of the sensitivity run results, we introduce two terms: "variation in runoff for varying spring conditions std" σ spring and "variation in runoff for varying summer weather std" σ summer .Imagine a matrix of all 36 × 36 model combinations, where the rows are the 36 summers and the columns are the 36 spring conditions.Each point in the matrix is the total runoff over a summer from April to October.For σ spring , a standard deviation across each row is calculated.Since the summer weather is constant within each row while the spring conditions varies, this parameter is used to determine the changes in runoff due to changing spring conditions.σ summer is similar, but this time a standard deviation across each column is calculated.In this case, the runoff variations due to changes in the summer weather can be assessed.Figure 7a shows the mean of σ spring .The largest sensitivity to spring snow thickness is found in the lower parts of the north-facing glaciers Brúarjökull and Dyngjujökull.The south and east facing glaciers generally have a smaller but still significant sensitivity to changes in snow thickness, while the west facing glaciers only appear to have a small sensitivity to the spring conditions.
Figure 7b shows the mean of σ summer .In many parts of the ablation zone, the summer weather clearly has a larger influence on the runoff variability than the spring snow thickness.There are several controlling factors for the variations in the summer weather.In addition to the amount of incoming solar radiation, an important factor for the summer melt on Vatnajökull is the amount of summer snow.This is an especially important factor for Brúarjökull, which is more likely to receive precipitation from cold northern winds and has a high-elevation ablation zone.This high dependence on summer snow could contribute to the high variability due to summer weather found in the simulations over the ablation zone of Brúarjökull.Figure 7c shows the difference between the two different maps (a,b).Except for parts of the north-facing glacier outlets, where the runoff variation due to changes in spring conditions and summer weather is approximately equal, and in a few points on two south-facing outlets, the runoff variability is much more affected by variations in summer weather than variations in spring snow thickness.For the whole ice cap, changing the spring snow thickness amounts to a maximum of 22% change between the minimum and maximum runoff over a summer.Individual outlet glaciers have a larger sensitivity to the spring snow thickness, and e.g., the runoff from Brúarjökull increases with a maximum of 72% compared to the minimum runoff.
On average, σ spring for the whole ice cap is equal to 31% of σ summer , with a maximum of 46%.The variability in runoff due to changes in summer weather is therefore much higher than the variability due to changing spring conditions.For some individual outlet glaciers, the contribution of changing spring conditions is larger.For Brúarjökull, the average σ spring is 50% of the σ summer , with a maximum of 79%.
These results are also illustrated in Figure 8, which shows the time series for the whole glacier and for selected outlet glaciers for all model runs.Relevant statistical values for this figure are given in Table 3.The grey area shows the span of the runoff between the different runs and the red line shows runoff when the spring conditions of the specific year are used.Three outlet glaciers are shown, Brúarjökull, due to its high sensitivity to spring conditions, Síðujökull, which has a slightly smaller sensitivity, and Skaftá GWC, which has a low sensitivity.To get a range on the effect of the spring conditions over the entire time series, the best conditions for runoff are chosen for every year, i.e., the years with the least spring snow, and the worst conditions for runoff are chosen for every year, i.e., the years with the most spring snow.For Vatnajökull, the total runoff when the best conditions are chosen is ∼19% higher than the runoff when the worst conditions are chosen.For Brúarjökull this difference is ∼53%, for Síðujökull it is ∼27%, and for Skaftá GWC only ∼14%.
Table 3.The mean runoff of all the sensitivity runs (SR), the mean runoff when the maximum and minimum runoff is chosen for the whole time series and the percentage increase in runoff between the maximum and the minimum.σ spring is the mean standard deviation in runoff when the spring conditions are varied and σ summer is the mean standard deviation in runoff when the summer weather is varied.The mean runoff of the thick snow (TS) run, and the difference in runoff between the original simulation and the thick snow simulation are also given.

Mean SR Min SR Max SR SR Diff σ spring σ summer
Mean TS TS Diff (m w.eq.) (m w.eq.) (m w.eq.) (%) (m w.eq.)The black lines in Figure 8 show the simulated amount of runoff for each year if there is a thick snow thickness that does not disappear during the summer.Comparing the thick snow run with the original runoff simulations (1980-2015 with accurate spring snow thickness), there is a 14% increase in runoff from Vatnajökull averaged over the entire period since 1980 when exposed ice is allowed, with a maximum amount of 31% in 2005.Individual outlet glaciers have a larger sensitivity to the spring conditions, as for Brúarjökull the simulations estimate that the runoff is increased with an average 27% over the entire period when ice is exposed compared to the thick snow run, with a maximum amount of 77% in 2005.
In addition to the sensitivity studies, the results of the experiment can be used as a statistical method to estimate the most likely runoff from the modeled outlet glaciers.Since snow thickness and summer weather from the last 36 years have been used, all these runs can provide an estimate of the most likely values of runoff based on present climate.Probability distributions for the whole ice cap and the three chosen outlet glaciers are shown in Figure 8b,d,f,h.These show the most likely runoff from the glaciers based on the present climate, as well as the lowest and highest expected values.However, it is important to note that the effect of volcanic eruptions and dust in the surface from sources outside the glacier [62,63], as well as glacial runoff in the form of jökulhlaups, are not considered with this evaluation.In these instances, the runoff can be much higher than the modeled distributions.This method therefore works best for an estimation of the lowest runoff from a glacier based on the current climate.
Another important aspect of hydrological monitoring is the timing and onset of melt. Figure 9 shows the percentage of the glacier area which is experiencing melt for the whole ice cap and the same three outlets discussed previously.As an example, the earliest and latest date that 20% of the area has experienced melt onset, defined as having experienced more than 0.1 mm w.eq. of melt over five consecutive days, is shown with a stippled black line.Figure 9a,c,e,g show the melt area when the spring conditions are held constant using the simulated conditions from 2000, where the colored lines are the results using 1980-2015 summer forcing.The difference between when 20% of the area has experienced melt onset is 56-62 days depending on the outlet, and the difference in melt area between the different summers can be up to 80-98% of the area, depending on the outlet.Figure 9b,d,f,h show the melt area of all the runs, i.e., with both changing summer forcing and spring conditions.The changing spring conditions shift the dates when 20% melt onset occurs, but once again the effect is not as large as that of changing summer forcing.The difference between when 20% of the area has experienced melt onset is increased with 7-17 days, depending on the outlet.
To validate if the area of the ablation zone is accurate in the model and the variation in ELA between the SR, we investigate how often each model point becomes snow-free.Figure 10 shows the fraction of the 1296 runs in each grid point which becomes snow-free at some time over the summer.The location of the ablation is realistic.The ablation zone points always become snow-free, while around the equilibrium line it depends on the weather and initial conditions.The location of the equilibrium line is generally consistent with observations.This can e.g., be seen for Brúarjökull, where the equilibrium line has been at approximately 1200 m since 1996 on the western part of the glacier, whereas on the eastern part it has been at approximately 1100 m in the same period.The simulations of the eastern part of Brúarjökull are consistent with this, while for about 60% of the model runs the ELA is slightly lower than observed on the western part.Having determined that the location of the snow-free area is realistic, we also investigate how much the varying spring conditions advanced or delayed the timing of exposure of ice in each grid point.Figure 11 shows the mean spread (max-min) in the day when the ice surface is exposed when the summer weather is held constant while the spring conditions are varied.The largest spread is found in the lower part of Brúarjökull and Dyngjujökull, as well as the middle of Skeiðarárjökull.These areas all become snow-free every year and have a large dependency on the spring snow thickness.There is a fairly large average deviation between the first ice-free day when the spring conditions are varied, as most points have an average spread of 1-3 months between the first day the surface becomes snow-free.

Modelling of Runoff
While the simulated runoff from Skaftá GWC and Brúarjökull was captured well by the HIRHAM5 snow pack model in combination with HARMONIE-AROME meteorological forcing, with correlations of 0.78 and 0.89 and RMSEs of 0.2 km 3 and 0.7 km 3 (Table 2), comparisons to observations on a daily time scale yield larger deviations.The HIRHAM5 snow pack model is not a hydrological model, and therefore does not have any scheme for routing of runoff.The only delay taken into account by the model is a runoff time scale based on the surface slope of the glacier.The structure of the path of the melt water determines how efficiently the water is expelled from the glacier.There are many processes which control the path, and the efficiency of water transport is very variable over the runoff season.For example, water flow through snow is significantly slower than flow over bare ice, and a thick snow pack therefore contributes to runoff delay (e.g., [64,65]).The presence and number of moulins and crevasses to transport water to the glacier bed (e.g., [65,66]) and the timing and shape of subglacial channels (e.g., [67]) also affect the timing of runoff to the rivers.
Since the runoff delay is not considered in the simulations, the daily runoff is more variable in the simulation than in the observations.Figure 12a shows an example from 2012 of the daily modeled and observed runoff.Compared to the seasonal means given in Table 2, the daily runoff has a much higher RMSE compared to the mean and a lower correlation (Table 4).To investigate the runoff delay, a 1-day shift in the observations was added to account for a runoff delay in the order of 12 h.This improved both the correlation, RMSE and NSE for both glacial rivers.Smoothing the time series but using a 7-day mean is another way to account for the delay.As shown in Figure 12b, a 7-day average improves the fit between observations and simulations.The Nash-Sutcliffe model efficiency coefficient (NSE) is also improved.The better fit for the 7-day average is also clearly showcased in the boxplot of the prediction error showed in Figure 12c.The 7-day average smoothes the time series which results in fewer outliers from the distribution.The quantiles and the median are also closer to zero in the weekly plot than in the daily.

Sensitivity of Runoff to Winter Snow Cover
One of the aims of this study was to quantify the sensitivity of runoff to the spring snow thickness for Vatnajökull and investigate whether the summer runoff could be inferred from the spring snow thickness.As shown in Figures 7 and 8, the summer runoff does not have a high dependency on the spring snow thickness, with the magnitude of the effect varying greatly between outlets.Although a large snow thickness will of course lead to less runoff than a small snow thickness, the weather conditions during the summer have a much larger effect on the runoff, as the effect of changing spring snow thickness is on average 31% of the variation due to changing weather conditions during the summer.
While the larger effect of changes in summer weather is not surprising, the advantage of the method used in this study is that we can quantify the importance of the spring snow thickness, and how much of the variability in runoff over the last 36 years is due to changes in summer conditions and how much is due to changes in winter accumulation.A study like this using over a thousand potential summer variations has not previously been done and provides a sound statistical background for assessing the sensitivity.For the entire glacier, we could estimate the amount of runoff due to the exposed dark ice (on average 14% for Vatnajökull since 1980) and how much of the runoff variability is due to changes in summer climate and how much due to changes in winter accumulation.The model runs also provide an estimate for individual outlet glaciers, as e.g., Brúarjökull has a large sensitivity to changes in spring snow thickness.According to the simulations, an estimated average of 27% of the runoff since 1980 is due to exposed dirty ice on Brúarjökull.In addition, the variation in runoff for Brúarjökull due to changing spring conditions is 50% of that due to changing summer weather.

Error Sources
HARMONIE-AROME meteorological forcing was chosen to drive the HIRHAM5 snow pack model to reduce the errors, but comparison with available AWS observations showed clear model biases.The incoming shortwave radiation still has an average overestimation of 2.4 W m −2 when compared to AWS observations, while the incoming longwave radiation and turbulent fluxes have an average underestimation of −13 W m −2 and −8.9 W m −2 , respectively.The outgoing longwave simulated by the snow pack scheme has an average underestimation of −1.3 W m −2 , while the outgoing shortwave is overestimated by 10 W m −2 .The deviation in outgoing shortwave radiation is partly due to an overestimation of incoming shortwave and precipitation, and partly due to the albedo scheme used in the snow pack model.The total deviation in energy available for melt is therefore approximately −28 W m −2 when compared to available AWS observations from May to September.This is equal to an underestimation in melt of approximately 0.7 cm/day.The meteorological input contributes a larger error than the snow pack simulations.This general underestimation occurs even though the snow pack model and input data have been tuned to try to best fit available observations, as they are introduced into the simulations e.g., due to uncertainties in model structure, chosen parameters, and input data.One way we tried to correct the input data was by precipitation scaling by assume that the precipitation bias is time-invariant, and the scaling therefore can be applied to the entire time period.However, this assumption is only valid for short-term simulations, and the correction should therefore not be used for model runs of longer time scales.It is only applied to get the best possible simulation of the winter balance for the present day, so that the simulated spring conditions are as realistic as possible.In the mid-1990s there was a shift in the specific mass balance of Vatnajökull from positive to negative due to a rapid increase in the ocean temperature south of Iceland [8].Few measurements of the mass balance before 1996 are available, and therefore the correction might not be applicable for the period before the shift.Nevertheless, the runoff from Skaftá GWC shows an equally good fit between model and measurements for the 1984-1996 period as for the 1996-present period.
Another error source is in the model setup.While running the snow pack scheme offline is a computationally cheap, fast, and flexible option, using this method ignores the feedbacks between the glacier conditions and the atmospheric circulation, such as e.g., the albedo and temperature.When forcing the snow pack scheme with HIRHAM5 incoming fluxes, the error due to neglected feedbacks is expected to be small.This is because only the refreezing [51] and the albedo [10] schemes have been altered from the online run, and changes in feedbacks due to these alterations should be small.On the other hand, running the snow pack scheme with a completely different model could lead to larger errors due to neglected feedbacks.For example, neither sublimation nor evaporation are calculated by the snow pack scheme, which affect the total water balance.However, since Vatnajökull is a temperate ice cap, sublimation and evaporation are generally assumed negligible.

Conclusions
In this study, we show that we can provide realistic estimates of the surface mass balance and runoff from Vatnajökull by combining the snow pack scheme from HIRHAM5 with the meteorological forcing from HARMONIE-AROME.HARMONIE-AROME was chosen as the forcing due to its complex atmospheric model using hydrostatic physics and prognostic precipitation, which simulate the incoming energy and mass fluxes needed to force the HIRHAM5 snow pack scheme better than atmospheric model in HIRHAM5.
The combined model can be used to simulate runoff on seasonal timescales, with high correlation (0.78 and 0.89) and small difference (16% and 29%) when compared to runoff observations from two rivers with a large glacial component.For this study, we were mostly interested in the total runoff from the glacier over the melt season and therefore not the routing of runoff.Coupling the snow pack scheme to a hydrological model which includes routing of runoff would be the next step to improve the comparison on short timescales.Additional work is also needed to simulate the runoff from ice-free zones, which would also improve the comparison.
The sensitivity of runoff to the spring snow thickness is evaluated by running a set of SR, using all the spring conditions from 1980 to 2015 as a starting point for runs over all the summers in the same period.We find that the variability in summer weather has a larger effect on the runoff than the variability in spring snow thickness, with an average variation in runoff due to changes in snow cover equal to only 31% of the variation due to summer weather for the whole ice cap.However, the spring snow thickness is still an important factor, especially for north-facing outlets, as e.g., Brúarjökull has an average variation in runoff due to changes in snow cover equal to 51% of the variation due to summer weather.
The range in values when the most extreme spring conditions are chosen are a good estimation of the importance of the spring conditions.The total runoff from Vatnajökull when the thinnest snow is chosen is ∼19% higher than the runoff when the thickest snow is chosen.For north-facing Brúarjökull this difference is ∼53%, for south-facing Síðujökull it is ∼27%, and for south-facing Skaftá GWC only ∼14%.Once again, we can conclude that while the spring snow cover is important for Brúarjökull, it is a less important factor for the south-facing outlets.
An additional simulation was conducted using thick enough snow so that the underlying ice surface is never exposed.This run was used to evaluate how much of the runoff from 1980 to 2015 was due to exposed dirty ice.We estimate that there was approximately a 14% increase in runoff from Vatnajökull due to the exposure of dirty ice from 1980 to 2015, a 27% increase from Brúarjökull, a 12% increase from Síðujökull, and a 7% increase from Skaftá GWC.As the climate warms and the dirty ice in the ablation zone is exposed earlier in the melt season, the amount of melt due to exposed ice is expected to increase.In these experiments, using only conditions simulated for the present climate, a maximum of 31% increase in runoff from Vatnajökull due to the exposure of dirty ice was found, while a maximum increase of 77% was simulated for Brúarjökull.The timing and duration of exposure of dirty ice can therefore have a large contribution to the amount of melt, and further investigations are needed to estimate the effect on the future mass balance.This could e.g., be done by using output from CORDEX simulations of the future climate to force the snow pack scheme.However, forcing the snow pack scheme with different models at lower resolution and not forced by reanalysis data will lead to new and larger biases that will need to be investigated and corrected before the model can be confidently applied.

Figure 1 .
Figure 1.Location of AWS stations (black triangles) and surface mass balance sites (black dots) used in this study.The red lines show the outlines of the ice divides, while the bluelines show the outlines of the Skaftá (west) and Hálslón (east) glacial water catchments.The symbols outside the ice cap show the water gauges (wg).

Figure 2 .
Figure 2. A schematic diagram of the coupling of HARMONIE-AROME and HIRHAM5.The precipitation scaling is discussed in Section 4.2.

Figure 3 .
Figure 3. (a) HARMONIE and (b) HIRHAM5 winter mass balance compared to measurements from 1992 to 2014.Color bar shows how many points are in the same spot and purple ellipsis outlines 68% of the points.Statistics on the regression lines are given in Table 1.Interpolated observed mean winter mass balance maps from 1992 to 2014 compared to (c) HARMONIE and (d) HIRHAM5 winter mass balance.

Figure 4 .
Figure 4. (a) The scaling matrix used for correction of each winter snowfall; (b) comparison between winter mass balance observations and the modeled corrected winter mass balance for 1992-2015.The RMSE is 0.4 m w.eq.and the correlation r = 0.89.Color bar shows how many points are in the same spot and purple ellipsis outlines 68% of the points.

Figure 7 .
Figure 7. Maps showing the effect of the summer weather and spring conditions on runoff from the ablation area (here set as below 1300 m).(a) the average varying spring conditions std (σ spring ); (b) the varying summer weather std (σ summer ); (c) the difference between the maps (a,b).

Figure 8 .
Figure 8.Time series of accumulated runoff for April-October from (a) Vatnajökull; (c) Brúarjökull; (e) Síðujökull; (g) Skaftá GWC for different spring conditions.The red line is the original run where the spring conditions of the specific year are used, the grey area shows the range in runoff when the spring conditions are changed, and the black line shows how much melt would occur if the underlying ice is never exposed.Probability distributions for the runoff created using the results of the sensitivity runs are shown for (b) Vatnajökull; (d) Brúarjökull; (f) Síðujökull; (h) Skaftá GWC.

Figure 9 .
Figure 9.Time series of melt area for (a,b) Vatnajökull; (c,d) Brúarjökul; (e,f) Skaftá GWC; (g,h) Síðujökull.The four left figures show the melt area when the spring conditions are held constant, using the conditions from April 2000.The right figures show all the simulated runs.The colored lines are different summers from 1980 to 2015.The stippled lines show the earliest and latest day that 20% of the area has had an onset in melt.

Figure 10 .
Figure 10.The fraction of the 1295 model runs where the grid points have exposed ice at some time during the summer.

Figure 11 .
Figure 11.The average spread (max-min) in DOY when each ice point becomes snow-free when the summer weather is kept constant and the spring conditions are altered.

Figure 12 .
Figure 12.(a) Daily and (b) weekly mean observed and simulated runoff from Skaftá GWC; (c) Box plots of the prediction error for daily and weekly comparisons for the entire time series.

Table 1 .
Evaluation of the simulated mass balance and energy flux components used to force the snow pack model, i.e., winter mass balance (a measure for precipitation), incoming short-and longwave radiation and turbulent fluxes, in HIRHAM5 and HARMONIE.The mean values of the average difference (model-observations), RMSE, and correlation (r) are given.The turbulent fluxes are a model-model comparison.

Table 2 .
Comparison statistics of the average summer runoff (April-October) from Skaftá GWC into

Table 4 .
Comparison statistics of the runoff (April-October) from Skaftá GWC into Skaftá river and from Brúarjökull into Jökulsá á Dal river at Hjarðarhagiand Hálslón reservoir(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)according to observations and the simulation.The mean difference, RMSE, correlation, and Nash-Sutcliffe model efficiency coefficient are given.The RMSE, correlation, and NSE are calculated both for daily means, daily means with a model shift of one day, and weekly means.