Evaluation of MODIS Albedo Product over Ice Caps in Iceland and Impact of Volcanic Eruptions on Their Albedo

Albedo is a key variable in the response of glaciers to climate. In Iceland, large albedo variations in the ice caps may be caused by the deposition of volcanic ash (tephra). Sparse in situ field measurements are insufficient to characterize the spatial variation of albedo over the ice caps. Here we evaluate the latest MCD43 MODIS albedo product (collection 6) to monitor albedo over the Icelandic ice caps using albedo from ten automatic weather stations in Vatnajökull and Langjökull as ground truth. We examine the influence of the albedo variability within MODIS pixels by comparing the results with a collection of Landsat scenes. The results indicate a good ability of the MODIS product to characterize the seasonal and interannual albedo changes with correlation coefficients ranging from 0.47 to 0.90 (median 0.84) and a small bias ranging from -0.07 to 0.09. The root-mean square errors (RMSE) ranging from 0.08 and 0.21, is larger than that from previous studies, but we did not discard the retrievals flagged as bad quality to maximize the amount of observations given the frequent cloud obstruction in Iceland. We find a positive but non-significant relationship between the RMSE and the subpixel variability as indicated by the standard deviation of the Landsat albedo within the MODIS pixel (R=0.48). The summer albedo maps and time series computed from the MODIS product show that the albedo decreased significantly after the Eyjafjallajökull and Grímsvötn eruptions in 2010 and 2011 in all the main ice caps (except the northernmost Drangajökull), with albedo reduction up to 0.6 over large regions of the accumulation areas. Following this validation, these data will be assimilated in an energy and mass balance model of to better understand the relative influence of the volcanic and climate forcing to the ongoing mass losses of Icelandic ice caps.


Introduction
The albedo of a surface is defined as the fraction of the incoming solar radiation that is reflected by this surface.It is a key variable for calculating a glacier's surface energy balance and thus its mass balance since it defines the amount of solar radiation that is absorbed at the glacier surface and available for melt [1,2].For instance, the role of albedo in the Greenland ice sheet mass balance has received much attention since its evolution could greatly affect Greenland´s contribution to sea level rise in addition to that caused by atmospheric warming alone [3][4][5][6].In contrast, glaciers in Iceland hold only a 1 cm potential contribution to sea level rise [7] but they can undergo very large albedo variations due to the deposition of volcanic ashes during volcanic eruptions or by the blowing of ashes from previous eruptions during storms [8].The volcanic ash is produced during explosive eruptions of the volcanoes in Iceland and can be transported by wind over tens of kilometers, depending on the atmospheric conditions and the particle properties [9][10][11].Also, buried layers of tephra from past eruptions may become exposed by glacier ablation.Since the ash is generally made of dark-colored basalt particles that strongly contrast with the snow and firn surfaces (at least in the visible portion of the spectra), this can cause strong seasonal and inter-annual variations in the albedo of Icelandic glaciers.Albedo observations are therefore important to understand and quantify a glacier's response to climate variability.
Recently, two major events caused significant ash deposition on Icelandic glaciers.The first event is a sequence of eruptions at Fimmvörðuháls from 20 March to 12 April 2010 followed by Eyjafjallajökull volcanoes from 14 April to 22 May 2010 [10].These eruptions ejected large amounts of fine dark trachyandesite tephra particles, leading to the closure of the airspace in many west European countries but also to the darkening of the Icelandic glaciers due to ash falls.The 2010 summer ablation of Langjökull and western Vatnajökull ice caps was up to three times higher than the average during the preceding warm decade [12].A year later, the Grímsvötn volcano erupted from 21 to 28 May 2011 and ejected basalt particles, also causing ash deposition and air travel disruptions, but mostly in Iceland due to the direction of the prevailing wind and the coarser size of the ash particles in comparison with the Eyjafjallajökull eruptions [13].
In collaboration with the National Power Company, the Institute of Earth Science at University of Iceland has installed a number of automatic weather stations (AWS) on the Langjökull and Vatnajökull ice caps.The AWS instruments include upward and downward looking shortwave radiometers on the ice caps that allow direct measurement of the albedo at the station locations (Section 2).However, the stations provide only point measurements that do not capture the spatial variability of the albedo over the whole extent of the ice caps.Optical satellite remote sensing techniques provide the opportunity to complement these data with spatially distributed observations.In particular mid-resolution sensors with a large swath, like the MODerate resolution Imaging Spectroradiometer (MODIS), have a daily repeat cycle that enables the reconstruction of the land surface albedo with a temporal resolution that is compatible with the rapid changes of the ice cap surfaces.Among the remote sensing albedo products, the MODIS Bidirectional Reflectance Distribution Function (BRDF)/albedo products are probably the most widely used [14].The MCD43A3 product [15] provides measurements of the global land surface spectral and broadband albedo on a 500 m resolution grid based on 16 days of multi-angular observations.MODIS data has been used to study the effect of the 2004 eruption in Grímsvötn on the albedo of Vatnajökull ice cap and how it evolves after the eruption [8].Previous evaluation studies of the MODIS albedo product in glacier environments focused on the Greenland ice sheet [16,17].The same MODIS product was also used to analyze inter-annual albedo changes over the Greenland ice sheet [5].However, as explained above, the albedo of Icelandic glaciers is expected to have larger temporal and spatial variability than that of the Greenland ice sheet.The daily revisit capability of MODIS is in theory appropriate to monitor the rapid changes of albedo due to ash deposition, snow fall, or snow melt [18], but the actual frequency of observations per pixel can be strongly reduced by cloud cover.Clear skies are rare in Iceland due to the persistent cyclonic activity in the North Atlantic known as the Icelandic low [19].The MODIS resolution at nadir is close to 500 m, which is considerably larger than the typical horizontal variability of albedo, especially in the ablation area (Figure 1).The previous version of the 16-days MODIS albedo (v005) albedo was evaluated in the period 2000-2008 using two stations on Tungnaárjökull, the western outlet glacier of Vatnajökull (AWSs tunef and tunmi, Figure 2) [8].The RMSEs with respect to the in situ albedo for the two stations were 0.12 and 0.09 with no systematic deviations.Since then, a new version of the MCD43 product (V006) has been generated based on the reprocessed collection 6 level 1 MODIS data (released in 2015).Collection 6 includes the latest improvements made by the MODIS science team since the release of collection 5 in 2008 [20].This update is of relevance for albedo studies since it is designed to remove a long-term drift in the visible and near-infrared bands caused by the aging of the MODIS sensor aboard the Terra satellite that was observed in the Collection 5 dataset [21,22].Also, relevant for this study, is the improvement of the geolocation grid based on a new terrain correction and updated sensor geometric parameters [23].The cloud masking algorithm was also improved, but the expected gain in the number of observations mainly affects desertic and very humid tropical areas [24].MCD43 version V006 achieved a stage 3 validation status according to the NASA MODIS land team, which means that "its accuracy is well less than 5% albedo at the majority of the validation sites studied thus far, and even those albedo values with low quality flags have been found to be primarily within 10% of the field data.Data for solar zenith angles greater than 70 degrees should be considered suspect."(Validation Status for MODIS BRDF/Albedo: MCD43, Product version: Collection 5/6, November 2015).
Here, we present an evaluation of the latest MODIS albedo product (MCD43A V006) in Iceland.Our first motivation is to assess the reliability of this product for future assimilation in an energy and mass balance model of the Icelandic ice caps.Therefore, given the persistent cloud coverage in Iceland, we started by analyzing the amount of missing data in the MCD43 albedo product (Section 3.1).Then, we used continuous incoming and outgoing shortwave radiation measurements made at 10 AWSs in the period 2001-2012 on Vatnajökull and Langjökull ice caps to assess the product accuracy (Section 3.2).In addition, higher resolution images from Landsat were used to assess whether subpixel variability can help explain the discrepancies between the in-situ and MODIS data (Section 3.3).Finally, albedo maps for Iceland are generated to quantify the effect of dust deposition on the albedo of the ice caps in the aftermath of the two aforementioned volcanic eruptions (Section 3.4).The previous version of the 16-days MODIS albedo (v005) albedo was evaluated in the period 2000-2008 using two stations on Tungnaárjökull, the western outlet glacier of Vatnajökull (AWSs tunef and tunmi, Figure 2) [8].The RMSEs with respect to the in situ albedo for the two stations were 0.12 and 0.09 with no systematic deviations.Since then, a new version of the MCD43 product (V006) has been generated based on the reprocessed collection 6 level 1 MODIS data (released in 2015).Collection 6 includes the latest improvements made by the MODIS science team since the release of collection 5 in 2008 [20].This update is of relevance for albedo studies since it is designed to remove a long-term drift in the visible and near-infrared bands caused by the aging of the MODIS sensor aboard the Terra satellite that was observed in the Collection 5 dataset [21,22].Also, relevant for this study, is the improvement of the geolocation grid based on a new terrain correction and updated sensor geometric parameters [23].The cloud masking algorithm was also improved, but the expected gain in the number of observations mainly affects desertic and very humid tropical areas [24].MCD43 version V006 achieved a stage 3 validation status according to the NASA MODIS land team, which means that "its accuracy is well less than 5% albedo at the majority of the validation sites studied thus far, and even those albedo values with low quality flags have been found to be primarily within 10% of the field data.Data for solar zenith angles greater than 70 degrees should be considered suspect."(Validation Status for MODIS BRDF/Albedo: MCD43, Product version: Collection 5/6, November 2015).
Here, we present an evaluation of the latest MODIS albedo product (MCD43A V006) in Iceland.Our first motivation is to assess the reliability of this product for future assimilation in an energy and mass balance model of the Icelandic ice caps.Therefore, given the persistent cloud coverage in Iceland, we started by analyzing the amount of missing data in the MCD43 albedo product (Section 3.1).Then, we used continuous incoming and outgoing shortwave radiation measurements made at 10 AWSs in the period 2001-2012 on Vatnajökull and Langjökull ice caps to assess the product accuracy (Section 3.2).In addition, higher resolution images from Landsat were used to assess whether subpixel variability can help explain the discrepancies between the in-situ and MODIS data (Section 3.3).Finally, albedo maps for Iceland are generated to quantify the effect of dust deposition on the albedo of the ice caps in the aftermath of the two aforementioned volcanic eruptions (Section 3.4).

MODIS Albedo
We used the latest version of the MODIS (Aqua/Terra combined) MCD43A albedo product V006 [15] (released in late 2015).Only MODIS tile h17v02 was used, as it covers the entire study area.We extracted the black-sky and white-sky broadband shortwave albedo values in every MODIS granule (one per day) between 1 January 2001 and 31 December 2012 at the 10 AWS locations using the gdallocationinfo utility of the Geospatial Data Abstraction Library [25].The black-sky albedo is the directional-hemispherical reflectance, i.e., albedo in the absence of diffuse illumination, while the white-sky albedo is the bi-hemispherical reflectance, i.e., albedo in the absence of direct illumination.We accounted for the horizontal displacement of the AWS due to glacier flow by using coordinates that were measured every year in summer at each AWS (see Section 2.2).Missing annual coordinates values (8%) were filled by the average of the other years' coordinates for each AWS.As recommended by the data provider [15,26], the data were discarded if the solar zenith angle is higher than 70° at the time of the latest MODIS acquisition.Before 02 July 2002, the MODIS sensor was only onboard the Terra platform, hence the acquisition time is around 10:30 in local solar time (LT).After that date the acquisition time is given by Aqua overpass time (13:30 LT).The data that were flagged as bad quality retrievals in the MODIS product ancillary files were not removed to maximize the number of observations.These quality flags are very conservative and lead to a strong reduction of the number of observations if flagged data are removed (about ten times less values, up to zero data at one station).This would make the data barely useful for our long-term goal, which is to provide the albedo maps in a surface energy balance model applied on Icelandic glaciers.For the comparison with AWS data, we computed the actual or "blue-sky" albedo as the average between the black-sky albedo and the white-sky albedo assuming a constant fraction of diffuse illumination as done by [8].In any case, the differences between the two albedo values are very low in high-latitude regions [8,17].In our study area the relative differences are lower than 5%.
Then, the broadband shortwave black-sky albedo product was used to study the effects of the Eyjafjallajökull 2010 and Grímsvötn 2011 eruptions on the summer ice caps surfaces (Section 3.4).Every granule of tile h17v02 for the summer months (1 July to 1 September) from 2001 to 2012 were re-projected from sinusoidal grid to the WGS84 UTM 28N coordinate system.This was done with the NASA's MODIS reprojection tool without changing the original resolution (approximately 500

MODIS Albedo
We used the latest version of the MODIS (Aqua/Terra combined) MCD43A albedo product V006 [15] (released in late 2015).Only MODIS tile h17v02 was used, as it covers the entire study area.We extracted the black-sky and white-sky broadband shortwave albedo values in every MODIS granule (one per day) between 1 January 2001 and 31 December 2012 at the 10 AWS locations using the gdallocationinfo utility of the Geospatial Data Abstraction Library [25].The black-sky albedo is the directional-hemispherical reflectance, i.e., albedo in the absence of diffuse illumination, while the white-sky albedo is the bi-hemispherical reflectance, i.e., albedo in the absence of direct illumination.We accounted for the horizontal displacement of the AWS due to glacier flow by using coordinates that were measured every year in summer at each AWS (see Section 2.2).Missing annual coordinates values (8%) were filled by the average of the other years' coordinates for each AWS.As recommended by the data provider [15,26], the data were discarded if the solar zenith angle is higher than 70 • at the time of the latest MODIS acquisition.Before 02 July 2002, the MODIS sensor was only onboard the Terra platform, hence the acquisition time is around 10:30 in local solar time (LT).After that date the acquisition time is given by Aqua overpass time (13:30 LT).The data that were flagged as bad quality retrievals in the MODIS product ancillary files were not removed to maximize the number of observations.These quality flags are very conservative and lead to a strong reduction of the number of observations if flagged data are removed (about ten times less values, up to zero data at one station).This would make the data barely useful for our long-term goal, which is to provide the albedo maps in a surface energy balance model applied on Icelandic glaciers.For the comparison with AWS data, we computed the actual or "blue-sky" albedo as the average between the black-sky albedo and the white-sky albedo assuming a constant fraction of diffuse illumination as done by [8].In any case, the differences between the two albedo values are very low in high-latitude regions [8,17].In our study area the relative differences are lower than 5%.
Then, the broadband shortwave black-sky albedo product was used to study the effects of the Eyjafjallajökull 2010 and Grímsvötn 2011 eruptions on the summer ice caps surfaces (Section 3.4).Every granule of tile h17v02 for the summer months (1 July to 1 September) from 2001 to 2012 were re-projected from sinusoidal grid to the WGS84 UTM 28N coordinate system.This was done with the NASA's MODIS reprojection tool without changing the original resolution (approximately 500 m) using the nearest neighbor resampling method.We did not use the bilinear or cubic resampling method to avoid smoothing the albedo transition between the snow cover and the bare ice surfaces.Also, the implementation of these methods increases the number of no-data values since a pixel will get no-data value if at least one adjacent pixel is no-data.
From these data we built annual summer albedo composite images by calculating, for each pixel and for each year, the mean values of all available observations between 1 July and 1 August (Section 3.4).During this period the glacier surfaces are the least affected by recent ash or snow falls.A two-month period of integration was needed to sufficiently reduce the large number of missing data due to cloud obstruction.The mean and standard deviation of the annual summer albedo values within the area of the six main ice caps was computed from each annual summer composites to study the temporal evolution of the albedo for each ice cap (Figure 2).The ice caps polygons were obtained from the Corine Land Cover 2006 in Iceland, which was performed by the National Land Survey of Iceland.We have used this database because it gives the ice caps area for the year 2006, which is in the center of our study period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).Last, we mapped the effects of the 2010 and 2011 eruption using the 2009 composite image as a reference.An albedo anomaly map was computed for both events by taking the difference between the reference and the post-eruption composites.

In Situ Data
The AWSs providing measurements to compute the full surface energy balance, have been operated on the Vatnajökull ice cap since 1996 and since 2001 on the Langjökull ice cap [27,28].The AWSs measure every 10 min the incoming and reflected solar radiation along with other meteorological variables.In spring the meteorological instruments are mounted on a mast that follows the surface while the ice melts throughout the summer.All instruments were compared during test setups in Reykjavík before installation.The AWSs were visited in summer to ensure that they are functioning properly.Here, we used the incoming and reflected shortwave measurements from ten AWSs during the period 2001-2012 (Figure 2, Table 1).Three of the ten stations are operated throughout the year; the others are installed in early spring and taken down in late autumn.The daily albedo was computed as the ratio of the upward to downward daily incoming and reflected shortwave radiation totals (daily integrated albedo as in [16]).Data in the late autumn and winter when the solar elevation is low cannot be used to calculate the surface albedo and were therefore excluded (Table 1).The radiative fluxes were measured using Kipp and Zonen CM14, CNR1 or CNR4 sensors (Table 1), all with relatively uniform spectral response over the range 0.3-2.8µm, and uncertainty that has been estimated to be 3-5% on an ice sheet surface [28,29].Exceptions are the stations hoffne, where the shortwave radiation was measured using LI-COR SPLite (incoming) and LI-COR (reflected) radiometers, and hoffef from 2001 to 2008 where LI-COR was used for both the incoming and reflected short wave radiation.The spectral response of the LI-COR instrument is not uniform over the full solar radiation range; it is very low at 0.4 µm, increases to 0.95 µm and decreases again to a cutoff at 1.2 µm [30].The accuracy is given as ±5% in daylight conditions when the glass is kept clean but it may cause large errors at a low solar elevation or in wet condition (e.g., raindrops will act as a lens).From the AWS setup and the given field of view of the instruments, the footprint of the in-situ albedo measurements is estimated to be within a circle of 10 m to 15 m in diameter.
Two of the 10 AWSs, tunef and bruef, are located high in the accumulation area (Ac) (Table 1).Normally the surface below these sensors (within the footprint of the albedo measurements) represents the surface conditions in the surrounding areas (confirmed in our early-, mid-and late-summer visits to the AWSs).The stations brumi, hoffef and langef are located close to the present ELA and the station tunmi is at an elevation of about 100 m below the average ELA of the last decade.The surrounding areas of tunmi and brumi are flat and therefore temporary water retention (slush) and water channels can form below the instruments and/or in the surrounding areas.This has been observed during some of the mid-summer visits to the stations (Figure 1).The stations hoffef and langef are located where the surrounding surface slopes are larger than at tunmi and brumi and there are crevasses that drain the surface runoff.It is not expected that the measurements at these stations are disturbed by water channels or slush.The four lowest stations are all close to the glacier terminus with low observed winter accumulation during the years 2001-2012 (within 0.5-1.5 m), and hence the previous year summer surfaces are normally exposed during the early melt season.The following should be noted for those stations:

•
AWS brune was located about 3.5-4 km from the terminus.In the observation period the surface has lowered by about 30 m.The surrounding area is characterized by volcanic ash layers emerging from the ice in the ablation zone as the surface melts.Generally, the ice under this AWS and within tens of meters from the station is covered with sand and dirt resulting in very low albedo values.

•
AWS hoffne was just above a steep glacier front, around 0.5-0.7 km from the terminus.In the observation period the surface has lowered about 40 m.The station and the nearest surrounding areas are assumed not to be greatly influenced by dirt blown from the nearby non-glaciated areas.However, the lowest 100-200 m of the glacier front is heavily debris covered.

•
AWS breid was above a steep glacier front, around 1-1.5 km from the terminus.From 2001 to 2012 the surface has lowered by about 180 m.The area around this station consists of dark dirt bands and is located close to a medial moraine.

•
AWS langne was moved in steps up the glacier during the observation period, it is located 3-4.5 km from the terminus, depending on the year.The lowest terminus area is heavily debris covered, as a consequence of sand and dirt blown from the outside areas.The ice areas above this dark terminus region are considerably brighter.This AWS was located within the debris region until 2006, with the debris reaching at most 100 m above the station according to a 2004 SPOT image.
As the glacier has been retreating, the station was moved higher up on the glacier in a few steps; it was located at the margin of the dark region and the cleaner ice above it in 2007, and then moved again in 2008 to a location of about 100 m above the terminus region.Hence, the in-situ albedo observations reflect a dirty surface until 2006/2007 and cleaner ice since 2008.

Landsat Data
To study the effect of the subpixel variability on the accuracy of the MODIS retrieval, we used the Landsat surface reflectance data from the United States Geological Survey.These images provide the surface reflectance at 30 m resolution for each Landsat spectral band after atmospheric correction.The surface reflectance values were converted to broadband albedo values using a narrowband-to-broadband albedo conversion for glacier ice and snow based on band 2 and band 4 [31].This empirical formula was obtained from airborne measurements performed over Vatnajökull and Greenland [31].
A collection of 15 Landsat 5 and Landsat 7 scenes of the study area, that have a limited cloud cover based on the inspection of their quicklooks in the USGS Earth Explorer portal [32], were downloaded.Then, color composites of bands 5, 3 and 2 (shortwave infrared, red and green) were generated at full resolution to visually determine, for every image, at which AWS it was safe to extract the reflectance values in a region of 500 m × 500 m without cloud contamination (Table 2).This band combination was efficient to distinguish between the clouds and snow surface.We found that this manual approach was more reliable than an automatic method that relied on the cloud mask provided with the Landsat data.Indeed, in this area the default cloud mask is inaccurate due to the mixing of snow, ice, fog and cloud surfaces.The computed Landsat albedo was sampled for each selected pair of AWS site and Landsat scene.As an extraction mask we used a square polygon of 16 × 16 Landsat pixels (480 m × 480 m) to simulate the approximate size of a MODIS pixel (500 m × 500 m).The standard deviation of the sampled values was computed as an indicator of the spatial variability.Pixels marked as saturated in band 2 or band 4 were considered as no-data in the calculation.The effective size of a MODIS pixel can be larger than 500 m depending on the date of the acquisition due to the bowtie effect of the MODIS sensor at high view angle [33].As a result, our method is expected to provide a lower bound of the albedo variability within a MODIS pixel.All these operations were made using the GDAL utilities [25] embedded in bash scripts.

Missing Data
The mean annual number of days without observation in the MCD43 product (including bad quality retrievals) during the period 1 January 2001 to 31 December 2012 was computed over the whole of Iceland (Figure 3).The average number of missing values is 49.9%, but it is higher over the ice caps, in particular in the central part of Vatnajökull, and in the north-west peninsula, Vestfirðir, where the snow probability is high [34].The number of missing values is also higher in the high-elevation areas like the Tröllaskagi peninsula in the north, where the glaciers are rather small but where the snow probability is also high.This suggests that the missing values are related to the snow cover.Another area with a high probability of no-data is Skeiðarársandur, the large outwash plain in front of the terminus of Skeiðarárjökull (main south flowing outlet glacier of Vatnajökull).masking algorithm.However, this should not affect the sandur area, which is made of numerous braided rivers and dark sand plain that have a very different spectral signature than the clouds.Rather, it may be due to a convergence failure in the MODIS albedo/BRDF algorithm in this area.In any case, these results suggest that the MODIS coverage is strongly reduced in Iceland and therefore further applications to data assimilation in a glacier model should rely on an algorithm that is robust to missing values.

Comparison within Situ Data
There is an overall good visual and statistical agreement between the MODIS shortwave broadband black-sky albedo and the in-situ albedo AWS measurements (Figures 4 and 5, Table 3).The seasonal variability in the in situ dataset is very high, ranging from 0.05 to 1 and is well captured by the MODIS product, as demonstrated by the high correlation coefficients (Table 3).When considering all the values, including low quality retrievals, the median of the correlation coefficients is 0.84, close to the value of 0.82 reported in Greenland [16].
The low consistency between the AWS and MODIS albedo at AWS hoffne and partly at AWS hoffef is probably partly due to the poor quality of the LI-COR radiometers used at those sites.Excluding AWS hoffne, there are still rather large differences during the summer at the lowest ablation stations breid, brune and langne.Field observations suggest that the differences at AWS brune and AWS breid might be explained with the resolution difference between the MODIS pixel and the AWS footprint (Section 3.3).When the winter accumulation has melted, the surrounding areas of AWS brune are characterized by dark volcanic fringes and at AWS breid by dirt bands and medial moraine for some years (Section 2.2).Apart from these explainable differences during the summer, a generally good agreement is obtained between the AWS and MODIS pixel albedo at breid, brune and langne (Table 3).The snow covered areas are prone to cloud misclassification.Therefore the high presence of no-data in the MODIS albedo product may be partly inherited from the upstream MODIS cloud masking algorithm.However, this should not affect the sandur area, which is made of numerous braided rivers and dark sand plain that have a very different spectral signature than the clouds.Rather, it may be due to a convergence failure in the MODIS albedo/BRDF algorithm in this area.In any case, these results suggest that the MODIS coverage is strongly reduced in Iceland and therefore further applications to data assimilation in a glacier model should rely on an algorithm that is robust to missing values.

Comparison within Situ Data
There is an overall good visual and statistical agreement between the MODIS shortwave broadband black-sky albedo and the in-situ albedo AWS measurements (Figures 4 and 5, Table 3).The seasonal variability in the in situ dataset is very high, ranging from 0.05 to 1 and is well captured by the MODIS product, as demonstrated by the high correlation coefficients (Table 3).When considering all the values, including low quality retrievals, the median of the correlation coefficients is 0.84, close to the value of 0.82 reported in Greenland [16].
The low consistency between the AWS and MODIS albedo at AWS hoffne and partly at AWS hoffef is probably partly due to the poor quality of the LI-COR radiometers used at those sites.Excluding AWS hoffne, there are still rather large differences during the summer at the lowest ablation stations breid, brune and langne.Field observations suggest that the differences at AWS brune and AWS breid might be explained with the resolution difference between the MODIS pixel and the AWS footprint (Section 3.3).When the winter accumulation has melted, the surrounding areas of AWS brune are characterized by dark volcanic fringes and at AWS breid by dirt bands and medial moraine for some years (Section 2.2).Apart from these explainable differences during the summer, a generally good agreement is obtained between the AWS and MODIS pixel albedo at breid, brune and langne (Table 3).The errors are not skewed, i.e., there is no clear tendency in the MODIS data to over or under-estimating the albedo.The mean absolute bias is 0.04 which is remarkably low given that the data span almost the entire possible range of albedo values (0-1).The RMSE by AWS ranges between 0.08 and 0.14 (median 0.12).These RMSE values are larger than that of previous studies (in Greenland RMSE is 0.07 and bias 0.022 [16]), presumably because only the high-quality retrievals were considered, whereas we only discarded values when solar zenith angle was lower than 70 • .The exclusion of low quality retrievals leads to lower albedo difference between MODIS and AWS albedo (RMSE range from 0.05 to 0.13, with a median value of 0.09) but also to a strong reduction of the number of data in the study area.High quality retrievals represent only 7% of the data that were used for validation at the AWS, ranging from 0% to 18% of the total number of values depending on the AWS.
In the accumulation area of Vatnajökull at AWSs bruef and tunef, where the spatial variability of the albedo is smaller, there are fewer available observations to evaluate the albedo errors (Section 3.1).However, the data obtained close to the ELA, at AWSs langef and brumi show a good correspondence over the whole study period and at AWS tunmi during most of the years, including the albedo drop following the eruptions in 2010 and 2011 and the albedo recovery in 2012 (Figures 4 and 5).Field observations suggest that the large discrepancy at AWS tunmi in 2001 to 2003 is partly explained by narrow water channels formed below the Kipp and Zonen radiometer before the melting exposes the previous year's summer surfaces, as well as localized sand and dirt exposed after melting through the winter accumulation (Figure 1).

Subpixel Variability
Some combinations of Landsat scenes and AWS site (Table 2) returned no-data in the computation of the standard deviation of variability within a MODIS sized pixel because all the pixels within the extraction mask were saturated for at least one of the bands (band 2 or 4) used to compute the broadband albedo (Figure 6).The number of values that were used to compute the standard deviation per AWS site ranged between 3 and 10.The RMSE of the difference between the MODIS and in-situ albedo (Section 3.2) was plotted against the mean of the standard deviations (STD) obtained from the extraction of the Landsat albedo (Figure 7).The results do not show a clear relationship (Spearman correlation is 0.48, p-value is 0.17), but the plot suggests that the RMSE of the difference between the in situ and MODIS data tends to increase when the corresponding Landsat subpixel spatial variability is higher.The lowest STD values are consistently obtained where the surface is less heterogeneous (Section 2.2; Figure 6), i.e., in the accumulation areas (AWSs bruef and tunef ) and at the stations hoffef and langef close to the current ELA, where surface slope does not allow water retention and slush.network is however too narrow and randomly spread to induce significant subpixel variability at the 30-m Landsat resolution.Hence, we believe that surface heterogeneity at a small scale not resolved by Landsat explains why high albedo RMSE are observed at tunmi and brumi despite the low STD values (Figure 7).On the other hand, the stratified dirt bands at AWS breid, the debris-covered terminus region at AWSs langne and hoffne, and the dark volcanic fringes at AWS brune (Section 2.2) are large enough features to be reflected in the high subpixel variability in the Landsat images (Figure 7).

Impact of Volcanic Eruptions
The evolution of the mean summer albedo shows that all the Icelandic ice caps except Drangajökull experienced significant albedo drops in the aftermath of the 2010 Eyjafjallajökull and 2011 Grímsvötn eruptions (Figure 8).Excluding Drangajökull, observed albedo decreases range from 0.12 to 0.27 with respect to the 2009 baseline (Table 4).Local reductions reaching 0.6 in absolute albedo value can be observed in the anomaly maps (Figure 8).The largest decrease is found in Eyjafjallajökull after the 2010 eruption (−0.27).In addition, a reduction of the spread after the eruptions is evident in every ice cap except Drangajökull, which means a reduction of the spatial variability in the summer albedo composite (Figure 8).This further indicates a widespread tephra deposition in the ablation and accumulation areas.This effect is also particularly marked in Eyjafjallajökull after the 2010 eruption.
Drangajökull was not affected by the 2011 eruptions and only marginally by the 2010 eruptions due to its remote location in northern Iceland (Table 4, Figure 9).MODIS data also indicate that the albedo decrease was less marked on the Langjökull ice cap in 2011 (Table 4).This is consistent with the fact that (i) the particles of the Grímsvötn tephra cloud were coarser and therefore heavier and did not get transported as far as the lighter Eyjafjallajökull tephra, and (ii) the dominant wind during the Grímsvötn eruption was blowing from the north (with a short period of SSW wind), leading to higher tephra deposition south of the eruption site [11].Hence, only small amounts of tephra were brought westwards to Langjökull.Deposited tephras were however redistributed during stormy dry periods in the summer of 2011 and reached Langjökull.In contrast, the albedo drop in Vatnajökull was larger in 2011 than in 2010 (Table 4), in particular in the central part of the ice cap (Figure 9).The flat surroundings at AWS tunmi and brumi (close to or at the current ELA) can facilitate narrow water channels or temporary slush below the instruments (Figures 1 and 7), explaining the high RMSE between the MODIS and in situ albedo values at those stations.Such surface drainage network is however too narrow and randomly spread to induce significant subpixel variability at the 30-m Landsat resolution.Hence, we believe that surface heterogeneity at a small scale not resolved by Landsat explains why high albedo RMSE are observed at tunmi and brumi despite the low STD values (Figure 7).On the other hand, the stratified dirt bands at AWS breid, the debris-covered terminus region at AWSs langne and hoffne, and the dark volcanic fringes at AWS brune (Section 2.2) are large enough features to be reflected in the high subpixel variability in the Landsat images (Figure 7).

Impact of Volcanic Eruptions
The evolution of the mean summer albedo shows that all the Icelandic ice caps except Drangajökull experienced significant albedo drops in the aftermath of the 2010 Eyjafjallajökull and 2011 Grímsvötn eruptions (Figure 8).Excluding Drangajökull, observed albedo decreases range from 0.12 to 0.27 with respect to the 2009 baseline (Table 4).Local reductions reaching 0.6 in absolute albedo value can be observed in the anomaly maps (Figure 8).The largest decrease is found in Eyjafjallajökull after the 2010 eruption (−0.27).In addition, a reduction of the spread after the eruptions is evident in every ice cap except Drangajökull, which means a reduction of the spatial variability in the summer albedo composite (Figure 8).This further indicates a widespread tephra deposition in the ablation and accumulation areas.This effect is also particularly marked in Eyjafjallajökull after the 2010 eruption.dust blown from the areas around the glaciers [35].
In 2012 the summer albedo returned to the pre-eruption value in Hofsjökull and Langjökull but not in the other ice caps.The persistence of a low albedo might be caused by the airborne redistribution of the volcanic dust in 2012 as already shown in the case of Vatnajökull [35].If so, wind transported dust may not have reached Hofsjökull and Langjökull, which are the farthest from the dust sources.Drangajökull was not affected by the 2011 eruptions and only marginally by the 2010 eruptions due to its remote location in northern Iceland (Table 4, Figure 9).MODIS data also indicate that the albedo decrease was less marked on the Langjökull ice cap in 2011 (Table 4).This is consistent with the fact that (i) the particles of the Grímsvötn tephra cloud were coarser and therefore heavier and did not get transported as far as the lighter Eyjafjallajökull tephra, and (ii) the dominant wind during the Grímsvötn eruption was blowing from the north (with a short period of SSW wind), leading to higher tephra deposition south of the eruption site [11].Hence, only small amounts of tephra were brought westwards to Langjökull.Deposited tephras were however redistributed during stormy dry periods in the summer of 2011 and reached Langjökull.In contrast, the albedo drop in Vatnajökull was larger in  4, Figure 9).

Conclusions
The large variability in the albedo of Icelandic ice caps is successfully captured with the low resolution MODIS albedo product MCD43A3, in spite of frequent cloud obstructions.We found a good agreement with in situ albedo data, even by including the retrievals flagged as low quality.A limitation of this comparison is the scale mismatch between a MODIS pixel and the AWS footprint.Field knowledge and Landsat data indeed suggest that there is a link between the variability of the albedo within a MODIS pixel and the discrepancy between AWS and MODIS data, however we did not identify a strong relationship.Higher resolution imagery might help resolve small scale  4, Figure 9).The albedo drop is larger in the accumulation areas where the "normal" summer albedo is typically higher than in the ablation areas due to the presence of snow or firn (Figure 9).In the ablation areas the albedo is already usually low due to the emergence of bare ice and the presence of light-absorbing impurities in the ice including old ash particles from past eruptions and mineral dust blown from the areas around the glaciers [35].
In 2012 the summer albedo returned to the pre-eruption value in Hofsjökull and Langjökull but not in the other ice caps.The persistence of a low albedo might be caused by the airborne redistribution of the volcanic dust in 2012 as already shown in the case of Vatnajökull [35].If so, wind transported dust may not have reached Hofsjökull and Langjökull, which are the farthest from the dust sources.

Conclusions
The large variability in the albedo of Icelandic ice caps is successfully captured with the low resolution MODIS albedo product MCD43A3, in spite of frequent cloud obstructions.We found a good agreement with in situ albedo data, even by including the retrievals flagged as low quality.A limitation of this comparison is the scale mismatch between a MODIS pixel and the AWS footprint.Field knowledge and Landsat data indeed suggest that there is a link between the variability of the albedo within a MODIS pixel and the discrepancy between AWS and MODIS data, however we did not identify a strong relationship.Higher resolution imagery might help resolve small scale variability around the station that is not resolved in Landsat 30-m data but that is known from the experience of the field visits.The AWS albedo measurements may also be affected by measurement biases due to the natural instability of the substratum on which the AWSs are placed, e.g., tilted surface slope below the downward pyranometer [36] or inaccurate horizontal leveling of the sensor.
The analysis of the MODIS albedo allowed us to characterize the large variability of the albedo in Icelandic glaciers over 2001-2012.For example the summer albedo decreased by up to 0.6 over large portions of the Icelandic ice caps after the two latest major eruptions (2010 and 2011) according to the MODIS product.Given the frequent occurrences of explosive eruptions in Iceland [37], it is recommended to consider the changed radiative forcing due to ash deposition (among other processes, like the insulating effect of the thickest tephra layers) to decipher the influence of climate change on glacier evolution in Iceland.Based on these encouraging results, our ongoing work will focus on the implementation of the MODIS albedo product within the HIRHAM5 regional climate model [38] to improve the computation of the surface energy and mass balance over the ice caps during the last decades.

Figure 1 .
Figure 1.Photographs of the Vatnajökull ice cap.(A) Volcanic tephra in the ablation area (photograph taken on 6 September 2014); (B) AWS tunmi in mid-summer 2001, showing water channels beneath the radiometer causing small scale variability in albedo.

Figure 1 .
Figure 1.Photographs of the Vatnajökull ice cap.(A) Volcanic tephra in the ablation area (photograph taken on 6 September 2014); (B) AWS tunmi in mid-summer 2001, showing water channels beneath the radiometer causing small scale variability in albedo.

Figure 3 .
Figure 3. Mean annual number of days with missing data during the period 2001-2012 in the MCD43A3 product.The black rectangle indicates the subset used in Section 3.4.The blue lines show the extent of the ice cap from the Corine Land Cover 2006 in Iceland.

Figure 3 .
Figure 3. Mean annual number of days with missing data during the period 2001-2012 in the MCD43A3 product.The black rectangle indicates the subset used in Section 3.4.The blue lines show the extent of the ice cap from the Corine Land Cover 2006 in Iceland.

Figure 4 .
Figure 4. Time series of daily in situ and MODIS albedo at AWS breid, brumi, hoffef and langef.Figure 4. Time series of daily in situ and MODIS albedo at AWS breid, brumi, hoffef and langef.

Figure 4 .
Figure 4. Time series of daily in situ and MODIS albedo at AWS breid, brumi, hoffef and langef.Figure 4. Time series of daily in situ and MODIS albedo at AWS breid, brumi, hoffef and langef.

Figure 6 .
Figure 6.An example of subpixel variability at each AWS site.Three summer Landsat scenes were used to compute albedo maps at 30 m resolution.A red square polygon corresponding to a MODIS pixel is drawn in each subset.Beside each subset is indicated the upper and lower values of the grayscale used to map the albedo.The beige areas are missing values either due to the Landsat-7 SLC-off issue or saturation (e.g., AWS breid).

Figure 6 .
Figure 6.An example of subpixel variability at each AWS site.Three summer Landsat scenes were used to compute albedo maps at 30 m resolution.A red square polygon corresponding to a MODIS pixel is drawn in each subset.Beside each subset is indicated the upper and lower values of the grayscale used to map the albedo.The beige areas are missing values either due to the Landsat-7 SLC-off issue or saturation (e.g., AWS breid).

Figure 7 .
Figure 7. Scatter plot of the RMSE obtained by the comparison of the MODIS albedo with the in situ data (Section 3.1) versus the mean of the standard deviations of the Landsat albedo that were sampled within the extent of a MODIS pixel.The mean STD was obtained by averaging the standard deviations computed for every combination of Landsat and AWS in Table 2.The color of the circles indicates if the AWS are in the accumulation area (red) or in the ablation area (blue).The photograph shows the brumi AWS during field maintenance on 24 Sepember 2001.

Figure 7 .
Figure 7. Scatter plot of the RMSE obtained by the comparison of the MODIS albedo with the in situ data (Section 3.1) versus the mean of the standard deviations of the Landsat albedo that were sampled within the extent of a MODIS pixel.The mean STD was obtained by averaging the standard deviations computed for every combination of Landsat and AWS in Table 2.The color of the circles indicates if the AWS are in the accumulation area (red) or in the ablation area (blue).The photograph shows the brumi AWS during field maintenance on 24 Sepember 2001.

Figure 8 .
Figure 8. Evolution of the mean summer albedo over the main ice caps (July-August) from 2001 to 2012.Each curve is centered on the mean annual summer albedo, which was computed by spatially averaging the annual summer albedo composites within the corresponding ice cap area (Section 2.1).The spread corresponds to the standard deviation of the same pixel values and therefore indicate the spatial variability of the summer albedo composites.The stars indicate the years of the Eyjafjallajökull and Grímsvötn eruptions.

Figure 8 .
Figure 8. Evolution of the mean summer albedo over the main ice caps (July-August) from 2001 to 2012.Each curve is centered on the mean annual summer albedo, which was computed by spatially averaging the annual summer albedo composites within the corresponding ice cap area (Section 2.1).The spread corresponds to the standard deviation of the same pixel values and therefore indicate the spatial variability of the summer albedo composites.The stars indicate the years of the Eyjafjallajökull and Grímsvötn eruptions.

Figure 9 .
Figure 9. Albedo of the main ice caps before and after the Eyjafjallajökull and Grímsvötn eruptions, and associated albedo changes.(a) Mean summer albedo in 2009 before the eruptions; (b) Mean albedo in summer 2010 after the Eyjafjallajökull eruption; (c) Albedo changes expressed as the difference between (b) and (a); (d) Mean albedo in summer 2011 after the Grímsvötn eruption; (e) Albedo changes expressed as the difference between (d) and (a).The volcanoes erupted in May 2010 and 2011 and are marked with red triangles in the panels (b) and (d).Here we excluded Drangajökull, which was not significantly affected by the eruptions (Table4, Figure9).

Figure 9 .
Figure 9. Albedo of the main ice caps before and after the Eyjafjallajökull and Grímsvötn eruptions, and associated albedo changes.(a) Mean summer albedo in 2009 before the eruptions; (b) Mean albedo in summer 2010 after the Eyjafjallajökull eruption; (c) Albedo changes expressed as the difference between (b) and (a); (d) Mean albedo in summer 2011 after the Grímsvötn eruption; (e) Albedo changes expressed as the difference between (d) and (a).The volcanoes erupted in May 2010 and 2011 and are marked with red triangles in the panels (b) and (d).Here we excluded Drangajökull, which was not significantly affected by the eruptions (Table4, Figure9).

Table 1 .
Average locations of the 10 AWSs (Figure2), time of operation and list of instruments.The Zone column indicates whether the AWS is located in the accumulation area (Ac) or in the ablation area (Ab), with ELA in parenthesis further indicating if it is close to the average equilibrium line altitude of the years 2001-2012.

Table 2 .
List of the Landsat scenes that were selected for this study.The AWS column indicates at which AWS location the extraction was performed.

Table 3 .
Statistics of the comparison between the MODIS shortwave broadband black-sky albedo and the in situ albedo AWS measurements.N is the sample size (number of matching dates between both datasets), RMSE is the root mean squared error, R is the Pearson's correlation coefficient, and B is the bias (mean difference).All correlations are significant at 1% level.

Table 4 .
Mean summer albedo changes for the main ice caps in Iceland after the volcanic eruptions of 2010 and 2011 with respect to 2009.

Table 4 .
Mean summer albedo changes for the main ice caps in Iceland after the volcanic eruptions of 2010 and 2011 with respect to 2009.