Assessing spatio-temporal variability and trends in modelled and measured Greenland Ice Sheet albedo ( 2000 – 2013 )

Accurate measurements and simulations of Greenland Ice Sheet (GrIS) surface albedo are essential, given the role of surface albedo in modulating the amount of absorbed solar radiation and meltwater production. In this study, we assess the spatio-temporal variability of GrIS albedo during June, July, and August (JJA) for the period 2000–2013. We use two remote sensing products derived from data collected by the Moderate Resolution Imaging Spectroradiometer (MODIS), as well as outputs from the Modèle Atmosphérique Régionale (MAR) regional climate model (RCM) and data from in situ automatic weather stations. Our results point to an overall consistency in spatio-temporal variability between remote sensing and RCM albedo, but reveal a difference in mean albedo of up to ∼ 0.08 between the two remote sensing products north of 70 N. At low elevations, albedo values simulated by the RCM are positively biased with respect to remote sensing products by up to ∼ 0.1 and exhibit low variability compared with observations. We infer that these differences are the result of a positive bias in simulated bare ice albedo. MODIS albedo, RCM outputs, and in situ observations consistently indicate a decrease in albedo of−0.03 to−0.06 per decade over the period 2003–2013 for the GrIS ablation area. Nevertheless, satellite products show a decline in JJA albedo of −0.03 to −0.04 per decade for regions within the accumulation area that is not confirmed by either the model or in situ observations. These findings appear to contradict a previous study that found an agreement between in situ and MODIS trends for individual months. The results indicate a need for further evaluation of high elevation albedo trends, a reconciliation of MODIS mean albedo at high latitudes, and the importance of accurately simulating bare ice albedo in RCMs.


Introduction
Over the past decade, the Greenland Ice Sheet (GrIS) has simultaneously experienced accelerating mass loss (van den Broeke et al., 2009;Rignot et al., 2011) and records for the extent and duration of melting (Tedesco et al., 2008(Tedesco et al., , 2011(Tedesco et al., , 2013;;Nghiem et al., 2012).Increased melt over Greenland has been associated with both changes in temperature and an amplifying ice-albedo feedback: increased melting and bare ice exposure reduce surface albedo, thereby increasing the amount of absorbed solar radiation and, in turn, further amplifying melting (Box et al., 2012;Tedesco et al., 2011).Recent studies (van den Broeke et al., 2011;Vernon et al., 2013) also indicate that albedo plays an essential role in the GrIS surface energy balance, and consequently, the surface mass balance (SMB) of regions where considerable melting occurs.Due to the impact of albedo on the surface energy balance, it is crucial to assess the performance of models that simulate albedo over the GrIS and the quality of albedo estimates from remote sensing or in situ observations.These assessments are pivotal for improving our understanding of the physical processes leading to accelerating mass loss, and for improving future projections.

P. M. Alexander et al.: Variability and trends in Greenland Ice Sheet albedo
Several studies investigating GrIS albedo trends and variability have primarily relied on satellite measurements, particularly those collected by the Moderate Resolution Imaging Spectroradiometer (MODIS; Stroeve et al., 2005Stroeve et al., , 2006Stroeve et al., , 2013;;Box et al., 2012).Remote sensing measurements can continuously capture changes at large spatial scales and for long periods, with the exception of cases when the surface is obscured by clouds.Previous studies have found MODIS albedo products to agree reasonably well with in situ data, especially with regards to capturing the seasonal albedo cycle and mean seasonal values in regions where variability is small (Stroeve et al., 2005(Stroeve et al., , 2006(Stroeve et al., , 2013)), but lower accuracy at high solar zenith angles has been identified (Stroeve et al., 2005(Stroeve et al., , 2006)), limiting the periods and locations for which these data can be used.Nevertheless, given their relatively high temporal and spatial resolution, these products are useful for evaluating albedo derived from regional climate models (RCMs).RCMs are an important tool for estimating both current and future changes in the GrIS SMB (Box and Rinke, 2002;Box et al., 2006;Ettema et al., 2009;Fettweis et al., 2007Fettweis et al., , 2011;;Rae et al., 2012;Tedesco and Fettweis, 2012), and the surface albedo schemes employed by these models have a substantial impact on their simulation of the SMB (Rae et al., 2012;van Angelen et al., 2012;Lefebre et al., 2005;Franco et al., 2012).
In this paper, we report the results of an assessment of GrIS albedo spatio-temporal variability and trends for the period 2000-2013.To our knowledge, this is the first time a multitool integrated assessment of albedo over Greenland is presented.We use (1) data from two remote sensing products from MODIS, the MOD10A1 daily albedo product (Hall et al., 2012) and MCD43A3 16-day albedo product (Schaaf et al., 2002), (2) in situ albedo data from the Greenland Climate Network (GC-Net, Steffen et al., 1996) and Kangerlussuaq-Transect (K-Transect;van de Wal et al., 2005), and (3) outputs from two versions (v2.0 and v3.2) of the Modèle Atmosphérique Régionale (MAR; Fettweis et al., 2013a, b).In order to carry out comparisons between products, MODIS data have been re-gridded to the MAR model grid and, where necessary, daily data have been averaged over 16-day periods.The role of potential errors associated with differences in spectral range between satellite and in situ data and cloud cover have been considered and corrected for where possible.

The MAR model
The MAR model (Gallée and Schayes, 1994;Gallée, 1997;Lefebre et al., 2003) is a coupled land-atmosphere RCM featuring the atmospheric model described by Gallée and Shayes (1994) and the Soil Ice Snow Vegetation Atmosphere Transfer scheme (SISVAT) surface model.SISVAT incorporates the multilayer snow model Crocus (Brun et al., 1992), ) and locations of all GC-Net and K-Transect weather stations.Pixels not defined as 100 % ice covered in MAR v3.2 are masked out.The dotted black line is the mean equilibrium line (where mean SMB is 0).The K-transect stations (S5, S6, S9, S10) are red, while GC-Net stations are black.Stations in grey are GC-Net stations not used in this study.Other contour lines indicate elevation in metres above sea level.The inset shows individual stations near the west coast ablation zone.
which simulates fluxes of mass and energy between snow layers, and reproduces snow grain properties and their effect on surface albedo.The model setup used here is described in detail by Fettweis (2007).We primarily used a recent version of MAR (v3.2), which features changes to the albedo scheme relative to previous versions (v1.0 and v2.0), detailed in Sect.2.2, but also examined differences between MAR v3.2 and a previous version, MAR v2.0.MAR v3.2 (v2.0) has been run at a 25 km horizontal resolution for the period 1958-2013 (1958-2012).Both model versions were forced at the lateral boundaries and ocean surface and initialized with 6-hourly reanalysis outputs from the European Centre for Medium-Range Weather Forecasts (ECMWF), using the ERA-40 reanalysis (Uppala et al., 2005) for the period 1958-1978 and the ERA-Interim reanalysis (Dee et al., 2011) for the period 1979-present.Here we focus on the 2000-2013 period for comparison with satellite data.The MAR v3.2 ice sheet mask (which gives the fraction covered by ice for each grid box) and surface elevation were defined using the Greenland digital elevation model of Bamber et al. (2013).MAR v2.0 uses the elevation model of Bamber et al. (2001), and the land surface classification mask from Jason Box (http://sites.google.com/site/jboxgreenland/datasets).
In contrast with MAR v2.0, MAR v3.2 sub-grid scale parameterizations make it possible to have fractions of different land cover types within a single grid box.Quantities were computed for the sectors within each grid box and a weighted average of these quantities was used to represent the average value for a grid box.
For convenience, mean SMB for September 2000-August 2013 from MAR v3.2 is shown in Fig. 1, along with the equilibrium line dividing positive and negative SMB, together with the locations of the weather stations used in this study.In this study, areas below the mean 2000-2013 equilibrium line as defined by MAR are collectively referred to as the "ablation area", while areas above this line are referred to as the "accumulation area".

The MAR albedo scheme
The basis for the MAR albedo scheme is described in detail by Brun et al. (1992) and Lefebre et al. (2003).MAR snow albedo (α) depends on the optical diameter of snow grains (d), which is in turn a function of other snow grain properties, such as grain size, sphericity and dendricity.In the model, the sphericity, dendricity, and size of snow grains are a function of snowpack temperature, temperature gradient, and liquid water content.Albedo is defined in MAR for three spectral intervals: Interval 1, visible (0.3 − 0.8 µm) : Interval 2, near infrared (0.8 − 1.5 µm) : Interval 3, far infrared (1.5 − 2.8 µm) : where α 1 , α 2 , and α 3 are wavelength dependent albedo values.The integrated snow albedo (α S ) for the range 0.3 to 2.8 µm is a weighted average of albedo over these intervals based on solar irradiance fractions: The minimum albedo of snow is set to 0.65.In MAR v2.0, bare ice albedo is simply assigned a fixed value.In MAR v3.2 (the version primarily used here), bare ice albedo is a function of accumulated surface water following the parameterizations of Lefebre et al. (2003), described below.In the case of bare ice, which occurs in MAR when the surface snow density is greater than 920 kg m −3 , ice albedo (a I ) is given by where α I,min and α I,max are the minimum and maximum bare ice albedo, K is a scale factor (set to 200 kg m −2 ), and M SW (t) is the time-dependent accumulated amount of excessive surface meltwater before run-off (in kg m −2 ).According to the parameterization of Zuo and Oerlemans (1996), there is delay in MAR v3.2 between the production of meltwater and evacuation towards the oceans (Lefebre et al., 2003), in order to account for the reduction of bare ice albedo due to the presence of surface water.The ice surface albedo (α I ) will therefore be lower if the melt rate is higher, asymptotically approaching the minimum bare ice albedo.Additionally, to ensure temporal continuity in simulated albedo, values of albedo between the maximum bare ice and minimum snow albedo are possible when the surface (or near-surface) snow density lies between 830 and 920 kg m −3 .In this case (which corresponds to the presence of firn), α I is a function of density as follows (Lefebre et al., 2003): where α S,min is the minimum albedo of snow, ρ I is the density of the upper firn layer, and ρ C is the density at which pores within firn close off (830 kg m −3 ).Table 1 provides the range of possible albedo values for ice, firn, and snow in MAR v2.0 and v3.2.
In cases where there is a snowpack with a thickness of < 10 cm overlaying ice or firn (with a density greater than 830 kg m −3 ), excluding the case of ice lenses, albedo is interpolated between the ice albedo and the surface snow albedo as a linear function of snowpack thickness, to produce an "integrated" surface albedo of snow, ice and water (α SI ): where H s is the snowpack height in metres.In cases where a snowpack thicker than 0.1 m lies above ice, or there is bare ice at the surface, α SI is simply set equal to α s or α I , respectively.

P. M. Alexander et al.: Variability and trends in Greenland Ice Sheet albedo
University (Schaaf et al., 2002; available at https://lpdaac.usgs.gov/),available for the same period.The MOD10A1 Version-5 product contains daily albedo (0.3-3 µm) based on the "best" daily MODIS observation, defined as the observation that covers the greatest percentage of a grid cell.Corrections are also applied to account for anisotropic scattering, for the influence of the atmosphere on surface albedo, and for the limited spectral range of MODIS bands (Klein and Stroeve, 2002;Stroeve et al., 2006).Here we used MODIS data from the TERRA satellite, as MODIS data from the AQUA satellite are less reliable due to an instrument failure in the near infrared band (Stroeve et al., 2006;Box et al., 2012).
The MCD43A3 Version-5 product makes use of all atmospherically-corrected MODIS reflectance measurements over 16-day periods to provide an integrated albedo measurement every eight days.A semi-empirical bidirectional reflectance function (BRDF) model is used to compute bi-hemispherical reflectance as a function of these reflectance measurements (Schaaf et al., 2002).The MCD43A3 product contains, in addition to albedo values for each MODIS instrument band, "short-wave" (SW) albedo values calculated over a wavelength interval of 0.3-5.0µm and "visible" albedo values for the 0.3-0.7 µm interval, calculated using the BRDF parameters.Here we primarily made use of SW MCD43A3 albedo, as its wavelength interval is consistent with those of MAR and MOD10A1, but briefly considered "visible" albedo as well.The MCD43A3 product provides, over each wavelength interval, an integrated diffuse white-sky albedo (WSA) and a direct black-sky albedo (BSA) for a specific viewing geometry (from above when the local solar zenith angle is at a maximum).A linear combination of WSA and BSA can be used to compute the true blue-sky albedo.Stroeve et al. (2005) suggest that there is little difference between BSA and WSA for typical summer midday solar zenith angles over Greenland.Simulation of blue-sky albedo requires models or observations of aerosol optical depth (Stroeve et al., 2013) that are not available for this study and, therefore, the following results consider BSA only.
Both MODIS products provide quality flags indicating "good quality" vs. "other quality" data.In the case of MCD43A3, "other quality" data were produced using a backup algorithm.When few observations were available, the backup algorithm was used to scale an archetypal BRDF based on past observations (Schaaf et al., 2002).In order to understand the influence of data quality on our results, we present results for both "all quality" as well as "good quality" data.

Weather station data
We used automatic weather station (AWS) data from two sources, GC-Net (Steffen et al., 1996) and the K-Transect (van de Wal et al., 2005).The locations of the weather sta-tions are shown in Fig. 1, and a list of the weather stations used and their period of coverage is provided in Table 2.We used all available GC-Net and K-transect June-July-August (JJA) data within the period 2000-2012 for comparison with MODIS and MAR albedo.GC-Net data for the summer of 2013 were not available when data analysis for this study was conducted.We followed a procedure similar to that used by Stroeve et al. (2005) to generate an albedo time series from GC-Net and K-Transect data.Mean daily albedo was computed as the sum of daily incident SW radiation divided by the sum of daily outgoing SW radiative flux.Instances where hourly upward SW radiative flux exceeded downward SW radiative flux were excluded.Upward and downward hourly radiative fluxes were excluded when downward fluxes were smaller than 250 W m −2 to reduce the impact of relative errors on measured albedo, especially during cases of low incident radiation (we investigated the sensitivity of our results to this threshold, and did not find a considerable effect on the results).Data from several locations and time periods were excluded from this analysis.These stations are listed in Table 2.In particular, measured albedo at Swiss Camp for the year 2012 appeared to be unrealistically high relative to previous years and was excluded (mean measured JJA albedo for 2012 was 0.99, 3.5 standard deviations above the mean 2000-2011 value of 0.64).Measured albedo at Crawford Point-2 underwent a step change after 2004 (mean albedo was 0.81 ± 0.03 for 2000-2004 and 0.90 ± 0.04 for 2005-2010) and, therefore, we excluded data after 2004, as was done by Stroeve et al. (2013).At stations NASA-U and NGRIP, levelling errors produced a low bias in upward radiation for all years (Stroeve et al., 2013), resulting in measured albedo values that were unrealistically low for snow outside of the ablation area (0.30±0.01 at NASA-U and 0.33±0.01 at NGRIP), and were excluded.At the Peterman Glacier and Peterman ELA, missing MODIS data prevented the inclusion of all weather station data in this analysis.

Corrections to MAR albedo
Snow albedo is generally higher during cloudy conditions due to the masking of a portion of the incoming solar spectrum by clouds (Greuell and Konzelman, 1994).Both MAR v3.2 and MAR v2.0 account for this by applying a correction to α SI as a function of cloud fraction, following Greuell and Konzelman (1994): where n is the cloud fraction computed by MAR, and α CL is the cloud-corrected albedo.Satellite data can only provide cloud-free measurements, and we therefore re-corrected MAR surface albedo to produce estimates of cloud-free surface albedo.This particular technique was used, rather than excluding pixels from MAR, because MAR does not  [2002][2003][2004][2007][2008][2009] Accumulation area, south of 70  necessarily replicate the actual cloud fraction observed by MODIS.The correction applied here reverses the correction applied in MAR, then corrects albedo for the case where there is a cloud fraction of 0: In this case, α MAR,daily is the daily mean MAR albedo, n MAR,daily is the daily mean cloud fraction from MAR, and α MAR,clear−sky is the daily mean clear-sky albedo.All analyses with MAR results were conducted using α MAR,clear−sky .

Aggregation of MODIS data to the MAR grid
For the purpose of comparing model results and satellite data, MODIS albedo products were re-gridded to the MAR 25 km resolution grid from the original 463 m spatial resolution at which they are distributed.Re-gridded values contain the median value of all the MODIS values falling within a MAR grid box.When comparing satellite data against model re-sults, our analysis was restricted to the GrIS.For all comparisons including MAR v3.2 results, areas where the MAR sub-grid level ice cover percentage was less than 100 % were excluded.For all comparisons including MAR v2.0 results, the same mask from MAR v3.2 was used, except pixels classified as 100 % ice covered in MAR v3.2, but classified as non-ice-covered in MAR 2.0 were also excluded from the analysis.

Comparisons at in situ stations
Comparisons at in situ stations were conducted between weather station data and data or outputs from the MODIS or MAR grid box that encompassed the in situ station.In this case, we used the original (463 × 463 m) MODIS grid box containing the station rather than the MODIS data aggregated to the encompassing 25×25 km MAR grid box, to reduce potential errors associated with spatial variations of albedo.In cases where an in situ station was contained within a MAR As in the case of the original MAR albedo outputs, in situ measurements also included measurements made during cloudy conditions while MODIS albedo data did not.Given a lack of available measurements, we did not explicitly correct in situ data for the presence of clouds in this study, but only considered data where coincident satellite and in situ measurements were available.Stroeve et al. (2013) applied a correction to GC-Net data using a radiative transfer model, but found that the correction did not significantly impact results.

Spectral differences
The GC-Net LI-COR sensors are sensitive within the 0.4-1.1 µm band, and K-Transect data are collected in the 0.3-2.8µm band.The GC-Net bands are narrower than the MOD10A1 interval of 0.3-3 µm and the MCD43 SW albedo interval of 0.3-5 µm, and the interval of 0.3-2.8µm over which albedo is calculated in the MAR model.GC-Net incoming and outgoing radiation values were calibrated to represent radiation for a spectral interval of 0.28-2.8µm (Wang and Zender, 2009).However, because snow has a high spectral reflectance over the 0.3 to 1.1 µm interval, and a much lower reflectance above 1.1 µm, measured albedo over the smaller interval will be higher for snow-covered areas (Stroeve et al., 2005).Stroeve et al. (2005) compared albedo derived with GC-Net LI-COR pyranometers to measurements from pyranometers with a larger spectral range, and found that the smaller wavelength interval results in a positive albedo bias of between 0.04 and 0.09 for GC-Net data relative to MODIS albedo, depending on the location and time period (Stroeve et al., 2005).This bias does not apply to K-Transect measurements, as the spectral sensitivity is comparable to the sensitivity of MODIS sensors.Because this bias may be smaller or larger depending on multiple factors, no correction has been applied here, however, we have provided an indication of spatial variability of this bias in Sect.4.2.1 by comparing MCD43A3 visible albedo (0.3-0.7 µm) with MCD43A3 SW (0.3-5.0 µm) albedo.

Calculation of bias, correlation, and trends
In the following analysis, we focus on the JJA period because MODIS data are less reliable during other months, when solar zenith angles are high, as discussed by Box et al. (2012), and because this is the period when surface albedo has the largest impact on SMB.
In order to compare spatial variations in albedo we calculated the mean 2000-2013 JJA MOD10A1, MCD43A3 SW BSA albedo, and MAR clear-sky albedo using all available measurements or model outputs over the specified period, ex-cluding cases where greater than 25 % of data were missing for a given pixel.When differences between data sets or between satellite data and model results were calculated, we only used measurements or results overlapping in time and space, to avoid the possibility of bias introduced by missing data.The mean difference between two samples for a given grid box was deemed to be statistically significant if the p value for a two-sample Student's t test was smaller than 0.05.Unless otherwise specified, only "good quality" MODIS data have been used in comparisons.
In some cases, observational data or model results have been spatially averaged or aggregated within the ablation and accumulation areas defined using MAR v3.2 or MAR v2.0.The ablation (accumulation) area is defined as the area that experienced a net loss (gain) of mass over the 2000-2013 period, as simulated by either version of the model.
For analyses of temporal variability, we considered daily variability, for which MOD10A1 data, in situ values, and MAR outputs were available, as well as variability over 16day MCD43A3 periods.In the case of the analysis of 16-day data, MOD10A1, MAR, and in situ daily data were averaged to produce a value for each overlapping MCD43A3 16-day period.The correlation between daily satellite data and between satellite data and model results was examined using Pearson's coefficient of determination (r 2 ).
To compare the distribution of ablation area albedo for satellite data and MAR model outputs, we produced frequency histograms for ablation area albedo using a bin width of 0.0099.Parameters for the best fit of a bimodal distribution to the histograms was obtained using the maximum likelihood estimation function in MATLAB, assuming a bimodal normal distribution for the fit.Box et al. (2012) investigated changes in GrIS albedo using the MOD10A1 albedo product, and found that between 2000 and 2012, surface albedo decreased over almost the entire ice sheet.Here, we built on the analysis of Box et al. (2012) and extended our analysis to include MCD43A3, MAR v3.2, and in situ JJA data for the period 2000-2013.Trends in albedo have been obtained by performing linear regression on 16-day albedo values for satellite products, in situ data, and model outputs, excluding albedo values outside of the JJA period.We have also computed trends for annual JJA average values.A trend was determined to be statistically different from 0 if the p value for a Student's t test was smaller than 0.05.For in situ stations, only stations with a record of at least nine years of data were included in the analysis, and only trends for albedo from the encompassing MAR v3.2 (25 × 25 km) grid box and MODIS (463 × 463 m) grid boxes over the same range of years were considered.

Albedo spatial variability
MAR v3.2 and the two MODIS data sets show coherent spatial patterns of JJA mean 2000-2013 albedo (Fig. 2) that are consistent with previous studies (e.g.Box et al., 2012), with low-elevation areas in the ablation area dominated by lower albedo values (< 0.7 on average, Table 3) due to the presence of meltwater and bare ice, and high elevation areas by relatively higher albedo (> 0.74).The most obvious discrepancy between the satellite products occurs north of 70 • N, where the MOD10A1 daily product exhibits an increase in albedo with latitude, while MCD43A3 points to the opposite.The difference between the two satellite products (Fig. 3a) is statistically significant (at the 95 % confidence level) above 70 • N, reaching ∼ 0.08 (for albedo ranging between 0 and 1) at the highest latitudes.The pattern of differences between MAR v3.2 and the two satellite products (Fig. 3b and Fig. 3c) appears to vary with both elevation and latitude, while the difference be- In each case, only coincident data for each of the two data sets being compared is used.MAR grid boxes where the difference is not statistically significant at the 95 % confidence level are marked with a grey "x".
tween the two satellite products varies primarily with latitude (Fig. 3a).Because any systematic biases in the satellite products are likely to be relatively consistent across space (at least as a function of longitude), it is likely that MAR v3.2 biases contribute to some of the elevational differences seen in Fig. 3b and Fig. 3c.Within the accumulation area south of 70 • N, MAR v3.2 albedo (0.77 on average) is comparable to MODIS albedo (average of 0.78 for MOD10A1 and 0.77 for MCD43A3).At low elevation areas, especially along the west coast ablation area, MAR v3.2 overestimates albedo (up to ∼ 0.1) relative to both satellite products.The mean ablation area albedo from MOD10A1 (0.68±0.07) is identical to MAR mean ablation area albedo (Table 3), despite the large positive bias in MAR albedo within the west coast ablation area that can be seen in Fig. 3.This is likely a result of a positive bias for MOD10A1 at high latitudes, as will be discussed further below.For areas north of 70 • N, the discrepancy between satellite products makes it impossible to determine the magnitude and direction of MAR biases.MAR v3.2, MOD10A1 and MCD43A3 mean 2000-2013 JJA albedo values show a similar logarithmic dependence of albedo with elevation (Fig. 4a); below 2000 m, albedo increases relatively rapidly with elevation (both MAR and the MODIS products show a statistically significant albedo increase of ∼ 0.01 to ∼ 0.02 per 100 m increase in elevation), while above 2000 m, the change is smaller (no statistically significant increase for MAR, and an increase of ∼ 0.002 to ∼ 0.003 per 100 m for both MODIS products).The discrepancies between MODIS products north of 70 • N are evident in Fig. 4b  For in situ stations in the ablation area (Table 4), in situ mean albedo (0.56 ± 0.08) is higher than coincident average MOD10A1 (0.51±0.09) and MCD43A3 (0.50±0.07) albedo values, and is comparable with MAR v3.2 clear-sky mean albedo for sectors classified as ice-covered (0.57 ± 0.07).Within the accumulation area, in situ albedo is larger by 0.01 to 0.06 relative to MAR and the MODIS products (Table 4).These results appear to be consistent with a positive bias in GC-Net measurements identified by Stroeve et al. (2005).However, we also find that the difference between in situ and satellite albedo is larger at K-Transect stations (+0.08) than at GC-Net sites (+0.04), and K-Transect data are not expected to exhibit the positive bias.It is likely that the high spatial variability of ablation area albedo contributes to the differences.Data from in situ stations may be positively biased relative to satellite data because of a bias introduced by station locations: locations are not chosen to be within streams, lakes, or crevasses, which have a lower albedo.In the accumulation zone, a lack of variation in surface features likely leads to smaller spatial variations in albedo.
Mean 2000-2012 JJA albedo values for ablation area GC-Net stations with a record of at least seven years do not appear to exhibit a clear variation with latitude when compared with satellite data and model results (Fig. 5).GC-Net albedo at stations north of 70 • N is on average larger by 0.02 relative to stations south of 70 • N (Table 4), suggesting that GC-Net albedo does not confirm the decrease in albedo with latitude indicated by MCD43A3.MOD10A1 accumulation area measurements are comparable (within 0.01 for aggregated station data) to uncorrected GC-Net data north of 70 • N (Fig. 5, Table 4).This suggests that the MOD10A1 may also be positively biased north of 70 • N.
It appears possible from Fig. 5 that the bias at GC-Net sites (between 0.04 and 0.09 according to Stroeve et al., 2005) could increase with latitude, rendering corrected GC-Net mean 2000-2013 albedo comparable to MCD43A3 albedo.In order to indicate how the GC-Net albedo bias is likely to vary spatially, the mean difference between MCD43A3 visible BSA (for the interval 0.3-0.7 µm) and MCD43A3 SW BSA (for the interval 0.3-5.0µm) was computed (Fig. 6).The difference is larger than the biases observed by Stroeve et al. (2005) at GC-Net stations, likely because the MCD43A3 visible wavelength interval is smaller than that for GC-Net stations.The difference does not vary with latitude.Rather, it is lowest in the ablation area where bare ice is exposed during summer months, and largest in regions where melting occurs but bare ice exposure is infrequent.The difference is relatively small at the highest elevations.
The spatial variability of the difference appears to be associated with the differences in spectral albedo between different materials.Because ice does not exhibit the spectral dependence of albedo that snow does (Hall and Martinec, 1985), the difference between MCD43A3 visible and SW albedo is lower in the ablation area where bare ice is exposed during summer.In locations where melting occurs, snow grains tend to be larger because of constructive metamorphism, reducing reflectance mostly in the near infrared band (Wiscombe and Warren, 1980), resulting in a larger difference between visible and near infrared reflectance.This suggests that in situ albedo values do not exhibit the decrease of albedo with latitude indicated by MCD43A3.

Albedo temporal variability
The standard deviation of an albedo time series provides information on the magnitude of its temporal variability.Within the low elevation ablation area of the ice sheet, both MAR and the MODIS products exhibit a relatively high standard deviation for the 2000-2013 period (0.07 on average for 16-day periods; Fig. 7, Table 3).At high elevations, variability is smaller (0.02 to 0.03 on average for 16-day periods).The MCD43A3 and MOD10A1 products show similar spatial patterns of standard deviation when the daily product is averaged over 16-day MODIS periods (Fig. 7a and Fig. 7c).Table 3 suggests  poral variability is identical to MODIS variability on average, but Fig. 7 shows there are locations, particularly within the west coast ablation area, where MODIS variability is considerably higher.MAR v3.2 albedo variability in low elevation areas reaches a maximum of 0.09, while MODIS variability for the same region is 0.15 at maximum.At a daily temporal resolution, MOD10A1 daily variability in the ablation area (0.17 maximum, 0.07 on average) is considerably larger than the variability of MAR v3.2 albedo (0.12 maximum, 0.04 on average).As will be discussed in Sect.4.2, this may be the result of a positive bias in bare ice albedo from MAR, but may also be associated with errors introduced by cloud artifacts in the MOD10A1 product.For the accumulation area, the standard deviation of albedo for MAR and MODIS generally falls within the 16-day uncertainty of 0.04 for MCD43A3 high-quality albedo and daily uncertainty of 0.067 for MOD10A1 albedo estimated by Stroeve et al. (2005Stroeve et al. ( , 2006)).This limits the comparison among MAR and the MODIS products for high elevations.
For areas south of 70 • N and in the ablation area north of 70 • N, the two MODIS products are highly correlated (for MCD43A3 16-day periods, r 2 > 0.5), but in the accumulation area north of 70 • N this correlation decreases (Fig. 8a).Poor correlation in this area is likely a result of the low standard deviation of albedo which falls within the uncertainty range for MODIS.Maps of the coefficient of determination between MAR and MODIS (Fig. 8b and Fig. 8c) indicate that MAR v3.2 captures more than 50 % of the ablation area variability detected by satellite products for 16-day periods and more than 25 % for daily periods.It is, however, important to note that the daily variability from MOD10A1 is partially driven by cloud artifacts retained in the MOD10A1 product (Box et al., 2012).Again, in the accumulation area, it is difficult to draw any conclusions regarding correlation, as the variability in albedo is smaller than the assumed uncertainty for the MODIS products.

Albedo spatio-temporal variability
Further insights into the consistency of spatio-temporal variations in albedo between MODIS products and between MAR and MODIS products can be drawn from scatter plots for all MCD43A3 vs. MOD10A1 2000-2013 JJA albedo values (Fig. 9a) and MAR vs. MODIS values (Fig. 9b-d).Figure 9a indicates that MCD43A3 albedo is lower (by 0.03 on average) compared to MOD10A1 albedo, consistent with the significant difference between the products at high latitudes seen in Fig. 3a.There is a fairly good correlation between MCD43A3 and MOD10A1 (r 2 = 0.66) and the slope of the best linear fit (0.83) is close to 1.
When MAR is compared with MCD43A3 and MOD10A1 over 16-day periods (Fig. 9b and Fig. 9c), the correlation between MAR and satellite data is as good or better than the correlation between MOD10A1 and MCD43A3 (r 2 = 0.66 vs. MOD10A1 and 0.81 vs. MCD43A3).However, there is less agreement about the 1 : 1 line; a linear fit reveals a slope of 0.58 for MAR vs. MCD43A3 and 0.51 for MAR vs. MOD10A1.MAR overestimates low values of albedo (below 0.6) relative to satellite data, which is consistent with the apparent positive MAR bias in the ablation area seen in Fig. 3b and Fig. 3c.On a daily basis, there is a poor agreement between MAR and MOD10A1 (Fig. 9d, r 2 = 0.35), consistent with the poor correlations observed in Fig. 8d.Note that MOD10A1 albedo is only accurate to two decimal places, resulting in the apparent vertical lines in Fig. 8d.
Scatter plots of 2000-2012 JJA albedo values for both satellite products and MAR v3.2 vs. all weather station measurements (Fig. 10) indicate a strong correlation between in situ data and the two satellite products over 16-day periods (Fig. 10a and Fig. 10b; r 2 = 0.80 for MOD10A1, r 2 = 0.81 for MCD43A3), as well as a good agreement about the 1 : 1 line (slope = 0.95 for MOD10A1 and 0.88 for MCD43A3).MAR agrees reasonably well with in situ data, but the correlation is lower (r 2 = 0.78), and the slope (0.66) is further The  from 1. Again, it appears that MAR also overestimates low albedo values relative to in situ measurements, in consistency with Fig. 9b and Fig. 9c.On a daily basis, MOD10A1 albedo exhibits a nearly 1 : 1 relationship with daily in situ albedo (Fig. 10d; slope = 0.99), although there is increased scatter (r 2 = 0.75) due to higher variability on daily timescales (Fig. 7).Similarly, when MAR is compared with daily in situ measurements, the correlation is lower relative to the 16-day comparison (r 2 = 0.74), while the slope of the best fit line does not change substantially (slope = 0.65).
In situ and satellite data and MAR v3.2 outputs all indicate that spatio-temporal variability of albedo is higher in the ablation area (where the standard deviation of albedo is ∼ 0.13) than in the accumulation area (standard deviation of ∼ 0.04).This is to be expected, given that the ablation area undergoes a substantial seasonal cycle in melting.

MAR v3.2 vs. MAR v2.0 albedo
In order to further examine some of the discrepancies between MAR and observations, it was useful to examine differences between MAR v3.2 and MAR v2.0.MAR v2.0 has been validated against satellite and in situ data (e.g.Fettweis et al., 2005Fettweis et al., , 2011) ) and used for making future projections (Fettweis et al., 2013b;Tedesco and Fettweis, 2012).A major difference between MAR v3.2 and MAR v2.0 is in the scheme for calculating the albedo of bare ice; MAR v2.0 bare ice albedo is set to 0.45, while in MAR v3.2 it ranges between 0.45 and 0.55 as a function of surface melt (Table 1).
Scatter plots for MAR vs. best fit curves of the distribution (Fig. 11), suggest there is a bimodal distribution of ablation area albedo, which we attribute to the presence of two main surface types: ice (and firn) and snow.Pixels classified by MAR as having bare ice (or firn, surface density > 830 kg m −3 ) for at least 8 days of each 16-day period coincide with one of the peaks in the bimodal distributions (Fig. 11).However, there are differences in the observed distributions.MAR v2.0 exhibits a clustering of albedo values above 0.65 and below 0.55 (Fig. 11a).MCD43A3 exhibits an overlap in the distribution of the two modes, and there is a wider range of low albedo values (σ = 0.10 for MCD43A3 and 0.05 for MAR for the best fit of the lower albedo peak; Table 5).The MAR v3.2 distribution exhibits a slightly wider range of low albedo values (σ = 0.06 for the low albedo peak) with a mean that is positively shifted relative to MAR v2.0 (µ = 0.61 vs. 0.50; Fig. 11b).MOD10A1 does not appear to exhibit a bimodal distribution with two distinct peaks, but the best-fit curve agrees qualitatively with the observed distribution (Fig. 11c).The higher uncertainty and, there-fore, increased variability for the MOD10A1 product (Fig. 6; Stroeve et al., 2006) may possibly mask the two peaks of the distribution.Indeed, the best-fit bimodal distribution from MOD10A1 has a higher standard deviation of albedo for the higher albedo peak (σ = 0.06 for MOD10A1 vs. 0.04 for MCD43A3; Table 5).
We compared MAR v2.0 mean 2000-2012 clear-sky JJA albedo with albedo from MAR v3.2 and MODIS in Fig. 12ac.MAR v3.2 albedo is significantly larger in the ablation area compared with MAR v2.0 (Fig. 12a).Rather than being positively biased relative to MODIS (as is the case for MAR v3.2; Fig. 3), MAR v2.0 albedo is either negatively biased or is not significantly different from MODIS data (Fig. 12b and Fig. 12c).The difference in albedo scheme is the major difference between MAR v3.2 and MAR v2.0, and it results in a significant difference in SMB (Fig. 12d).The average ablation area JJA SMB for MAR v3.2 is higher by 0.53 mWE yr −1 compared with the average for MAR v2.0, a considerable fraction (roughly 25 %) of the mean ablation area JJA SMB from MAR v3.2, which is on average

Greenland Ice Sheet albedo trends
MAR v3.2, MCD43A3, and MOD10A1 consistently agree that there has been a significant decrease in albedo within the ablation area over 2000-2013, and that the largest decreases in albedo have occurred below 2000 m a.s.l.(Fig. 13 and Fig. 14).MCD43A3 shows a decrease of up to −0.1 per decade for pixels in the ablation area, as does MOD10A1 (both products show a decrease of −0.06 per decade for the entire area).MAR v3.2 agrees with these trends, but the overall magnitude is smaller (−0.03 per decade for the entire area).Within the accumulation area, MAR v3.2 disagrees with the two MODIS products as to the direction and magnitude of trends.MCD43A3 shows a decrease of −0.03 per decade on average, and MOD10A1 trends are somewhat larger (−0.04 per decade on average), while for MAR v3.2, trends are generally not statistically significant at the 95 % level for grid boxes above 2500 m a.s.l., and are slightly positive in some high-elevation areas.
For locations within the GrIS ablation area, trends at GC-Net stations with a record of at least nine years are consistent with significant decreases in albedo, indicated by MODIS and MAR for the periods covered (2000-2012or 2004-2012;Table 6).The magnitude of the trends varies between MAR v3.2, MODIS, and in situ data at individual stations.These differences can be attributed in part to the high spatiotemporal variability of albedo within the ablation area.This can potentially lead to trends at a weather station that are substantially different from trends within a 500 m MODIS grid box containing the location of that weather station.At higher elevations, this factor is less important as there is less spatiotemporal variability in albedo (Fig. 9 and Fig. 10  ) for the same period.Note that in the ablation area, where net SMB is negative (Fig. 1), a positive SMB bias indicates a net mass loss that is reduced in magnitude.Grid boxes where differences are not significant at the 95 % confidence level are marked with a black "x".6).

Albedo properties of the GrIS common to all observations and model results
The results presented above highlight certain features of GrIS albedo variability that are common to in situ, satellite data, and model results.MAR, MODIS, and in situ data capture general spatial patterns of low albedo in the ablation area, which increases with increasing elevation below ∼ 2000 m a.s.l and is relatively insensitive to higher elevations (Fig. 2 and Fig. 4a, Table 4).This spatial variability is consistent with the presence of meltwater and bare ice exposure at low elevations, which are a function of surface air temperatures, and therefore elevation (Tedesco et al., 2011;Fettweis et al., 2011).Bare ice in the ablation area is often covered with dust, further reducing low elevation albedo (Bøggild et al., 2010;Wientjes and Oerlemans, 2010;Wientjes et al., 2010).At high-elevation areas that are permanently snow covered, particularly at northern sites, melting is infrequent (Nghiem et al., 2012;Tedesco et al., 2011) and albedo variability is primarily associated with accumulation, subsequent dry snow grain size metamorphism (Wiscombe and Warren, 1980), and possibly impurities (Dumont et al., 2014).Low elevation melting and bare ice exposure during warm summer months reduce surface albedo relative to snow albedo, resulting in a seasonal cycle that increases local variability.As Fig. 7 and Tables 3 and 4 show, there is higher variability in ablation area albedo (where the mean standard deviation of albedo at in situ stations ranges between ±0.06 and ±0.09) relative to the accumulation area (where standard deviations range between ±0.02 and ±0.04).
As noted in Sect.3.5, MAR, MODIS, and in situ data agree there has been a significant decline in ablation area albedo between 2000 and 2013.These trends in surface albedo are associated with increased melting and bare ice exposure resulting in a decline in ablation area SMB, captured by models (Fettweis et al., 2011;Ettema et al., 2009) and in situ observations (van de Wal et al., 2012).Increased melting has been linked to higher regional atmospheric air temperatures, asso-  (Fettweis et al., 2013a;Häkkinen et al., 2014).

Variation of albedo with latitude
Results from Sect.3.1 indicate that above 70 • N, MOD10A1 shows an increase in albedo with latitude, MCD43A3 exhibits a decrease, MAR shows little change, and a small increase with latitude at local weather stations is noted (Figs. 2, 4b, 5 and Table 4).The increase with latitude at local stations is likely unaffected by differences in spectral range between MODIS and in situ sensors (Fig. 6).
Theoretically, snow albedo is expected to increase with increasing solar zenith angle, particularly for high solar zenith angles (Wiscombe and Warren, 1980), and, therefore, will increase slightly at high latitudes, as long as other factors do not contribute to lower albedo values.Wang and Zender (2009) compared 16-day MCD43C3 albedo with GC-Net measurements and suggest that the MCD43C3 product is unrealistic at higher latitudes, in particular for solar zenith angles > 55 • (the MCD43C3 product differs from the MCD43A3 product used here only in its grid).Schaaf et al. (2011) and Stroeve et al. (2013) suggest that the findings of Wang and Zender (2009) are inaccurate, partially because they did not separate results for high-vs.lowquality albedo.We have considered this in our study: results for all MCD43A3 data are shown along with good quality MCD43A3 data in Fig. 4b.While the use of only good quality data increases MCD43A3 albedo above 70 • N, it does not fundamentally change the dependency of MCD43A3 albedo on latitude.For MOD10A1, excluding low quality data has little effect on the binned values.
It should also be noted that the MOD10A1 product may be positively biased above 70 • N, given that it is comparable with uncorrected GC-Net data, which are likely to be positively biased (Fig. 5).We do not have a reasonable explanation for this potential bias, but as noted by Box et al. (2012), the MOD10A1 product contains artifacts that have not been removed during quality control, even for "good quality" data.The in situ observations of Konzelmann and Omura (1995) also suggest that values of albedo above 0.84 are unrealistic for snow under clear-sky conditions.Part of the reason for discrepancies in the latitudinal dependence of albedo may be associated with biases resulting from viewing geometry or sun angle, which vary with latitude, making it difficult to draw conclusions from the various observational data sets as to "true" variations in albedo with latitude.

Differences between MAR v3.2, MAR v2.0, and observed albedo
The major difference between MAR v3.2 albedo and observed albedo is an overall positive bias in the ablation area.
This bias can be seen most clearly as a difference of ∼ 0.1 between MAR v3.2 and the two MODIS products along the west coast ablation area in Fig. 3b and Fig. 3c, and in a difference between MAR v3.2 and both MODIS products of 0.06 at in situ stations (with low elevation stations mostly located in the west coast ablation area).Mean ablation area albedo from local stations is also comparable with coincident MAR v3.2 albedo (Table 4), but local station measurements are likely positively biased, further confirming a positive MAR bias in this area.Scatter plots of ablation area albedo appear to confirm this: when MAR v3.2 is compared with both MODIS data and in situ measurements (Figs.9b, 9c and 10c) the result is a best fit line with a slope smaller than one.Additionally in the same area where MAR v3.2 appears positively biased in the west coast ablation area, MODIS exhibits relatively high variability compared with MAR v3.2 (as discussed in Sect.3.2; Fig. 7).
Biases in MAR ablation area albedo are related to its ability to capture the observed bimodal distribution in ablation zone albedo (Fig. 11) associated with two main surface types, ice and snow.The positive bias from MAR v3.2, as well as the relatively low modelled variability in the ablation area is the result of the albedo values set for bare ice in MAR v3.2 (Table 1) that may be too high on average.MAR v2.0 albedo, by contrast, which has a fixed bare ice albedo of 0.45, generally exhibits a negative bias in most portions of the ablation zone.A bare ice albedo that is too high will also lead to a smaller difference between the albedo values of melting snow and bare ice, reducing temporal variations in ablation area albedo, resulting in the relatively low variability from MAR v3.2 (Fig. 7).
An examination of Fig. 11 indicates that the low albedo peak for MAR v3.2 is closer to being normally distributed compared with the peak for MAR v2.0, and is, therefore, a better match to the distribution from MCD43A3.However, MAR v3.2 overestimates bare ice albedo, as already discussed, and still does not fully capture variability in the low albedo peak for MODIS albedo (σ = 0.06 for MAR v3.2 and σ = 0.10 for both MODIS products).Although MAR v3.2 appears to correct a low albedo bias present in MAR v2.0, and introduces a somewhat more realistic distribution of albedo in the ablation area, the results suggest it also introduces a positive albedo bias, particularly along the west coast ablation zone, where impurities are numerous.
MAR albedo is only a function of accumulated meltwater and does not explicitly take into account the presence of dust, surface lakes and surface streams, including the West Greenland "dark zone" (van de Wal and Oerlemans, 1994;Wientjes and Oerlemans, 2010), which reduces bare ice albedo and likely introduces increased ablation area albedo variability.Assigning a wider range of MAR albedo values for bare ice (which has been implemented in the most recent release of MAR, v3.4) may improve its representation of the distribution of bare ice albedo, but may not necessarily improve its ability to capture the spatial distribution of ablation area albedo.This could potentially be achieved through the inclusion of an explicit representation of dust and sub-grid-scale hydrology in the model.

Discrepancies in accumulation area trends
As noted in Sect.3.5, there is a discrepancy between the satellite products, in situ data, and model results regarding albedo trends in the accumulation area of the ice sheet.MOD10A1 and MCD43A3 show significant decreases in accumulation area albedo (−0.04 to −0.03 per decade) while MAR v3.2 trends are generally not statistically significant, and in situ trends are generally small (not larger than −0.01 per decade) or not significant.A possible explanation for this discrepancy is that MODIS trends are negatively biased as a result of declining instrument sensitivity of the MODIS sensors (Wang et al., 2012).In particular, a larger degradation has been observed for the MODIS Terra satellite (Wang et al., 2012).The MCD43A3 product uses data from both the Terra and Aqua satellites, while MOD10A1 only uses data from Terra.This could potentially explain the larger trends for MOD10A1 relative to MCD43A3 (Table 6, Fig. 13).Box et al. (2012) conclude that declining instrument sensitivity does not substantially affect GrIS albedo trends, as larger trends are found in GC-Net data relative to MOD10A1 for 70 % of cases where trends are deemed to be significant.In contrast to the findings of Box et al. (2012), we did not find JJA GC-Net trends larger than those of MODIS, except in some instances within the ablation area, where there is high local variability in surface properties.The analysis performed here differs from that employed by Box et al. (2012).Differences in trends found in this study may have resulted from a focus on trends for the entire JJA period rather than on monthly trends, and calculated trends for 16-day albedo values rather than calculated monthly albedo from integrated fluxes over a 1-month period, as was done by Box et al. (2012).We also investigated the possibility that the smaller spectral interval of GC-Net data influences trends by comparing MCD43A3 visible vs. SW albedo trends, but did not find the trends to be significantly different from each other.We are not able to confirm that the larger trends from MODIS are associated with declining instrument sensitivity, as this analysis is outside the scope of this study.However, the findings of this study seem consistent with this possibility and it is suggested as a topic for future research.

Conclusions
We have examined spatio-temporal variability and trends in GrIS albedo using in situ measurements, satellite products obtained from MODIS data, and outputs of two versions of the MAR RCM.The results presented here reveal areas of agreement as well as discrepancies between observational The results presented here show that albedo varies spatially as a function primarily of surface properties, in particular melting and bare ice exposure in the ablation area.These factors are also associated with temporal variations in albedo, resulting in high variability in low elevation regions.The differences in variations with latitude indicated by satellite products appear likely to be a function of inaccuracies associated with the products themselves, rather than a record of actual variations in surface albedo, particularly as the two products are derived from the same MODIS sensors.
Both satellite products and MAR model outputs (for v2.0 and v3.2) suggest a bimodal distribution of surface albedo within the ablation area of the ice sheet.Based on model results, we infer that this distribution is associated with the presence of two primary surface types within the ablation area, snow and bare ice.The model's inability to capture the full range of low elevation albedo leads to inaccuracies in the representation of spatio-temporal variations in albedo, which can substantially impact the representation of SMB.The MAR version examined here (v3.2) appears to better represent the full range of bare ice albedo in the ablation area relative to a previous version (v2.0), but a lower minimum bare ice albedo value (as is implemented in the next version of MAR, v3.4), may produce results that are more consistent with observations.Even so, it may be necessary to account for the presence of impurities and sub-grid-scale hydrology in order to fully capture spatial variations in albedo.
The analysis performed here indicates a statistically significant decrease in ablation area albedo over the period 2000-2013 and is consistent with previous studies (Box et al., 2012;Tedesco et al., 2011Tedesco et al., , 2013;;Stroeve et al., 2013).This decrease is consistent with a coincident decline in ablation area SMB recorded by both models and observations (e.g.Fettweis et al., 2011;Ettema et al., 2009).Our results are inconclusive regarding high elevation trends in albedo; we observe inconsistencies between satellite-derived trends and trends obtained from in situ measurements and MAR v3.2 results.We are therefore unable to confirm previously reported decreases in surface albedo at high elevations.
Future research should be directed towards understanding the reasons for discrepancies between satellite products, in situ data and model results, in order to better understand changes in GrIS albedo.This includes resolving discrepancies regarding high-elevation trends, and discrepancies in mean satellite-derived surface albedo at high latitudes.Models such as MAR appear to be effective at capturing surface albedo, but refinements are necessary for representation of surface albedo in low elevation areas.In particular, the representation of bare ice albedo is critical.Sensitivity studies, such as those performed by van Angelen et al. (2012) of the impact of surface albedo on SMB variability, may help to quantify the accuracy with which surface albedo must be modelled for a given region.Analysis of spatio-temporal variations in albedo across different spatial scales (including at a higher spatial resolution than has been examined here) may also become increasingly important as models operate at higher spatial resolutions, and as we seek to understand the GrIS surface mass and energy budget in greater detail.Given the strong relationship between surface albedo and SMB, future studies are crucial for efforts aimed at estimating and predicting the impact of current and future climate change on GrIS SMB.

Figure 1 .
Figure 1.MAR v3.2 mean September 2000-August 2013 SMB (mWE yr −1 ) and locations of all GC-Net and K-Transect weather stations.Pixels not defined as 100 % ice covered in MAR v3.2 are masked out.The dotted black line is the mean equilibrium line (where mean SMB is 0).The K-transect stations (S5, S6, S9, S10) are red, while GC-Net stations are black.Stations in grey are GC-Net stations not used in this study.Other contour lines indicate elevation in metres above sea level.The inset shows individual stations near the west coast ablation zone.

Figure 2 .
Figure 2. Mean 2000-2013 JJA albedo (unitless) for (a) the MCD43A3 BSA SW product (on the MAR grid) (b) MOD10A1 product (on the MAR grid), and (c) MAR v3.2 clear-sky albedo.Only good quality data MODIS data are used here.

Figure 3 .
Figure 3. Mean difference in JJA albedo (unitless) for the 2000-2013 period: (a) MCD43A3 BSA SW minus MOD10A1 (b) MAR v3.2 clear-sky minus MOD10A1, and (c) MAR v3.2 clear-sky minus MCD43A3.In each case, only coincident data for each of the two data sets being compared is used.MAR grid boxes where the difference is not statistically significant at the 95 % confidence level are marked with a grey "x".

Figure 4 .
Figure 4. (a) Mean 2000-2013 JJA MOD10A1, MCD43A3 BSA SW, and MAR v3.2 clear-sky GrIS albedo (unitless) as a function of elevation divided into 150 m elevation bands.Error bars indicate standard deviation within each elevation band.(b) The same as (a) but for albedo as a function of latitude, divided into 2 • Latitude bands."Good qual."indicates results obtained by only using "good quality" MODIS data."All qual."indicates that all available MODIS observations have been used.

Figure 5 .
Figure 5. 2000-2012 mean JJA albedo (unitless) for the MAR accumulation zone vs. latitude, for MOD10A1, MCD43A3 BSA SW, MAR v3.2 clear-sky, and GC-Net station data (black circles) for stations with a record spanning at least seven years of the 2000-2012 period.Only MODIS data flagged as "good quality" were used here.The error bars for GC-Net stations indicate the range of corrections to GC-Net data (between 0.04 and 0.09) employed by Stroeve et al. (2005).

Figure 11 .Figure 11 .
Figure 11.Scatter plots and histograms for JJA 2000-2012 albedo [unitless] within the MAR 4 v3.2-definedGrIS ablation zone, for (a) MAR v2.0 clear sky (16 day avg.) vs. MCD43A3 5 BSA shortwave.(b) The same as (a), but for MAR v3.2.(c) The same as (a) but for 6 MOD10A1 albedo (averaged to 16 day periods).(d) The same as (c) but for MAR v3.2.7Points where there is snow or firn (surface snowpack density > 830 kg/m 3 ) for more than 8 8 days of a 16 day period are shown in red.Light blue curves show the best fit to each 9 distribution obtained using maximum likelihood estimation.10 11

Figure 12 .
Figure 12.(a) MAR v3.2 clear-sky minus MAR v2.0 clear-sky mean JJA albedo (b) MAR clear-sky v2.0 minus MOD10A1 2000-2012 mean JJA albedo, MAR clear-sky v2.0 minus MCD43A3 BSA SW 2000-2012 mean JJA albedo, and (d) MAR v3.2 minus MAR v2.0 mean JJA SMB (mWE yr −1) for the same period.Note that in the ablation area, where net SMB is negative (Fig.1), a positive SMB bias indicates a net mass loss that is reduced in magnitude.Grid boxes where differences are not significant at the 95 % confidence level are marked with a black "x".

Figure 13 .
Figure 13.JJA mean albedo trends (2000-2013) in units of fraction per decade for (a) MCD43A3 BSA SW albedo, (b) MOD10A1 albedo, and (c) MAR clear-sky albedo.Grid boxes where trends are not significant at the 95 % confidence level are marked with a black "x".

Figure 14 .
Figure 14.Mean annual JJA ice sheet albedo (solid lines) simulated by MAR v3.2 (clear-sky; blue), MOD10A1 (black) and MCD43A3 BSA SW (orange) for 2000-2013 and best linear fit (dashed lines) for (a) the entire ice sheet, (b) the accumulation zone, and (c) the ablation zone defined with MAR v3.2.The trends shown are statistically significant at the 95 % confidence level for the MODIS products, but are not statistically significant for MAR.Shaded areas show annual JJA standard deviation of albedo for 16-day periods from each data set.Note that the y axis interval is the same for all graphs, but is shifted by 0.1 for (c).

P
. M. Alexander et al.: Variability and trends in Greenland Ice Sheet albedo and model estimates of GrIS albedo spatio-temporal variability.Examining local measurements, satellite data, and model results concurrently reveals information about the GrIS albedo and potential biases that would not be revealed by examining observational data sets or model results individually.

Table 1 .
Range of possible albedo values for different surface types in MAR v2.0 and v3.2

Table 2 .
GC-Net and K-Transect weather stations used in this study and years of coverage.

Alexander et al.: Variability and trends in Greenland Ice Sheet albedo grid
box classified as less than 100 % ice-covered in MAR v3.2, we compared in situ data to MAR v3.2 data from the ice-covered sector of that grid box rather than data from the entire grid box.

Table 3 .
Mean 2000-2012 JJA GrIS albedo, for MOD10A1, MCD43A3 BSA SW, and MAR clear-sky albedo, averaged within the mass balance areas shown in Fig.1and Table2.Only good quality MODIS data are used here.All data have been averaged over the same 16-day period of the MCD43A3 product.Only periods when coincident data for all data sets were available have been included.

Table 4 .
Same as Table3, but for the average of 16-day data for all in situ stations within each region.

Table 6 .
Trends (and 95 % confidence intervals) in JJA albedo (fraction per decade) at GC-Net and K-Transect weather stations and the nearest MOD10A1, MCD43A3, and MAR pixels.In this case, MODIS data flagged as "other quality" have been included.Only 16 day periods when coincident estimates are available for all data sets have been used.Values in bold indicate trends significant at the 95 % confidence level.