Dark Glacier Surface of Greenland’s Largest Floating Tongue Governed by High Local Deposition of Dust

: Surface melt, driven by atmospheric temperatures and albedo, is a strong contribution of mass loss of the Greenland Ice Sheet. In the past, black carbon, algae and other light-absorbing impurities were suggested to govern albedo in Greenland’s ablation zone. Here we combine optical (MODIS/Sentinel-2) and radar (Sentinel-1) remote sensing data with airborne radar and laser scanner data, and engage ﬁrn modelling to identify the governing factors leading to dark glacier surfaces in Northeast Greenland. After the drainage of supraglacial lakes, the former lake ground is a clean surface represented by a high reﬂectance in Sentinel-2 data and aerial photography. These bright spots move with the ice ﬂow and darken by more than 20% over only two years. In contrast, sites further inland do not exhibit this effect. This ﬁnding suggests that local deposition of dust, rather than black carbon or cryoconite formation, is the governing factor of albedo of fast-moving outlet glaciers. This is in agreement with a previous ﬁeld study in the area which ﬁnds the mineralogical composition and grain size of the dust comparable with that of the surrounding soils. Our analysis shows a decline in penetration depth with the transition into the impurity covered area. In the accumulation zone the stratigraphy can be resolved down to 13 m, with a decline to 1.5 m at the onset of the ablation zone. Further downstream no reﬂections of the last summer were found at all.


Introduction
Surface melt of the Greenland Ice Sheet contributes to about 60% of the mass loss of the ice sheet in the past decade [1,2]. Strong surface melt appears naturally in the lower elevated zones around the margins of the ice sheet, that are prone to higher summer air temperatures and lower albedo. Albedo itself is driven by light absorbing impurities and grain size. This effect has been important for the paleo-evolution of the Greenland Ice Sheet [3] and will remain a governing factor in future [4].
The surface of the ice sheet may consist of several constituents beside snow, firn or bare ice of different grain size: empty or water filled crevasses and thin fractures [5], surface melt water lakes [6], rivers [7], saturated patches [8] and impurities. Impurities may range from insoluble light absorbing impurities (LIA) as black carbon (BC, soot), dust arising from local sources or remote deserts, algae, Direct BC deposition measurements and documented behaviour in the snow pack are still limited despite previous efforts to include microphysical properties and chemistry related to aging in order to provide a successful simulation of BC concentrations deposited over the Arctic [25,27]. In addition, the few available measurements of the dry deposition of submicron particles over snow display a large variability (0.02 to 0.33 cm s −1 [28,29]). For areas with a high ratio of dust to BC, the scavenging effect will enforce an increasingly dominance of dust over BC [30]. This is what we presume to take place in our study area.
Nearby our study site, a field experiment aiming at the characterization of ice types, albedo and impurity content was conducted in 2006 [31]. The team of Bøggild et al. [31] found, that the mineralogical composition of most impurities matches that of the surrounding soils, suggesting that most of the material is derived from local sources. If equally distributed, these impurities are found to reduce the albedo of Holocene ice to 0.2-0.4. Interestingly, ice with more impurities per unit area, but aggregated in cryoconite holes, reaches an albedo of 0.6 as the impurities are shaded by the walls of the hole.
Shimada et al. [13] investigated a prominent dark area of the Greenland Ice Sheet by means of optical remote sensing using the MODIS sensors. They introduced a distinction between so called bare ice and dark ice using the reflectance at different wavelengths. Bare ice is defined at λ = 860 nm (band 2) with a reflectance of less than 0.6, while dark is defined at λ = 660 nm (band 1) with a reflectance of less than 0.4. Shimada et al. [13] find a maximum of 1.5% of the ice sheet to consist of dark ice. In the following we are referring to this classification, using Sentinel-2 bands of similar wavelength.
It has also been found that wind-blown sedimentary rock particles that form a layer of debris are the dominant trigger for surface melt ponds on the George VI Ice Shelf in Antarctica [32]. The setting of George VI Ice Shelf with the nearby mountains of Alexander Land is similar to our study area at the margin of the North-East Greenland Ice Stream (NEGIS) and in particular the lower range of 79 • Glacier (79 NG).
In this study, we take advantage of changes of drainage of supraglacial lakes leaving a bright surface and measure changes in reflectance over few years and set this into the context of surface properties along a fast flowing outlet glacier and a nearby slow moving area with exposed ice age ice.

Study Site
In this study we focus on the outlet region of the North East Greenland Ice Stream (NEGIS, Figure 1). This region is of high importance as 11% of the Greenland Ice Sheet is drained through the two main outlet glaciers Nioghalvfjerdsbrae (79 NG in the following) and Zachariae Isstrøm [33]. Next to these fast flowing outlet glaciers which reach maximum velocities of >2 km a −1 moderate velocities are prevailing allowing older ice to reach the surface of the ice sheet. The study region is further characterized by many supraglacial lakes which fill and drain frequently [5]. For a detailed analysis, we have selected three flow lines that encompasses the accumulation to the ablation zone as shown in Figure 1. (i) a profile along a slow moving area that exhibits Holocene and Pleistocene ice at the surface allowing for a comparison with Bøggild et al. [31], (ii) a profile along the northern branch of 79 NG, following an airborne survey crossing supraglacial lakes and (iii) a profile along the southern branch of 79 NG including a supraglacial lake that has annually drained for some time, allowing to assess the change in reflectance due to in-situ deposition of impurities. In the following we address these profiles with the names profile "Old Ice", "Northern Branch" and "Southern Branch" respectively. Figure 1 shows the study area including the locations of all three profiles.

Optical Remote Sensing: MODIS
The evolution of the snow cover on the glaciated and not-glaciated areas has been described using high-time resolved optical data based on the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard on the Terra & Aqua platforms. While the M.D10A.006 product can provide the Terra/Aqua Snow Cover Daily Global 500 m [34,35], the M.D09GA.006 can support the estimation of the Terra/Aqua Surface Reflectance Daily Global 1 km and 500 m [36,37]. The first product is characterized by the information concerning the cloud cover for each pixel, useful for performing cloud masking, and the estimation of the Normalized Difference Snow Index (NDSI), suitable for discriminating snow-covered pixels from other types of surfaces (especially soil and rocks). Since the two platforms provide two overpass per day, it is possible to have a daily revisit of the study site even if the spatial resolution is coarser than other platforms. The second product provides the daily spectral reflectance, atmospherically corrected, for 7 bands ranging from 400 nm to 2150 nm. We considered for each pixel a time-series of 246 layers (two layers per day) having 123 days every melting season: from the Julian day 122 (1st May) to 246 (31st August). All of the cloud-covered pixels can be discarded and a daily check between the two available platforms can confirm the final estimation both in terms of NDSI and spectral reflectance at 500 m resolution. The analysis has been focused on characterizing firstly the different melting seasons, from 2016 to 2018, looking at the Snow Cover Area (SCA) on bare soil/rock areas. Similarly the analysis has been concerned on estimating the occurrence of snow cover on glaciated areas, estimating the first day without snow, the "First Snow-Free Day" (FSFD) and the lowest reflectance value at 660 nm (band 1). The first discrimination, between snow/firn and ice surfaces, has been approached following the threshold proposed by Shimada et al. [13], based on reflectance at 860 nm (band 2) below 0.6. The estimation of the lowest value at 660 nm has been based on defining for each pixel the time window where snow cannot be detected and calculating the reflectance median in that time interval. The importance of λ 660 derives from Shimada et al. [13] that proposed a threshold of 0.4 for discriminating between "clean" ice and "dirty" ice.

Optical Remote Sensing: Sentinel-2
The detailed analysis with high-resolution data has been approached using Sentinel-2 images. This satellite constellation is composed by two platforms that are active since 2015 and 2017 in the framework of the European Copernicus EO Programme [38]. The Copernicus Open Access Hub (https://scihub.copernicus.eu/) provides L1C data that must be converted from top-of-atmosphere to bottom-of-atmosphere reflectance. The analysis has been developed selecting cloud-free images and processing those data with the Sen2Cor toolbox [39]. The atmospheric correction process required some input data that we retrieved from the AERONET network Giles et al. [40]. Considering the available values of the 500 nm aerosol optical depth (<0.1) at Thule and Ittoqqortoormiit (since there are no available station in the study area), we used as input data a visibility of 40 km. The total columnar ozone content was selected referencing to Bahramvash Shams et al. [41] and we used as input values 400 and 300 DU, respectively for April 2017 and July 2017, The aerosol type was selected to be maritime and the atmosphere type as mid latitude winter, following the indications of König M and N [42]. The preparation of a L2A final product, with 10 m resolution and surface reflectance, has provided information about the spectral optical behaviour (12 bands with wavelengths ranging from 440 nm to 2200 nm) of the surface at selected dates characterized by a low cloud coverage. The Sentinel-2 bands are distributed as follows: B1 at 443 nm, B2 at 490 nm, B3 at 560 nm, B4 at 1665 nm, B5 at 705 nm, B6 at 740 nm, B7 at 783 nm, B8 at 842 nm, B9 at 945 nm, B11 at 1614 nm, B12 at 2202 nm. We defined for the next analysis two normalized difference indexes useful for enhancing the different contribution of visible, near-infrared and short-wave infrared wavelength domains.

Radar Remote Sensing: Sentinel-1 Polarimetry
Dual-polarisation SAR radar imagery can be used to distinguish between different types of reflection. A flat surface (as, e.g., a water body with no significant waves) causes a specular reflection, which means that the side-looking radar system hardly receives any return power at all. On a rough surface the reflection is dominated by Bragg-like scattering [43] and the incident horizontal polarisation is mainly reflected back as a horizontal polarised wave (HH). In contrast, volume scattering (e.g., in the firn pack) causes many random reflections within the volume and a part of the incident horizontally polarised wave is reflected in vertical polarisation (HV). The change in polarisation can therefore be used to distinguish between the scattering mechanism, that can subsequently be used for inferring material properties. For this study we are only using Sentinel-1 scenes in HH and HV polarisation for two different time periods, one representing frozen conditions and one representing summer, melt conditions. The Single Look Complex (SLC) S1 scenes (IW mode) are processed using GAMMA Software [44] with the following steps: (1) calibration for instrumental noise, (2) mosaicing, (3) multilooking for speckle reduction and (4) georeferencing the scenes to EPSG:3413. The resulting products has a ground pixel resolution of 50 × 50 m. Topographic effects in the radar backscatter are corrected with a line-of-sight area normalization that transfers σ 0 to γ 0 using the ArcticDEM [45]. The details of polarimetric processing is presented in [5].
We selected two periods for S1 acquisition, one close to the acquisition of our airborne data in April 2018 and one during the melt season.

Roughness
To estimate surface roughness along the profiles we analyzed laser scanner data conducted with a RIEGL-LMSQ560 during the airborne survey in April 2018. The geo-located laser scanner point cloud (Latitude, Longitude, Elevation referenced to WGS84) was obtained using laser range and angle measurements, post processed GNSS trajectory (PPP solution using Waypoint 8.4) and attitude information (Inertial Measuring Unit: Honeywell Laseref V) within our processing scheme, following Hofton et al. [46].
Surface roughness, here defined as standard deviation of a high pass filtered laser scanner elevation model (ALS-DEM), was estimated at 1 km along track spacing. A series of 20% overlapping ALS-DEMs, each with a pixel resolution of 1 m and dimension of 1200 m × 1200 m, were conducted from the geo-located laser scanner point cloud data using an inverse distance weighting interpolation scheme. The final ALS-DEMs were referenced to the EPSG:3413 coordinate system.
High pass filtering (subtracting a smoothed version of the ALS-DEM using a 300 m kernel) was applied to remove large scale topographic features, as we were only interested in small scale roughness in the order of meters (e.g., sastrugis).

Firn Modelling
To simulate the density and age of the firn at the ice sheet surface, we relate spatial points along various flow lines, including the three selected sites, to points in time. This is done using the present time surface velocity field. For every flowline a one dimensional firn densification model is forced over time with corresponding mean annual values of the surface temperature and accumulation rate. The data is obtained from the regional climate model RACMO [47].
The model follows an approach of continuum mechanics solving conservation of mass, linear momentum and energy. The needed constitutive equations are derived from sintering mechanics, as the densification of firn can be compared to the process of hot isostatic pressing, a common process in sintering. This physics based approach of firn densification modelling provides the possibility to simulate not only the density of the firn profile but also its age, temperature, grain size and various other properties [48].
The implementation of the model allows to simulate varying accumulation and ablation rates at high temporal resolution. For accumulation a constant density of ρ 0 = 360 kgm −3 is chosen to represent recent deposited firn. In case of a negative accumulation rate the ablated mass is removed from the top of the simulated firn profile, uncovering older and more dense firn.

Annual Deposition of Impurities
At many locations of 79 NG we find the glacier surface to be clean after the drainage of supraglacial lakes. Some examples are shown in Figure 2. For quantifying the change in reflectance over time, we found unfortunately only one appropriate lake, due to constraints such as S2 launch, cloud free coverage and snow free coverage with frozen surface. This lake is shown in Figure 3. We named this lake 'Lake iCUPE'. Lake iCUPE has drained annually in the past years on: 21 August 2015, 17 July 2016, 16 July 2017, 1 August 2018, 17 July 2019 and 22 June 2020. After the drainage the surface appears bright, while over time this bright spots darkens, as can be seen in the right panel of Figure 3. As in 2019 major cracks further south appeared, which reached into our lake ground, the undisturbed time period ends 2018. In 2015 there has been no S2 coverage, thus we end up with the time period 2016-2018 to quantify the change in reflectance.
From a Sentinel-2 scene acquired on 28 July 2017 three brighter spots were digitized manually and the timing of their center positions was quantified from a flowline analysis (Figure 3). By projecting the locations of the bright spots back in time we found the centers of the two former lake grounds to be at the actual lake position on 11 August 2015 and 24 July 2016. As the glacier moves with a speed of about 1000 m/a, the lake ground (∼800-900 m length) of year 2 is separated from the year 1's lake ground surface by about 100-150 m. A comparison of the reflectance of the lake ground over four subsequent patches of former lake ground indicates the in-situ darkening by deposition: the reflectance drops within the first year after drainage from 0.8 to 0.68 and further to 0.62 two years after drainage (Figure 3). After that, the former lake ground cannot be distinguished from surrounding areas. Here we assume that snow albedo is connected to reflectance as analogous but we are only making qualitative use of the terms.   We infer that the change in reflectance is due to local deposition of dust, as we can rule out lag-deposit to govern this effect: a comparison with a further upstream situated lake, Lake CE (Figure 1), shows that in that area, changes of reflectance with distance (time) is not significant (see Figure S1). Ablation rates and flow speed together amounts to roughly the same amount of melt at the surface, leaving to similar amounts of lag-deposit. This is in conflict to the strong change in reflectance at the floating tongue (Lake iCUPE), that can only be achieved by deposition in addition to lag-deposits.
We hypothesize that this deposition is caused by surrounding soil and mountains (despite that we have not analysed any samples), as Bøggild et al. [31] have shown dust to come indeed from nearby mountains. Although Lake iCUPE is somewhat away from the former field site of the study of [31], it is very likely that at about the same elevation the source of the dust is similar. Fresh wet and dry deposition from non-local sources is anticipated to have orders of magnitude lower impact and is thus unlikely to cause such severe change in reflectance. Ref. [18] have concluded that annual deposited BC was causing a marginal reduction of the snow albedo, while post-depositional processes may cause a stronger direct and indirect impact on the snow albedo via enhanced metamorphism processes but in any case it has several order of magnitude lower impact than Holocene ice.
In addition to that, grain-size increase due to fluid sintering is similarly unlikely to be the cause, as the time scale for grain-growth is by far larger, consequently such a change in reflectance due to change in specific surface area would take considerably more time than just one year. Having ruled out these two causes only a local source is plausible. In order to check if this is changing with altitude, we analyse a second lake location further upstream in an area where lakes drain, too.
In order to broaden this, we are analysing profiles along flow lines of the glacier, that are shown in Figure 1. One profile is along the southern branch, crossing Lake iCUPE, another profile in the northern branch of 79 NG, and a third one is ending at the field site of the study of Bøggild et al. [31]. We start with the first, which we name ?Southern Branch?, picking up the findings of Bøggild et al. [31], who found Holocene and Pleistocene ice at the glacier surface.

Profile "Southern Branch"
The last along flow profile addressed as "Southern Branch" and situated furthest south (see Figure 1) is dedicated to cross "Lake iCUPE". This supraglacial lake is used to study the effect of in-situ deposition of impurities in detail as described in the following section. Figure 4 is showing a similar figure than for profiles "Old Ice" and "Northern Branch". Also there the transition into an impermeable surface layer (pore-close off density is ∼830 kg/m 3 ) is at an elevation of about 900 to 1000 m a.s.l. The backscatter shows more variation in the upstream area in this profile. Again the crevassed area is coincident with higher backscatter, both in summer and winter conditions. Similar to profile "Northern Branch", the downstream part shows significant lower reflectance values, indicating a surface with high impurity content ( Figure 3).
The average for the areas with elevation above 1000 m a.s.l. (distance from the ice margin > 20 km) showed on April about 0.99, 1.00 and 0.066 for bands 3, 8 and 11, respectively. The values in areas closer to the ice margin (<1000 m a.s.l.) were 0.95, 0.96 and 0.10 for the above mentioned bands and same season. The indexes were also in this case close to each other highlighting the occurrence of similar surface types (Index1: −0.017 and −0.004: Index2: 0.88 and 0.81). The same analysis evidences on July reflectance values significantly different: 0.94, 0.78 and 0.034 far from the ice margin; 0.78, 0.60 and 0.031 close to the ice margin. The Vis-NIR index (Index 1) highlighted limited differences (0.10 against 0.13) as well as the NIR-SWIR index (Index 2) (0.92 against 0.91).
In addition to these profiles we were using AWI's snow radar (UWB-M) to detect reflections in the upper layers. Normally, summer layers are forming a stratigraphy that is also often used to determine accumulation rates over past decades to centuries [49]. Figure S4 displays radargrams along some short transects. Our analysis shows a decline in penetration depth with the transition into the impurity covered area. In the accumulation zone the stratigraphy can be resolved down to 13 m, with a decline to 1.5 m at the onset of the ablation zone. Further downstream no reflections of the last summer were found at all.

Profile "Northern Branch"
The data along profile "Northern Branch", representing the fast moving part of 79 NG (shown in Figure 1), is presented in Figure 5. Firn modelling indicates a transition of a permeable to impermeable surface layer at around 120 km in an upstream located topographic plateau at about 1000 m a.s.l. (see also Figures S2 and S3). The profile crosses supraglacial lakes (grey bars in Figure 5), over which HH and HV backscatter drops, both in summer and winter conditions [5]. In winter the difference of HH and HV backscatter remains nearly constant down to the crevassed area, while in summer the difference increases towards lower elevations, indicating a decreasing volume scattering due to wet surface conditions. Starting from around 140 km surface roughness increases and becomes the dominating effect on radar backscatter and optical reflectance. With the onset of high roughness more of the radar signal is reflected back to the sensor. Moreover, a part of this increase in backscatter can also be explained by side reflections from sharp crevasses that are dominant in this area.
Most striking is the difference in backscatter between upstream and downstream areas indicating differences in surface characteristics: In the upstream (i.e., higher elevated) region, the backscatter is dominated by the ratio firn pack penetration while further downstream the crevassed surface roughness dominates the signal. Upstream in winter, the backscatter is nearly constant for elevations above 1000 m and decreases quite slowly in lower elevations. In contrast in summer, the backscatter drops very fast around 1500 m a.s.l., indicating a significantly wetter surface below.
Looking at reflectance values the average for the areas with elevation above 1000 m a.s.l. (distance from the ice margin > 20 km) showed on April about 0.97, 1.00 and 0.073 for bands 3,   Figure 6 presents various data sets along the flow line of profile "Old Ice", located in a slow moving area outside the main trunk of 79 NG (see also Figure 1). As the flowline is heading towards the area studied by Bøggild et al. [31] we mark in Figure 1) their study site, allowing for a direct comparison. Visual inspection of neither optical nor radar satellite imagery indicate any heavily crevassed region in this profile, nor lakes, but numerous rivers. Our firn modelling results indicate that the ablation zone is consisting predominantly of an ice surface with a density above pore-close off (see also Figures S2 and S3), blocking surface melt water to penetrate into the glacier, fostering surface run-off. The steep slope in this area leads is further intensifying run-off through rivers. In-situ deposition of impurities from nearby soil, as observed by [31] is likely to either remain on the surface or to be washed off with the run-off. Sentinel-2 data indicates a change from clean late Holocene ice to early Holocene ice and also the transition to Pleistocene ice, as also reported by Bøggild et al. [31]. With the onset of dirty Holocene ice, both, winter and summer periods show a smooth decline in HH and HV backscatter. Also the transition in reflectance of various bands in S2 is smooth, with a drop of about 0.1 in the visible range over the last 60 km towards the ice margin. In details, referring to areas in Figure 6, the average reflectance value for the areas with elevation above 1000 m a.s.l.  Profile "Old Ice" with nearly same panels as Figure 5, except that no laser scanner based roughness is available here. The glacier flow direction is from left (upstream) to right (downstream). The grey bars indicate here Holocene ice (wide bar) and Pleistocene ice (narrow bar), that is detected in S2 imagery following [31].

Snow, Bare and Dark Ice Surface
The evolution of the snow cover can only be assessed with remote sensing data with a short repetition time, as provided by the MODIS platform. While this platform provides lower spatial resolution, it ensures higher time resolution that can be used to obtain a time series of spectral reflectance during each melting season. This approach was followed by Shimada et al. [13] with the aim of detecting bare ice and dark ice on the surface. We performed such an analysis in the 79 NG area for three consecutive years 2016, 2017 and 2018, focusing particularly on the relationship with the mountains surrounding the floating tongue of the glacier. The reduction of the fractional snow cover occurring on the outcropping rocks showed (Figure 7 Figure 1a. The dashed line represents the reflectance threshold defined by [13] for discriminating between snow cover and ice. All the considered data were obtained analysing MODIS data as described in Section 3.1.

Discussion
We find that the backscatter in the two profiles in the main trunk of the glacier has its lowest values in the low elevated part in which also the optical reflectance is low. In addition to that, the drop of backscatter in wet, summer, conditions is in the order of 10 dB, whereas the drop in the areas with low reflectance is only about 5 dB (Figures 4 and 5). This means that the water content does not affect the absorption in the same way in areas of high or low impurity content at the surface. This needs to be interpreted in relation to the change in backscatter from low to high impurity content of about 10 dB. Once the impurities at the surface reach a significant level, the backscatter is more dominated by impurity content than by moisture. This is further supported by the 2 m air temperature (see Figure S5), which has the same range at the lower elevated floating tongue than in the upper areas (km 70-140 in Figures 4 and 5).
Non local sources annual rate of change is currently stabilized or only slightly increasing compared to past decades. For long term trends in BC, global emissions were estimated to grow by about 15% by 2010 compared to 1990 [50]. The aerosol light absorption levels do not show a significant trend in the Arctic as recently displayed in the analysis by Collaud Coen et al. [22] for long term data from Arctic stations Alert, Summit and Zeppelin-Ny Ålesund.
How will this develop in future? With rising temperatures the area that may undergo melt in future is likely to increase. Given that the floating tongue is experiencing now already strong melt during summer, we consider two scenarios: (i) with increasing melt the number of ponds on the floating tongue is increasing and/or the size of the ponds is increasing. In both cases the albedo becomes closer to that of open ocean. (ii) Increasing melt will lead to more surface water streams that flush the surface around the melt and sweep more dirt into the lakes. The consequence may be that the ice surface is coming to an increased albedo that is predominantly driven by coarsening of the upper layer. Also the area upstream the floating tongue that is currently prone to in-situ deposition from nearby mountains may undergo change. With increasing melt the cycle of lake filling/drainage can increase and lake size may increase, as well as rivers may increase in size and number. All that would lead to a more regular cleaning of the surface and consequently increase in albedo.
Taking the findings of Willeit and Ganopolski [3] into account, we recognise that there is a need of representation of the scavenging effect in firn hydrology modelling, which itself requires an evolution equation of grain size in wet conditions. Similar to that, using an ice sheet model tracer for dust concentration for accounting for albedo of early Holocene and Pleistocene ice in the ablations zones may become an important factor for projections of Greenland's contribution to sea level rise.
Satellite remote sensing approaches to detect albedo in ablation zones may provide valuable information for assessing the influence over time span of decades and to validate simulated albedo. Approaches like the one by Kokhanovsky et al. [51] that employ a sophisticated atmospheric correction to Sentinel-3's OLCI instrument obtaining an albedo estimation of high accuracy can potentially be combined with a debris cover ratio, as proposed by Azzoni et al. [52] allowing for separating in future better between debris/dust influence and grain size.

Conclusions
Melt in the ablation zone of fast moving outlet glaciers in Greenland that are situated in vicinity to surrounding mountains is an important factor in ice sheet mass loss. We find, that within two years time, the reflectance drops by about 25%, turning a clean ice surface after supraglacial lake drainage into a dirty surface. A change of this order of magnitude in such short time scale is unlikely to be facilitated by grain growth, black carbon or cryoconite formation, but points towards high local dust deposition. Given the magnitude of this effect we suggest that dust deposition and scavenging should be taken into account for simulation of albedo in regional climate models, despite being challenge to model adequately.