The darkening of the Greenland ice sheet: trends, drivers, and projections (1981–2100)

The surface energy balance and meltwater production of the Greenland ice sheet (GrIS) are modulated by snow and ice albedo through the amount of absorbed solar radiation. Here we show, using space-borne multispectral data collected during the 3 decades from 1981 to 2012, that summertime surface albedo over the GrIS decreased at a statistically significant (99 %) rate of 0.02 decade between 1996 and 2012. Over the same period, albedo modelled by the Modèle Atmosphérique Régionale (MAR) also shows a decrease, though at a lower rate (∼−0.01 decade) than that obtained from space-borne data. We suggest that the discrepancy between modelled and measured albedo trends can be explained by the absence in the model of processes associated with the presence of light-absorbing impurities. The negative trend in observed albedo is confined to the regions of the GrIS that undergo melting in summer, with the drysnow zone showing no trend. The period 1981–1996 also showed no statistically significant trend over the whole GrIS. Analysis of MAR outputs indicates that the observed albedo decrease is attributable to the combined effects of increased near-surface air temperatures, which enhanced melt and promoted growth in snow grain size and the expansion of bare ice areas, and to trends in light-absorbing impurities (LAI) on the snow and ice surfaces. Neither aerosol models nor in situ and remote sensing observations indicate increasing trends in LAI in the atmosphere over Greenland. Similarly, an analysis of the number of fires and BC emissions from fires points to the absence of trends for such quantities. This suggests that the apparent increase of LAI in snow and ice might be related to the exposure of a “dark band” of dirty ice and to increased consolidation of LAI at the surface with melt, not to increased aerosol deposition. Albedo projections through to the end of the century under different warming scenarios consistently point to continued darkening, with albedo anomalies averaged over the whole ice sheet lower by 0.08 in 2100 than in 2000, driven solely by a warming climate. Future darkening is likely underestimated because of known underestimates in modelled melting (as seen in hindcasts) and because the model albedo scheme does not currently include the effects of LAI, which have a positive feedback on albedo decline through increased melting, grain growth, and darkening.

and a record of 570 ± 100 Gt in total mass loss, doubling the average annual loss rate of 260 ± 100 Gt for the period 2003-2012 .
Net solar radiation is the most significant driver of summer surface melt over the GrIS (van den Broeke et al., 2011;Tedesco et al., 2011), and is determined by the combination of the amount of incoming solar radiation and surface albedo. Variations in snow albedo are driven principally by changes in snow grain size and by the presence of light-absorbing impurities (LAI, Warren and Wiscombe, 1980). Generally, snow albedo is highest immediately following new snowfall. In the normal course of destructive metamorphism, the snow grains become rounded, and large grains grow at the expense of small grains, so the average grain radius r increases with time (LaChapelle, 1969). Subsequently, warming and melt/freeze cycles catalyse grain growth, decreasing albedo mostly in the near-infrared (NIR) region (Warren, 1982). The absorbed solar radiation associated with this albedo reduction promotes additional grain growth, further reducing albedo, potentially accelerating melting. The presence of LAI such as soot (black carbon, BC), dust, organic matter, algae, and other biological material in snow or ice also reduces the albedo, mostly in the visible and ultraviolet regions (Warren, 1982). Such impurities are deposited through dry and wet deposition, and their mixing ratios are enhanced through snow water loss in sublimation and melting (Conway et al., 1996;. Besides grain growth and LAI, another cause of albedo reduction over the GrIS is the exposure of bare ice: once layers of snow or firn are removed through ablation, the exposure of the underlying bare ice will further reduce surface albedo, as does the presence of melt pools on the ice surface (e.g. Tedesco et al., 2011).
Most of the studies examining albedo over the whole GrIS have focused on data collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) starting in 2000 (e.g. Box et al., 2012;. At the same time, regional climate models (RCMs) have been employed to simulate the evolution and trends of surface quantities over the GrIS back to the 1960s using reanalysis data for forcing (e.g. Fettweis et al., 2013). Despite the increased complexity of models, and their inclusion of increasingly sophisticated physics parameterizations, RCMs still suffer from incomplete representation of processes that drive snow albedo changes, such as the spatial and temporal distribution of LAI, and from the absence of in situ grain size measurement to validate modelled snow grain size evolution. In this study, we first report the results from an analysis of summer albedo over the whole GrIS from satellite for the period 1980-2012, hence expanding the temporal coverage with respect to previous studies. Then, we combine the outputs of an RCM and in situ observations with the satellite albedo estimates to identify those processes responsible for the observed albedo trends. The model, Modèle Atmosphérique Régionale (MAR), is used to simulate surface temperature, grain size, exposed ice area, and surface albedo over Greenland at large spatial scales. MAR-simulated surface albedo is tested against surface albedo retrieved under the Global LAnd Surface Satellite (GLASS) project, and it is used to attribute trends in GLASS albedo. Lastly, we project the evolution of mean summer albedo over Greenland using the MAR model forced with the outputs of different Earth system models (ESMs) under different CO 2 scenarios. Discussion and conclusions follow the presentation of the methods and results.

The MAR regional climate model and its albedo scheme
Simulations of surface energy balance quantities over the GrIS are performed using the Modèle Atmosphérique Régionale (MAR; e.g. Fettweis et al., 2005Fettweis et al., , 2013. MAR is a modular atmospheric model that uses the sigma-vertical coordinate to simulate airflow over complex terrain and the Soil Ice Snow Vegetation Atmosphere Transfer scheme (SIS-VAT, e.g. De Ridder and Gallée, 1998) as the surface model. MAR outputs have been assessed over Greenland in several studies (e.g. Tedesco et al., 2011;Fettweis et al., 2005;Vernon et al., 2013;Rae et al., 2012;van Angelen et al., 2012), with recent work specifically focusing on assessing simulated albedo over Greenland (Alexander et al., 2014). A discussion of this evaluation is presented later in the manuscript. The snow model in MAR is the CROCUS model of Brun et al. (1992), which calculates albedo for snow and ice as a function of snow grain properties, which in turn are dependent on energy and mass fluxes within the snowpack. The model configuration used here has 25 terrainfollowing sigma layers between the Earth's surface and the 5 hPa model top. The spatial configuration of the model uses the 25 km horizontal resolution computational domain over Greenland described in Fettweis et al. (2005). The lateral and lower boundary conditions are prescribed from meteorological fields modelled by the global European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis (ERA-Interim, http://www.ecmwf.int/en/research/ climate-reanalysis/era-interim). Sea-surface temperature and sea-ice cover are also prescribed in the model using the same reanalysis data. The atmospheric model within MAR interacts with the CROCUS model, which provides the state of the snowpack and associated quantities (e.g. albedo, grain size). No nudging or interactive nesting was used in any of the experiments. The MAR albedo scheme is summarized below. Surface albedo is expressed as a function of the optical properties of snow, the presence of bare ice, whether snow is overlying ice (and whether the surface is waterlogged), and the presence of clouds. In the version used here (MARv 3.5.1), the broad-M. Tedesco et al.: The darkening of the Greenland ice sheet: trends, drivers, and projections (1981-2100) 479 band albedo (α s , 0.3-2.8 µm) of snow is a weighted average (Eq. 1) of the albedo in three spectral bands, α 1 , α 2 , and α 3 , which are functions of the optical diameter of snow grains (d, in metres), as modified from equations by Brun et al. (1992;e.g. Lefebre et al., 2003;Alexander et al., 2014): The optical diameter d is, in turn, a function of snow grain properties and it evolves as described in Brun et al. (1992). In MAR, the albedo of snow is calculated by Eqs.
(1)-(4), but it is not permitted to drop below 0.65. For the transition from snow to ice, MAR makes the albedo an explicit function of density. On a polar ice sheet, densification of snow/firn/ice occurs in three stages, with a different physical process responsible for the densification in each stage (Herron and Langway Jr., 1980;Arnaud et al., 2000). Newly fallen snow can have density in the range 50-200 kg m −3 . After then, densification can occur due to wind processes, which break and round grains, forming windslab of density typically around 300-400 kg m −3 . The remaining densification happens by grain-boundary sliding, attaining a maximum density of ∼ 550 kg m −3 at the surface. Old melting snow at the surface in late summer typically has this density, but does not exceed it, because this is the maximum density that can be attained by grain-boundary sliding and corresponds to the density of random-packing of spheres (Benson, 1962, p. 77). Further increases of density (the second stage) occur in firn under the weight of overlying snow, by grain deformation (pressure-sintering). In this case the density range is 550-830 kg m −3 . At a density of 830 kg m −3 the air becomes closed off into bubbles and the material is called ice. In the third stage, the density of ice increases from 830 to 917 kg m −3 by shrinkage of air bubbles under pressure. Moving down the slope along the surface of the GrIS, at the transition between the accumulation area and the ablation area, the snow melts away, exposing firn. Continuing farther down, the firn melts away, exposing ice. The albedo of firn may be approximated as a function of its density, ρ, interpolating between the minimum albedo of snow and the maximum albedo of ice. In MAR these values of albedo are set to 0.65 and 0.55, respectively. We would then have for the density range of firn (550-830 kg m −3 ) The MARv3.5.1 version used here maintains a minimum albedo of 0.65 for any density up to 830 kg m −3 , and specifies the gradual transition from snow albedo to ice albedo across the density range 830-920 kg m −3 . This means that the albedo of exposed firn is not allowed to drop below 0.65, with the result that the positive feedbacks of snow/firn/ice albedo will be muted in MAR. This aspect is being addressed in future versions of MAR (MAR v3.6) and a sensitivity analysis is being conducted to evaluate the impact of the changes on the albedo values when snow is transitioning from firn to ice. Such analysis is computationally expensive and preliminary outputs will be published once available. In MAR, the albedo for bare ice is a function of the accumulated surface meltwater preceding runoff and specified minimum (α i,min ) and maximum (α i,max ) bare ice values: Here α i,min and α i,max are set, respectively, to 0.4 and 0.55, K is a scale factor set to 200 kg m −2 , and M SW(t) is the time-dependent accumulated excess surface meltwater before runoff (in kg m −2 ). When a snowpack with depth less than 10 cm is overlying a layer with a density exceeding 830 kg m −3 (i.e. ice), the albedo in MAR is a weighted, vertically averaged value of snow albedo (α S ) and ice albedo (α I ; e.g. if snow depth is 3 cm then albedo is obtained by multiplying the snow albedo by 0.3 and adding the ice albedo multiplied by 0.7). When the snowpack depth exceeds 10 cm, the value is set to α S . The presence of clouds can increase snow albedo because they absorb at the same NIR wavelengths where snow also absorbs, skewing the incident solar spectrum to wavelengths for which snow has higher albedo ( Fig. 5 of Grenfell et al., 1981;Fig. 13 of Warren, 1982;Greuell and Konzelman, 1994), in which case the albedos of snow and ice are adjusted based on the cloud fraction modelled by MAR (Greuell and Konzelman, 1994).

The GLASS albedo product
The GLASS surface albedo product (http://glcf.umd.edu/ data/abd/) is derived from a combination of data collected by the Advanced Very High Resolution Radiometer (AVHRR) and the Moderate Resolution Imaging Spectroradiometer (MODIS, Liang et al., 2013). Shortwave broadband albedo (0.3-3 µm) is provided every 8 days at a spatial resolution of 0.05 • (∼ 56 km in latitude) for the period 1981-2012. GLASS albedo data with a resolution of 1 km are also available from 2000 to 2012 but these data are not used here for consistency with the data available before 2000. There have been several efforts to make the AVHRR and MODIS albedo products consistent within the GLASS product, including the use of the same surface albedo spectra to train the regression and the use of a temporal filter and climatological background data to fill data gaps . Monthly averaged broadband albedos from GLASS-AVHRR and GLASS-MODIS were cross-compared over Greenland for those months when there was overlap (July 2000(July , 2003(July , and 2004, revealing consistency in GLASS-retrieved albedo from the two sensors (He et al., 2013). More information on the GLASS data processing algorithm and product is available in Zhao et al. (2013) and Ying et al. (2014).
The GLASS product provides both black-sky albedo (i.e. albedo in the absence of a diffuse component of the incident radiation) and white-sky albedo (albedo in the absence of a direct component, with an isotropic diffuse component). The actual albedo is a value interpolated between these two according to the fraction of diffuse sunlight, which is a function of the aerosol optical depth (AOD) and cloud cover fraction. In the absence of the full information needed to properly reconstruct the actual albedo, here we use in our analysis the black-sky albedo, because we focus mostly on albedo retrieved under clear-sky conditions. Our analysis using the white-sky albedo (not shown here) is fully consistent with the results obtained using the black-sky albedo and reported in the following. A full description of the GLASS retrieval process and available products can be found in Liang et al. (2013) and references therein. An assessment of the GLASS product complementing existing studies is reported below.
Data collected by the MODIS Terra and Aqua sensors are used in the GLASS albedo retrieval for the period 2000(2000for Terra and 2002 for Aqua, respectively). Wang et al. (2012) have shown that the MODIS Terra sensor has been degrading at a pace that can be approximated by a second-order polynomial, with the coefficients being spectrally dependent. Over Greenland, the impact of sensor degradation on albedo trends has been estimated at −0.0059 decade −1 (Stroeve et al., 2013). Polashenski et al. (2015) found a much greater impact on retrieved broadband albedo from Terra sensor degradation (−0.03 decade −1 ). However, Polashenski et al. (2015) use a daily product (MOD10A1) rather than a 16-day integrated product as in the case of GLASS (e.g. Ying et al., 2014), which does account for the bidirectional reflectance distribution function (BRDF) at high solar zenith angles. The performance of the MODIS daily product has been shown to deteriorate with latitude (e.g. Alexander et al., 2015). On the other hand, the use of the BRDF (as in the case of the GLASS product) improves the performance of the product at high latitudes (Alexander et al., 2015). This, together with the good agreement between the MCD43 albedo product and the surface station albedo data (Alexander et al., 2014), gives us confidence in the GLASS trends.
We complement previous assessments of the MODIS and GLASS albedo, evaluating the absolute accuracy of the GLASS retrievals by comparing monthly GLASS albedo to in situ measurements of albedo collected at automatic weather stations of the Greenland climate network (GC-Net, Steffen and Box, 2001). GC-Net data are distributed at hourly temporal resolution and were temporally averaged to match the temporal window used in the GLASS product data. The root mean square error (RMSE), percentage RMSE (RM-SEp), and the slope of a linear fit between GLASS and in situ measured albedos for 12 stations are given in Table 1. The number of available years used for the statistics is also reported for each station. We considered only stations for which at least 10 years were available for the analysis in at least one of the months. Our results are consistent with the findings reported by Alexander et al. (2014) and Stroeve et al. (2005Stroeve et al. ( , 2013 concerning the assessment of the MODIS albedo products over the GrIS. The mean value of the RMSE for all stations is 0.04-0.05 in all months, with individual station values as high as 0.15 for station JAR1 in August and as low as 0.01 for Summit and Saddle stations in June. The relatively large RMSE value for JAR1 (and other stations located within the ablation zone) is probably due to the heterogeneity of albedo values within the pixel containing the location of the station and to the point-scale nature of the in situ observations. At Summit, where spatial inhomogeneity on the surface is small, it is reasonable to assume that the effect of the spatial scale and heterogeneity on the comparison is smaller.

Albedo trends
The time series of the mean summer GLASS albedo values between 1981 and 2012 over Greenland can be separated into two distinct periods (Fig. 1a): the period 1981-1996, when albedo shows no trend, and a second period, 1996-2012, when a statistically significant trend (99 %) is detected. The year 1996 was identified as yielding the highest value of the coefficient of determination when fitting the albedo time series with two linear functions using a variable breaking point.
The GLASS albedo shows significant darkening (p < 0.01) of the surface of the GrIS for the 1996-2012 period, with the summer (JJA) albedo declining at a rate of 0.02 ± 0.004 decade −1 (Fig. 1a). About 25 % of this decline might be attributed to sensor degradation, per the analysis of Stroeve et al. (2013). However, the Terra sensor degradation is spectrally dependant and temporally non-linear (Wang et al., 2012). This, together with the fact that the GLASS product uses a combination of both Terra and Aqua data (which reduces the impact of the Terra sensor degradation) indicates that the impact of the sensor degradation on the observed decline is much smaller than 25 %. Over the same period, MAR-simulated summer near-surface temperature increased at a rate of 0.74 ± 0.5 • C decade −1 (Fig. 1b, p < 0.05), consistent with observed enhanced surface melting (e.g. Fettweis et al., 2013). MAR simulations also point to positive trends between 1996 and 2012 in summer surface grain radius (0.12 ± 0.03 mm decade −1 , p < 0.01, Fig. 1c) and the extent of those regions where bare ice is exposed during summer (380 ± 190 km 2 decade −1 , p < 0.01, Fig. 1d). There is no statistically significant trend in GLASS summer albedo or MAR-simulated surface grain size and bare ice extent for the 1981-1996 period. Simulated  the 1981-2012 mean. We suggest that for 2010-2012, in addition to surface melting, reduced summer snowfall might have played a key role in the accelerated decline in summer albedo.

Drivers: surface grain size and bare ice
Interannual variability in the mean summer GLASS albedo is captured by the MAR albedo simulations (Fig. 1a). For the period when the darkening has been identified, MAR albedo values explain ∼ 90 % (de-trended) of the spaceborne-derived summer albedo interannual variability. A multilinear regression analysis indicates that, over the same period, the interannual variability of summer values of surface grain size and bare ice extent simulated by MAR explain, respectively, 54 % (grain size) and 65 % (bare ice) of the in-terannual variability of GLASS albedo when considered separately. When linearly combined, grain size, bare ice extent, and snowfall explain ∼ 85 % of the GLASS interannual variability, with the influence of summer new snowfall alone explaining only 44 % of the GLASS summer albedo variability. The spatial distribution of observed summer albedo trends from space shows that the largest trends (in magnitude) occur over those regions where surface temperature, grain size, and bare ice exposure have also changed the most (Fig. 2). In particular, darkening observed from space is most pronounced at lower elevations in southwest Greenland, with trends as large as −0.20 ± 0.07 decade −1 (Fig. 2a; note that the colour bar only goes down to −0.06 decade −1 for graphical purposes), where trends in the number of days when simulated surface temperature exceeds 0 • C (Fig. 2b)  ( Fig. 2c), and the number of summer days when bare ice is exposed (Fig. 2d) are the largest. While MAR is able to capture a large component of the observed variability in albedo retrieved by GLASS, the simulated albedo trend is smaller in magnitude than that estimated using the GLASS product. The largest differences occur along the southwest margin of the ice sheet (Fig. 3), where a "dark band" of outcropping layers of ice containing large concentrations of LAI is known to be present on the surface (Wientjes et al., 2011). In this region the number of days when surface temperature exceeds 0 • C has increased, with trends of up to more than 20 days decade −1 along the margins of the GrIS (Fig. 1b). During this time period, GLASS albedo values are as low as 0.30, lower than that of bare ice (i.e. 0.45), consistent with in situ measured values of dirty ice (Wientjes and Oerlemans, 2010;Bøggild et al., 2010). Figure 4 shows the spatial distribution of MAR and GLASS mean JJA albedo for the year 2010 over an area centred on the dark band in southwest Greenland, as well as the time series of GLASS albedo averaged over the same ice-covered area contained within the region identified by the black rectangle in Fig. 4a. The black line in Fig. 4c shows the GLASS spatially averaged albedo within this region, with the top and the bottom margins of the grey area indicating, respectively, the maximum and minimum albedo values within that area. Note that we included only pixels that contained 100 % ice in all years (i.e. coloured areas in Fig. 4a and b) in the calculation shown in Fig. 4c, so trends are not driven by exposure of underlying land surface. Mean summer albedo from GLASS decreased over this area between 2005 and 2012 from ∼ 0.6 to ∼ 0.45 (vs. a decrease simulated by MAR of 0.075). Minimum summer albedo across all years averaged over the region is ∼ 0.4, but dips close to ∼ 0.3 in 2010, a value consistent with dirty bare ice, as shown in previous studies (Wientjes and Oerlemans, 2010;Wientjes et al., 2011;Bøggild et al., 2010). We hypothesize that the discrepancy along this dark band between MAR and GLASS albedo values is likely due to trends in the concentrations of LAI in the snow and ice in this region, which are not currently captured by the model.

Drivers: light-absorbing impurities on the surface of the GrIS
MAR simulations of albedo in different spectral bands (see Eqs. 1-4) point to comparable trends in the visible (0.3-0.8 µm; −0.009 ± 0.005 decade −1 , p < 0.05) and nearinfrared (0.8-1.5 µm; −0.010 ± 0.004 decade −1 , p < 0.05) bands (Fig. 5a) and to a much smaller trend in the shortwaveinfrared band that is not statistically significant (1.   the 16-day MODIS MCD43A3 product (Stroeve et al., 2013), which also has a visible albedo product (Fig. 5b). The MODIS albedo product we used is distributed by Boston University (https://lpdaac.usgs.gov/) and makes use of all atmospherically corrected MODIS reflectance measurements over 16-day periods, to provide an averaged albedo every 8 days. A semi-empirical bidirectional reflectance distribution function (BRDF) model is used to compute bihemispherical reflectance from these reflectance measurements (Schaaf et al., 2002). The comparison between the GLASS-and MODIS-retrieved visible albedo anomalies is shown in Fig. 5b, indicating that the two visible albedo anomalies are highly consistent, with a mean absolute error of 0.01 and a standard deviation of 0.005. There are differences in the estimated summer albedo trends from  MAR model. The underestimated darkening by MAR relative to GLASS can be attributed to several factors, including the modelled spatial and temporal variability of the exposed bare ice area and the concentration of surface LAI on the ice surface, which is currently not included in the MAR albedo scheme. A lack of impurities in the MAR albedo scheme can affect simulated albedo trends in at least two ways: first, the concentration of impurities over bare ice areas could be increasing, which would not be captured by MAR; second, the lack of impurities in the MAR albedo scheme causes bare ice areas to have an overestimated albedo. More frequent exposure of bare ice would lead to a decline in annual average albedo over time, but if the underlying bare ice is darker, such a trend would be larger. Thus, the difference in trends could result solely from an overestimation of the bare ice albedo by MAR. We are not able to discern the degree to which the difference is due to (a) errors in the area and frequency of bare ice exposure from MAR, (b) increasing concentration of impurities not captured by MAR, or (c) overestimation of albedo of an unchanging impurity-covered bare ice surface. The study by Alexander et al. (2015) suggests that bare ice albedo is, indeed, overestimated in MAR. To test the impact of a fixed bare ice albedo on the simulated albedo trend, we performed a sensitivity experiment in which daily albedo for those pixels showing bare ice exposure is reduced by a fixed value of 0.1. The magnitude of the difference in trends between the original MAR simulation (with no change on the bare ice albedo) and the one with a modified albedo (Fig. 6) being comparable to the difference between the MAR and GLASS trends (Fig. 3a) suggests that this factor alone could explain the difference. To further investigate this aspect, we test the hypothesis of increased concentration of LAI on the snow and ice surface. The concentrations of LAI in surface snow and ice can increase either because of increased atmospheric deposition or because of post-depositional processes, including (a) loss of snow water to sublimation and melt, resulting in impurities accumulating at the surface as a lag deposit (e.g. Doherty et al., 2013), and (b) the outcropping of "dirty" underlying ice associated with snow/firn removal due to ablation. These processes are themselves driven by warming, and therefore constitute positive feedbacks. Quantifying the contribution of surface LAI to GLASS summer albedo trends is a challenging task because of the relatively low impurity concentrations over most of the GrIS Bond et al., 2013), and because of known limitations related to remote sensing estimates of LAI from space (Warren, 2013). Moreover, quantifying the causes of potential increased impurity concentrations on the surface (atmospheric deposition vs. other factors) is also challenging, if not prohibitive, given the current state-of-the-art spaceborne measurements (e.g. accuracy of the satellite products) and the scarcity of in situ data. Therefore, in the next section, we look for trends in forest fires and the emissions of BC from forest fires in the main source regions for aerosols over the GrIS and assess whether atmospheric aerosol concentrations over the GrIS have increased (as a proxy for whether the deposition of aerosol has increased). . These records show that snow at these locations was significantly more polluted in the first half of the twentieth century than presently. Both these records and in situ measurements at Summit (Cachier and Pertuisot, 1994;Chýlek et al., 1995;Hagler et al., 2007;) also indicate that in recent decades, the snow in central Greenland has been relatively clean, with concentrations smaller than 4 ng g −1 for BC. This amount of BC could lower snow albedo by only 0.002 for r = 100 µm, or 0.005 for r = 500 µm (Fig. 5a of Dang et al., 2015). More recently, Polashenski et al. (2015) analysed BC and dust concentrations in 2012-2014 snowfall along a transect in northwest Greenland. They found similarly low concentrations of BC and concluded that albedo decreases in their study region are unlikely to be attributable to increases in BC or dust. Black carbon measurements from a high snowfall region of west central Greenland made on an ice core collected in 2003 show that black carbon concentrations varied significantly during the previous 215 years, with an average annual concentration of 2.3 ng g −1 during the period 1952-2002, characterized by high year-to-year variability in summer and a gradual decline in winter BC concentrations through the end of the century (McConnell et al., 2007). Snow sampled in 1983 at Dye-3 had a median of 2 ng g −1 (Clarke and Noone, 1985). In 2008 and 2010, measurements 160 km away at Dye-2, using the same method, had medians of 4 ng g −1 in spring and 1 ng g −1 in summer (Table 9 of Doherty et al., 2010).

Ice
In the absence of in situ measurements of impurity concentration trends over Greenland more broadly, or of trends in aerosol deposition rates (which are absent entirely), we investigate trends in emissions from key sources of aerosols deposited to the GrIS and trends in atmospheric aerosol optical depth (AOD) over GrIS.

Trends in fire count and BC emissions
Biomass burning in North America and Siberia is a significant source of combustion aerosol (BC and associated organics) to the GrIS (Hegg et al., 2009(Hegg et al., , 2010. Therefore, we investigated trends in the number of active fires in these two source regions, as well as BC emissions from fires in subregions within the Northern Hemisphere. For fire counts we used the MODIS monthly active fire products produced by the Terra (MOD14CMH) and Aqua sensors (MYD14CMH) generated at 0.5 • spatial resolution and distributed by the University of Maryland via anonymous ftp (http://www.fao.org/fileadmin/ templates/gfims/docs/MODIS_Fire_Users_Guide_2.4.pdf, http://modis.gsfc.nasa.gov/data/dataprod/dataproducts.php? MOD_NUMBER=_14). The results of our analysis are summarized in Fig. 7  In addition to number of fires we looked for trends specifically in BC emissions from fires in potential source regions for GrIS, using estimates from the Global Fire Emissions Database (GFED version 4.1, http://www.globalfiredata. org/). There is a great deal of interannual variability in annual BC emissions from fires in all regions (Fig. 8), with no statistically significant increase during the 1997-2012 or 1997-2014 periods from either of the Boreal source regions or from central Asia or Europe. BC emissions from fires in temperate North America increased by, on average, 0.35 × 10 9 g yr −1 during 1997-2014 and by 0.52 × 10 9 g yr −1 during 1997-2012 (p < 0.1 in both cases). However, the total BC emissions from fires in this region constitute a small fraction of that from the Boreal regions. In addition, the only statistically significant trend in regional BC emissions is a decrease in central Asia (112.6 × 10 9 g yr −1 ; p − 0.02), when GrIS albedo has declined most precipitously. Xing et al. (2013Xing et al. ( , 2015 point out that direct anthropogenic emissions have also been decreasing across almost all of the mid-to high-latitude Northern Hemisphere.

Trends in AOD over Greenland
To investigate trends in AOD over GrIS we look at AOD as simulated by models and as measured at ground-based stations at several locations around the GrIS. AOD is a measure of the total extinction (omnidirectional scattering plus absorption) of sunlight as it passes through the atmosphere,   and is related to atmospheric aerosol abundance. Thus, it is a metric for the mass of aerosol available to be potentially deposited onto the GrIS surface. In the aerosol models, we are able to examine trends in total AOD as well as in aerosol components: BC, dust, and organic matter. In addition, we examined trends in modelled deposition fluxes of these species to the GrIS.
For our analysis, we used model results from the Aerosol Comparisons between Observations and Models (AeroCom) project, an open international initiative aimed at understanding the global aerosol and its impact on climate (Samset et al., 2014;Myhre et al., 2013;Jiao et al., 2014;Tsigaridis et al., 2014). The project combines a large number of observations and outputs from fourteen global models to test, document, and compare state-of-the-art modelling of the global aerosol. We specifically show standardized (i.e. subtracting the mean and then dividing by the standard deviation) deposition fluxes of BC, dust, and organic aerosol (OA) from the GISS modelE contribution to the AeroCom phase II series of model runs (http://aerocom. met.no/aerocomhome.html). The runs used here took as input the decadal emission data from the Coupled Model Intercomparison Project Phase 5 (CMIP5). In this case, we report the outputs of the NASA GISS ModelE obtained from the AeroCom. In particular, Figs. 9 and 10 show modelled deposition fluxes at the two locations of Kangerlussuaq (Fig. 9, 67 • 00 31 N, 50 • 41 21 W) and Summit (Fig. 10, 72 • 34 47 N, 38 • 27 33 W) for the months of June, July, and August and aerosol components (BC, dust, and organic matter). These locations were selected as representative of the ablation zone (Kangerlussuaq) and the dry-snow zone (Summit). The analysis of the NASA GISS ModelE AeroCom outputs shows no statistically significant trend in the modelled fluxes for either location, consistent with the results recently reported by Polashenski et al. (2015) for the dry-snow zone. Results of the analysis of fluxes over different areas point to similar conclusions. Similar results are obtained when considering the months of January, February, and March, when aerosol concentration is expected to be higher. The results here presented complement other studies (e.g. Stone et al., 2014), indicating that, since the 1980s, atmospheric concentrations of BC measured at surface stations in the Arctic have decreased, with variations attributed to changes in both anthropogenic and natural aerosol and aerosol precursor emissions.  show statistically significant trends in AOD, consistent with the results of the analysis of the modelled deposition fluxes.
A recent study (Dumont et al., 2014) concluded that dust deposition has been increasing over much of the GrIS and that this is driving lowered albedo across the ice sheet. That conclusion was based on trends of an "impurity index", which is the ratio of the logarithm of albedo in the 545-565 nm MODIS band (where LAI affect albedo) to the logarithm of albedo in the 841-876 nm band (where they do not). In the MODIS product used in the study by Dumont The Cryosphere, 10, 477-496, 2016 www.the-cryosphere.net/10/477/2016/ et al. (2014), albedo values rely on removal of the effects of aerosols in the atmosphere. In the Dumont et al. (2014) study, this correction was made using simulations of atmospheric aerosols by the Monitoring Atmospheric Composition and Climate (MACC) model. Their resulting impurity index shows positive trends, and these are attributed in part (up to 30 %) to increases in atmospheric aerosol not accounted for by the model, and the remainder to increases in snow LAI. The latter is consistent with our findings herein: that GrIS darkening is in part attributable to an increase in the impurity content of surface snow. However, Dumont et al. (2014) assume that this increase in surface snow LAI is a result of enhanced deposition from the atmosphere. They do not account for the possibility that positive trends in impurity content may instead be a result of warming-driven insnow processes. Indeed, their own table shows variable AOD at AERONET stations in Greenland, but no trend over the period studied (2007)(2008)(2009)(2010)(2011)(2012).
The results of the analysis discussed above reinforce our argument that the decline in the visible albedo over Greenland is probably not due to an increase in the rate of deposition of LAI from the atmosphere, but instead is due to the consolidation of LAI at the snow surface with warmingdriven increases in melt and/or sublimation and with the increased exposure of underlying dirty ice.

Albedo projections through to 2100
We estimated future projections of summer albedo over the GrIS using MAR forced with the outputs of three different Earth system models (ESMs) from CMIP5, driven by two radiative forcing scenarios (Meinshausen et al., 2011) over the 120-year period 1980-2100. The first scenario corresponds to an increase in the atmospheric greenhouse gas concentration to a level of 850 ppm CO 2 equivalent (RCP45); the second scenario increases CO 2 equivalent to > 1370 ppm in 2100 (RCP85) (Moss et al., 2010;Meinshausen et al., 2011). The three ESMs used are the second generation of the Canadian Earth System Model (CanESM2), the Norwegian Community Earth System Model (NorESM1), and the Model for Interdisciplinary Research on Climate (MIROC5) of the University of Tokyo, Japan. More information is available in Tedesco and Fettweis (2012). The ESMs are used to generate MAR outputs for the historical period    Figure 11. Projections of broadband albedo anomaly (with respect to the year 2000) averaged over the whole GrIS for 1990-2012 from MAR simulations and GLASS retrievals (black and red lines, respectively), and as projected by 2100. Future projections are simulated with MAR forced at its boundaries with the outputs of three ESMs under two warming scenarios, with the first scenario (RCP45) corresponding to an increase in the atmospheric greenhouse gas concentration to a level of 850 ppm CO 2 equivalent by 2100 and the second (RCP85) to > 1370 pm CO 2 equivalent. The top and the bottom of the coloured area plots represent the results concerning the RCP45 (top panels) and RCP85 (bottom panels) scenarios. Semi-transparent colours are used to allow overlapping data to be viewed. Dark green corresponds to the case where MIROC5 and CANESM2 results overlap and brown to the case when the results from the three ESMs overlap. jor difference from the standard CESM configuration concerns a modification to the treatment of atmospheric chemistry, aerosols, and clouds (Seland et al., 2008) and the ocean component. Lastly, MIROC5 is a coupled general circulation model developed at the Center for Climate System Research (CCSR) of the University of Tokyo, composed of the CCSR/NIES (National Institute of Environmental Studies) atmospheric general circulation model (AGCM 5.5) and the CCSR Ocean Component Model, including a dynamicthermodynamic sea-ice model (e.g. Watanabe et al., 2010Watanabe et al., , 2011. We refer to Tedesco and Fettweis (2012) for the evaluation of the outputs of MAR when forced with the outputs of the ESMs during the historical period . All simulations consistently point to darkening accelerating through the end of the century (Fig. 11), with summer albedo anomalies (relative to the year 2000) as large as −0.08 by the end of the century over the whole ice sheet, and even greater (−0.1) over the western portion of the ice sheet (Fig. 12). The magnitude of the projected albedo anomalies by 2100, however, is probably underestimated by our simulations, because (a) the model tends to underestimate melting when forced with the ESMs , and therefore underestimates grain size growth, and (b) the model currently does not account for the presence of LAI in the snow or on the ice surface, nor for the positive feedback between LAI and snow/ice melt.

Discussion
Our results show a darkening of the GrIS 1996-2012, and indicate that this darkening is associated with increased surface snow grain size, an expansion in the area and persistence of bare ice, and by an increase in surface snow light-absorbing impurity (LAI) concentrations. We find no evidence for general increases in the deposition of LAI across the GrIS, so we associate the higher surface snow impurity concentrations predominantly with the appearance of underlying dirty ice and the consolidation of LAI in surface snow resulting from snow melt. Interannual variability in the JJA GLASS albedo is captured by the MAR albedo simulations, with the latter explaining ∼ 90 % of the space-borne-derived albedo interannual variations for the period 1996-2012. The strong correlation between MAR and GLASS albedo time series for this period suggests that MAR is capturing the processes driving most of the albedo interannual variability (grain size metamorphism and bare ice exposure) and that these processes have more influence than those associated with the spatial and temporal variability of surface impurity concentrations at seasonal timescales (currently not included in the MAR albedo scheme). This is reinforced by the fact that the range of snow grain size found across the GrIS produces larger changes in albedo than does the range of LAI concentrations measured over the GrIS, at least in the cold-snow and percolation zones of the ice sheet. As pointed out by Tedesco et al. (2015), for pure snow, grain growth from new snow (with r = 100 µm) to old melting snow (r = 1000 µm) can reduce broadband albedo by ∼ 10 %. By comparison, adding 20 ng g −1 of BC, which has been found in the top layer of melting GrIS snow, reduces albedo by only 1-2 %, consistent with the results reported by Polashenski et al. (2015).
Modelled (MAR) and retrieved (GLASS) albedo are compared, with the latter showing stronger declines in GrIS albedo, particularly over the ablation zone. Based on our analysis, we suggest that the difference between MAR and GLASS trends cannot be driven solely by the MODIS sensor degradation on the Terra satellite (also used in the GLASS product), because the estimated impact of sensor degradation on the albedo trend is much smaller than the difference between the MAR and GLASS trends, and because the GLASS product is obtained by combining data from both Terra and Aqua satellites, hence likely reducing the impact of the Terra sensor degradation on the trends. This is especially true over the dark zone, where substantial melting occurs and where the albedo decline is pronounced. As mentioned, a lack of impurities in the MAR albedo scheme can affect simulated albedo trends in at least two ways: first, the concentration of impurities over bare ice areas could be increasing or/and the lack of impurities in the MAR albedo scheme causes bare ice areas to have an overestimated albedo. Moreover, more frequent exposure of bare ice would lead to a decline in annual average albedo over time, with such a trend being larger in the case of the presence of impurity concentrations on the ice surface. Our sensitivity analysis of the simulated trends on the bare ice albedo value indicates that the difference between MAR and GLASS estimated trends is consistent with a relatively darker (e.g. containing LAI) bare ice. Since MAR does not account for the presence of surface LAI, and because the impact of LAI is mostly in the UV and visible portion of the spectrum, we suggest that another mechanism explaining the difference of −0.017 decade −1 between the MAR and GLASS visible albedo trends is associated with increasing mixing ratios of LAI in surface snow and ice on some parts of the GrIS. As we pointed out, this could be due to a combination of increased exposure of dirty ice with ablation (Wientjes and Oerlemans, 2010;Bøggild et al., 2010), to enhanced melt consolidation with warming (e.g. Doherty et al., 2013), or to increased deposition of LAI from the atmosphere. The absence of in situ spatially distributed measurements to separate these processes means that we cannot quantify their relative contributions to the darkening in the visible region. Based on our analysis of trends in AOD over Greenland and the lack of a trend in forest-fire counts and BC in North America and Eurasia, we argue that increased deposition of LAI is not a large driver for the observed negative trends in Greenland surface albedo. An exception could be an increase in the deposition of locally transported dust near the glacial margins, which would primarily affect the ablation zone. In particular, locally lofted dust may be playing a substantial role in the southwest GrIS ablation zone. However, we note that increased deposition is not needed in order to have an increase in the concentration of LAI at the GrIS surface. As noted above, indeed, temperatures and melt rates have been accelerating over the GrIS during the past decades (e.g. Tedesco et al., 2014). When snow melts, snow water is removed from the surface more efficiently than particulate impurities; the result is an increase in impurity concentrations in surface snow (e.g. Flanner et al., 2007;Doherty et al., 2013). Large particles, such as dust, in particular, will have poor mobility through the snowpack (Conway et al., 1996) so their concentration at the surface is expected to increase with snowmelt. This effect may be especially amplifying snow impurity content in the low-altitude ablation zone of the GrIS, where enhanced melting has been occurring (e.g. Tedesco et al., 2014). Further, the albedo reduction for a given concentration of an absorbing impurity in snow is greater in large-grained snow than in small-grained snow (Fig. 7 of Warren and Wiscombe, 1980;, so climate warming itself will amplify the effect of LAI on surface albedo. Warming may also lead to increased sublimation, removing snow water but not particles from the snow surface, again increasing concentrations of LAI in surface snow. Snow and ice warmed by increased temperatures and higher LAI concentrations also promotes darkening via socalled "bio albedo", with biological growth on the surface depressing the albedo. Green, pink, purple, brown, and black pigmented algae, indeed, occur in melting snow and ice. Microbes can bind to particulates, including BC, retaining them at the surface in higher concentrations than in the parent snow and ice. The magnitude of this source of darkening is currently unquantified, but as the climate warms and melt seasons lengthen, biological habitats are expected to expand, with their contribution to darkening likely increasing (Benning et al., 2014).
Quantifying the impact of aerosols on Greenland darkening is also made difficult by the large disagreements among models in their predicted aerosol deposition rates over the GrIS. We examine the contrast between AOD trends from the MACC model used by Dumont et al. (2014) and the Goddard Chemistry Aerosol Radiation and Transport model (GO-CART). The GOCART model simulates major tropospheric aerosol components, including sulfate, dust, BC, organic carbon, and sea-salt aerosols using assimilated meteorological fields of the Goddard Earth Observing System Data Assimilation System (GEOS DAS), generated by the Goddard Global Modeling and Assimilation Office. Figure 13 compares results for AOD at 550 nm from MACC and GOCART for dust, organic matter, and black carbon for the domain bounded by 75 to 80 • N and 30 to 50 • W (the same area considered by Dumont et al., 2014). The MACC model shows statistically significant trends for dust (p < 0.01) and for total aerosols (p < 0.05). All remaining trends are not statistically significant for both MACC and GOCART outputs (Fig. 13).
Neither model represents the process of increased exposed silt/dust as Greenland glaciers recede; therefore, we would not expect them to capture trends in dust from this source. The inconsistency between the MACC and GOCART values and trends is puzzling, and indicates that the simulation of aerosol deposition rates over Greenland needs improvement.

Summary and conclusions
We studied the mean summer broadband albedo over the Greenland ice sheet between 1981 and 2012 as estimated from space-borne measurements and found that summer albedo decreased at a rate of 0.02 decade −1 between 1996 and 2012. The analysis of the outputs of the MAR regional climate model indicates that the observed darkening is associated with increasing temperatures and enhanced melting occurring during the same period, which in turn promote increased surface snow grain size as well as the expansion and persistency of areas with exposed bare ice. The MAR model simulates the interannual variability in the retrieved GLASS albedo well, but the albedo trend is larger in the GLASS albedo product than in MAR, indicating that processes not represented in the MAR physics account for some of the declining albedo. Specifically, we suggest that the absence of the effects of light-absorbing impurities in MAR could account for the difference. We also suggest that this hypothesis is supported by the trends observed along the ablation zone, where the differences between observed and modelled trends are more pronounced and the effect of the Terra sensor degradation plays a relatively small role. On the other hand, over the dry-snow zone, our hypothesis requires further testing, in view of the potentially higher impact of the sensor degradation on the observed albedo trend. The analysis of modelled fields and in situ data indicated an absence of trends in aerosol optical depth over Greenland, as well as no significant trend in particulate light-absorbing emissions (e.g. BC) from fires in likely source regions. This is consistent with the absence of trends in surface aerosol concentrations measured around the Arctic. Consequently, we suggest that the increased surface concentrations of LAI associated with the darkening are not related to increased deposition of LAI, but rather to post-depositional processes, including increased loss of snow water to sublimation and melt and the outcropping of "dirty" underlying ice associated with snow/firn removal due to ablation. Future projections of GrIS albedo obtained from MAR forced under different warming scenarios point to continued darkening through the end of the century, with regions along the edges of the ice sheet subject to the largest decrease, driven solely by warming-driven changes in snow grain size, exposure of bare ice, and melt pool formation. We hypothesize that projected darkening trends would be even greater in view of the underestimated projected melting (and effect on albedo) and in view of the fact that the current version of the MAR model does not account for the presence of surface LAI and the associated positive direct and indirect impact on lowered albedo.
The drivers we identified to be responsible for the observed darkening are related to endogenous processes rather than exogenous ones and are strongly driven by melting. Because melting is projected to increase over the next decades, it is crucial to assess our capability of studying, quantifying, and projecting these processes as they will inevitably impact, and be impacted by, future scenarios. Intrinsic limitations of current observational tools and techniques, the scarcity of in situ observations, and the albedo schemes currently used in existing models of surface energy balance and mass balance limit our ability to separate the contributions to darkening by the different processes, especially with regard to the cause and evolution of surface impurity concentrations. Moreover, as with all instruments, sensors undergo deterioration, and it can be difficult to separate an albedo trend from sensor drift. This is especially true in the dry-snow zone, where impurity concentrations are extremely low (only a few parts per billion (ppb) in the case of BC). In this regard, a recent study by Polashenski et al. (2015) suggests that the decline and spectral shift in dry-snow albedo over Greenland contains important contributions from uncorrected Terra sensor degradation when using the MODIS data collection C5. The new MODIS Terra version (accounting for the sensor degradation) does not appear to show any trend (C. Polashenski, personal communication, 2015), hence supporting the hypothesis of the absence of trends of LAI deposition over the dry zone.
Remote sensing and in situ observations should be complemented with models that simulate the surface energy balance to account for the evolution of the snowpack, in particular changes in surface grain size and exposure of bare ice. Simulations with regional climate models can provide such quantities, but they do not currently account for the transport and deposition of LAI to Greenland, the post-depositional evolution of impurities in the snowpack, and the synergism between surface LAI and grain growth (whereby a given impurity content causes more albedo reduction in coarsegrained snow than in fine-grained snow). In this regard, the current parameterization for snow albedo in MAR is based on that of Brun et al. (1992), as part of an avalanche-forecasting model. As a consequence of the results of this study, we began evaluating an alternative albedo scheme using a parameterization that can also account for the albedo reduction by absorptive impurities (e.g. Dang et al., 2015). Moreover, we are also considering using the firn/ice albedo parameterization of Dadic et al. (2013), based on measurements covering the range of densities from 400 to 900 kg m −3 . Surface-based measurements are needed to test satelliteretrieved albedo and to quantify the drivers behind albedo changes in different areas of Greenland. To date, most surface-based observations have been made in the dry-snow zone or the percolation zone, and they have generally focused on measuring the mixing ratios of BC (Hagler et al., 2007;McConnell et al., 2007;Polashenski et al., 2015) or of the spectral light absorption by all particulate components collectively Hegg et al., 2009Hegg et al., , 2010. The regions of Greenland that are darkening the most rapidly are within the ablation zone. Here, there is no direct evidence that the rate of atmospheric deposition of LAI has been increasing. In view of the cumulative effect of snowmelt leaving impurities at the surface, the intra-seasonal variation of deposition may not be as important as the exposure of LAI by melting. Changes in the abundances of light-absorbing algae and other organic material with warmer temperatures may also be contributing to declining albedo, particularly for the ice, but this is an essentially unstudied source of darkening. Until measurements are made that quantify and distinguish the relative roles of each of these factors in the darkening of the GrIS, it is not possible to reduce the uncertainty in their contributions to the acceleration of surface melt. In addition to the need for targeted ground observations, it is necessary for the models that simulate and project the evolution of surface conditions over Greenland to start including the contribution of surface LAI, their processes, and their impact on albedo, as well as aerosol models that account for their deposition. Concurrently, space-borne sensors or missions capable of separating the contributions from the different processes (with increased spatial, spectral, and radiometric resolution) should be planned for remote sensing to become a more valuable tool in this regard.