MODIS snow albedo bias at high solar zenith angles relative to theory and to in situ observations in Greenland Remote Sensing of Environment

In situ measurements of snow albedo at ﬁ ve stations along a north – south transect in the dry-snow facies of the interior of Greenland follow the theoretically expected dependence of snow albedo with solar zenith angle (SZA). Greenland Climate Network (GC-Net) measurements from 1997 through 2007 exhibit the trend of modest surface brightening with increasing SZA on both diurnal and seasonal timescales. SZA explains up to 50% of seasonal albedo variability. The two other environmental factors considered, temperature and cloudiness, play much less signi ﬁ cant roles in seasonal albedo variability at the ﬁ ve stations studied. Compared to the 10-year record of these GC-Net measurements, the ﬁ ve-year record of MODIS satellite- retrieved snow albedo shows a systematic negative bias for SZA larger than about 55°. Larger bias of MODIS snow albedo exists at more northerly stations. MODIS albedos successfully capture the snow albedo dependence on SZA and have a relatively good agreement with GC-Net measurements for SZA<55°. The discrepancy of MODIS albedo with in situ albedo and with theory is determined mainly by two related factors, SZA and retrieval quality. While the spatiotemporal structure, especially zonal features, of the MODIS-retrieved albedo may be correct for large SZA, the accuracy deteriorates for SZA>55° and often becomes physically unrealistic for SZA>65°. This unphysical behavior biases parameterizations of surface albedo and restricts the range of usefulness of the MODIS albedo products. Seasonal-to-interannual trends in surface brightness in Greenland, and in polar (i.e., large SZA) regions in general, and model simulations of these trends, should be evaluated in light of these limitations.


Introduction
Albedo is the fraction of incident solar energy reflected by the surface over the entire solar spectrum (Dickinson, 1983). The land surface albedo is one of the key parameters and a driver in climate and weather models since it regulates the shortwave radiation absorbed by the Earth surface (Lucht, 1998;Wang et al., 2005). Snow or ice covers about 40% of northern hemisphere land in winter and causes the greatest seasonal surface albedo variability (Qu & Hall, 2005). Snow and ice cover are important components of the Earth's energy balance because of their high reflectance in the visible bands and great seasonal variation. Fresh snow reflects more than 80% of incident solar energy. Snow albedo decreases as grain size increases (due to aging and/or warming) and due to accumulating impurities. The snowalbedo feedback amplifies the sensitivity of snow and ice to small changes in albedo (Stroeve et al., 2005;Flanner & Zender, 2005;Flanner et al., 2007). Accurate determination of snow albedo is therefore critical for understanding and predicting cryospheric climate sensitivity.
Snow bi-directional reflectance varies strongly with solar zenith angle (SZA) and viewing geometry Jin et al., 2003a,b;Salomon et al., 2006), and snow directionalhemispherical reflectance (DHR) has a far greater magnitude of increase with SZA in infrared wavelengths (1.03 µm) than in visible (0.55 µm) wavelengths (Schaepman-Strub et al., 2006). The hemispherically integrated flux reflectance, or spectral albedo, integrates this angular variation at a given wavelength. The scalar of greatest interest to climate studies is the spectrally integrated broadband solar albedo which describes the net flux reflectance of the land surface. Near-global maps of surface albedo retrieved from satellite measurements are available from multiple instruments that include the Advanced Very High Resolution Radiometer (AVHRR), the Earth Radiation Budget Experiment (ERBE), Polarization and Directionality of the Earth's Reflectances (POLDER), Clouds and the Earth's Radiant Energy System (CERES), and the MODerate resolution Imaging Spectroradiometer (MODIS) (Schaaf et al., 2008a). Accurate MODIS snow albedos are of particular value to climate research because the two platforms (Terra and Aqua) sample the cryosphere, which is rapidly changing, at high temporal and spatial resolutions relative to the other platforms.
Greenland's snow and ice play pivotal roles in the global climate because of their year-long high reflectivity (or albedo), large area and substantial volume of fresh water stored (Steffen & Box, 2001). The high SZA poses serious challenges not only to in situ solar radiation measurements (Augustine et al., 2000), but also to the consistency and accuracy of the MODIS radiative transfer models applied to polar regions (Lucht, 1998;Liu et al., 2009). Assuming that snow surface structure and conditions (related to snow grain size, age, contamination, etc.) do not change, theory shows that snow surface albedo increases with SZA ( Fig. 1) because the increased path over which obliquely-incident photons interact with snow grains allows more multiple scattering and less penetration or absorption by the snow surface. This results in a larger fraction of solar radiation reflected by snow Wang et al., 2005;Lucht, 1998). Several field studies show albedo increases with SZA not only over snow, but also over desert, vegetated and ocean surfaces as well Jin et al., 2003b;Painter & Dozier, 2004;Wang et al., 2005;Liu et al., 2009).
The existence of a MODIS snow albedo bias at high SZA is also noted by Hall et al. (2009a). We will show that this bias leads to substantial underestimates of climatological albedo. The standard MODIS BRDF/Albedo algorithm estimates the albedo at local noon (Lucht, 1998;Schaaf et al., 2002;Jin et al., 2003a). Jin et al. (2003b) show (1) that in situ albedo measurements do follow, and that MODIS standard albedo retrievals do not follow, the expected SZA-albedo relation at large zenith angles. The overall accuracy of MODIS albedo is within 0.05 as compared to the ground observations at the Oklahoma CART site and has an increasing negative bias when SZA > 70° (Liu et al., 2009). High quality MODIS snow albedo retrievals can be obtained over the homogeneous snow surfaces on Greenland, usually for SZA < 55-60°, and larger errors exist at high SZA (Stroeve et al., 2005). These large errors may bias parameterizations of surface albedo and restrict the range of usefulness of the MODIS product. Moreover, these errors also undermine measurement-based evaluation of climate model albedo estimates over Greenland, Antarctica, and other high latitude snow-covered regions (Zhou et al., 2003;Oleson et al., 2003).
This study therefore aims to characterize, understand, and address the influence of SZA on both in situ and remotely sensed snow albedo and their discrepancy. Our strategy adopts the in situ GC-Net solar radiation observations at the perennial snow-covered stations on Greenland to evaluate the MODIS albedo product. First we analyze diurnal and seasonal cycles of albedo and solar radiation. We then identify the SZA dependence of, and air temperature impacts on snow albedo using the in situ measurements. Finally we compare the MODIS albedo with in situ measurements and quantify the bias of MODIS snow albedo at high SZA from several aspects.

MODIS albedo products
MODIS is used in two distinct sets of 16-day BRDF/albedo products. One is the Terra-only 1 km MOD43B in an Integrated Sinusoidal Grid (ISG) projection with standard tiles representing 1200 × 1200 gridpoints and 0.05°MOD43C in the Climate Modeling Grid (CMG) (Schaaf et al., 2002). The other is the Terra and Aqua MODIS combined products (MCD), which includes 500 m MCD43A and 1 km MCD43B in ISG projection, and the global 0.05°MCD43C in the CMG projection. Both Terra and Aqua data are used to generate this product. The combination provides the highest probability for quality input data, and increases full retrievals across the globe by 50% (Salomon et al., 2006). The 500 m product is first produced at the latest version 5 in the combined product (Schaaf et al., 2008b). Version-5 MODIS/Terra +Aqua BRDF/Albedo products are validated stage 1 and its accuracy has been estimated using a small number of independent measurements obtained from selected locations and time periods and groundtruth/field program efforts. In addition, there is another daily albedo product developed by Klein and Stroeve (2002) and Stroeve et al. (2006), discussed in Hall et al. (2009a), and used over the Greenland Ice Sheet by Hall et al. (2009b).
Each of the five products (two Terra MODIS-only and three Terra and Aqua MODIS combined) consists of four components. The first two components are seven spectral bands (MODIS bands 1-7, e.g., MCD43C1) and three broadbands' (0.3-0.7, 0.7-5.0, and 0.3-5.0 µm, e.g., MCD43C2) BRDF model parameters from which users can reconstruct the entire BRDF and compute the directional reflectance at any view or desired sun zenith angle. Thus, the directional hemispheric reflectance (black-sky albedo, BSA) at any angle and further the bihemispheric reflectance (white-sky albedo, WSA) can be generated. For consistency, we use the reflectance terminology employed in much, but not all, of the relevant MODIS literature (Schaaf et al., 2002;Stroeve et al., 2005), such as MODIS BRDF/albedo product, white-sky albedo and black-sky albedo. Here BRDF, whitesky albedo and black-sky albedo respectively represent the bidirectional reflectance distribution function (BRDF; Case 1), directional hemispheric reflectance (DHR; Case 3) and bihemispheric reflectance (BHR; Case 9) as documented in Schaepman-Strub et al. (2006). The third component (e.g., MCD43C3) is the standard suite of albedo values for black-sky and white-sky of the seven MODIS spectral bands and the three broadbands. The black-sky albedos are reported for the local noon SZA at each gridpoint. The final component (e.g., MCD43C4) is the nadir viewing BRDF-Adjusted Reflectance (NBAR)  Zender (1999) and Flanner and Zender (2006). for the seven spectral bands (Schaaf et al., 2002). Besides the standard suite of albedo values in MCD43C3, it also contains the snow cover fraction (SCF), local noon solar zenith angle (SZA), and BRDF quality code which ranges from 0 (the best, above 75% of full retrievals) to 4 (the worst, backup retrievals).
Following Wang et al. (2005), this study uses the global 0.05°s tandard suite of albedo values (MCD43C3). It is produced every 8 days within 16 days of acquisition. For example, production of period 2002001 includes acquisition from days 001 to 016, production period 2002009 includes acquisition from days 009 to 024. In the 16 days of sequential multi-angle observations after cloud screening and atmospheric corrections, if the majority of observations are recorded as snow-covered, then the algorithm uses only snowcovered observations for the parameter retrievals; otherwise, the algorithms conservatively use snow-free observations for parameter retrievals (Schaaf et al., 2002;Salomon et al., 2006). Meanwhile, if at least seven cloud-free observations (there are several observations per day in Greenland, cloud-cover permitting) during the 16 days, a full model inversion is attempted; otherwise, a magnitude inversion or a back up algorithm is performed (Schaaf et al., 2002). The MODIS albedo retrieval strategy, including radiative transfer methods and bidirectional models, full model inversions or full retrievals, backup algorithms, a priori estimate and fill values, is documented in Lucht (1998), Lucht et al. (2000), Schaaf et al. (2002) and Salomon et al. (2006).

GC-Net measurements
The GC-Net was established in spring 1995 and consisted of 21 Automatic Weather Stations (AWS) by fall 2000 (when MODIS became operational) distributed over the Greenland ice sheet (Fig. 2). Each AWS is equipped with a number of meteorological and glaciological instruments to measure multiple parameters, such as air and snow layer temperature, humidity, wind speed, snow height, shortwave incoming and reflected radiation (Steffen & Box, 2001). Among the 21 stations, we select five stations ( Fig. 2 and Table 1) meeting this study's selection criteria: (1) air temperature is always below zero, and (2) they are 100% covered by (perennial) snow through an entire year. These criteria ensure with high confidence that wet processes such as snow-melt do not confound our study. The five stations selected are located along the crest of the ice sheet and cross Greenland from south to north (South Dome-SDM, Saddle-SDL, Summit-SMM, NGRIP-NGR, Homboldt GI.-HMG) at elevation ∼ 2000 m and above, while most of the melt supra-glacier lakes occur below 1600 m in the southwest of Greenland (Sundal et al., 2009) and there is no wet snow melting on this high ice sheet facies of Greenland even in the most extensive melting period of August in 2007 (Hall et al., 2009b). Restricting our study to dry snow regions ensures that snow-melt does not contribute to the discrepancy between in situ and satellite-retrieved albedos. In addition, there is no nearby vegetation and the study area is far away from desert, boreal forest and extensive human activities. The influences of soil type, vegetation cover/type, and contamination by aerosol (wild fire, dust storm) and/or other extensive human activity are minimized though it is impossible to completely exclude Flanner et al., 2007). Thus, our selected GC-Net stations are ideal for examining the influences of SZA on the in situ and MODIS snow albedo.
The in situ measurements used in this study are air temperature, and shortwave solar downwelling and upwelling irradiance. Air temperatures are sampled every 60 s by a Type-E Thermocouple at 2 m height and averaged over an hour. Shortwave fluxes are measured by pairs of LI-COR 200SZ pyranometers with spectral band width from 0.4 to 1.1 µm. The LI-COR 200SZ photodiode sensor does not cover the entire solar radiation range, but is factorycalibrated against an Eppley Precision Spectral Pyranometer (PSP) under natural daylight conditions to an equal spectral response from 0.28 to 2.80 µm, the standard optical black thermopile instruments' spectral sensitivity (LI-COR, 2005;Stroeve et al., 2005), which represents over 98% of solar radiation. The LI-COR 200SZ  Table 1 Greenland Climate Network (GC-Net) Automatic Weather Stations (AWS) used in this study (Steffen & Box, 2001 pyranometers are horizontally leveled to measure incident and reflected hemispheric radiant flux density (irradiance) and to provide hourly averages from 15-s sample intervals. The downwelling and upwelling shortwave radiative fluxes have 5% uncertainties (Box & Steffen, 2000;LI-COR, 2005). Systematic positive biases may exist in the reflected irradiance since snow has high reflectance in the 0.3-0.7 µm wavelength range and low reflectance for wavelengths larger than 1.1 µm and the LI-COR pyranometer is calibrated based on the downward solar radiation against Eppley PSP measurements. Thus, the GC-Net snow albedo data may have a positive bias especially under clear-sky condition because the snow surface depletes most of the solar radiation for wavelengths larger than 1.1 µm and the reflected irradiance might be overcalibrated. Under cloudy sky, the potential over-calibration for the upward reflected irradiance is negligible. Compared to Eppley PSP measurements at Swiss Camp in Greenland, the upward and downward LI-COR-measured shortwave irradiance errors did not exceed 2.7%. To compensate for the bias in the reflected LI-COR irradiances, Stroeve et al. (2005) apply a site-specific correction to GC-Net albedo in the form of an offset value ranging from 0.01 to 0.09. Because the snow albedo discrepancies against the precise pyranometer measurements vary greatly among different calibration sites reported in Stroeve et al. (2005) and our uncorrected GC-Net measured snow albedo values at the five stations that we consider fall in a similar range as those precise pyranometer measurements, and in order to avoid additional site-specific uncertainties due to correction to the GC-Net snow albedo, none of the downward and upward irradiance data are adjusted in this study. The potential positive bias in the GC-Net snow albedo must be borne in mind when comparing to the MODIS-retrieved albedo. Moreover, the relatively consistent GC-Net snow albedo uncertainties do not prevent us from assessing the overall MODIS snow albedo bias at the high SZAs in different seasons.

Spectral range difference between MODIS and GC-Net Albedo
The MODIS shortwave albedo is obtained from converting MODIS narrowband albedo for spectral channels (e.g. 1, 2, 3, 5, and 7 with actual wavelength from 0.46 to 2.15 µm) via a narrow-to-broadband (NTB) conversion coefficients (Liang et al., 1999;Stroeve et al., 2005) and represents spectral sensitivity from 0.3 to 5.0 µm. The LI-COR 200SZ pyranometer used in GC-Net has narrower spectral band width from 0.4 to 1.1 µm and is calibrated to equal spectral response from 0.28 to 2.8 µm. The spectral sensitivity difference between GC-Net and MODIS albedo is negligible (∼1.5% of downwelling solar irradiance and 0.01 snow albedo difference). We will examine both the absolute albedos at five GC-Net stations and their changes with high SZA. The albedo changes (especially their signs) are expected to be relatively insensitive to, and thus more robust than, the absolute albedo comparisons.

MODIS data
Five-year MCD43C3 data from 2003 to 2007 are used in this study. They originate in Hierarchical Data Format (HDF) on the global CMG at 0.05°resolution, and we convert them into the network Common Data Form (netCDF) for regridding and analysis. For zonal analysis on Greenland, the area-weighted mean albedo is computed for each 1°l atitude by 10°longitude (310°-320°E) from latitude 66°N to 82°N. At each latitude band, there is a total of 4000 gridpoints (20 × 200). Only those gridpoints for which the SCF is 100% are counted. The white-sky and black-sky albedos are similar on near homogeneous snow surfaces (Stroeve et al., 2005). Each can closely represent the bluesky albedo, which is proportionally combined from the white-sky and black-sky albedos. The shortwave albedo is inferred from the seven single band spectral albedo (Schaaf et al., 2002;Jin et al., 2003a,b;Wang et al., 2005). Therefore, only white-sky shortwave (0.3-5.0 µm) and one single band (band4, 0.545-0.565 µm) albedos are used as examples in this study although similar analyses are conducted and results found for other single band and broadband albedos.
A square of 3 × 3 MODIS gridpoints (∼17 km × 17cos(β) km, β = stations' latitude, from N63°to N78°) centered on each of the five stations is used to compare the MCD43C3 product with GC-Net observations. In order to reduce the effect of missing data, the albedo average of the 9 gridpoints is compared with the in situ measurements since the albedo difference between the central gridpoint and the average albedo of the 9 gridpoints is negligible at the near homogeneous snow surface. The point-scale in situ measurements are assumed to represent the ground truth of areal value over the 3 × 3 grid patch. Although the (assumed to be homogeneous) snow cover and relatively flat surface may minimize the point-areal representative discrepancy (Jin et al., 2003b), any conclusion about the accuracy of MODIS BRDF/albedo must consider the scale mis-match as well as the uncertainties of in situ shortwave downwelling and upwelling irradiance measurements.
Besides the MCD43C3 albedo product, the cloud fraction value contained in the daily snow cover product (MOD10C1) (Hall & Riggs, 2007) is also used to examine its seasonal pattern and to define the clear sky when in situ solar irradiance measurements are converted into a 16-day average.

GC-Net data
The calibrated data obtained from GC-Net are averaged hourly. We compare the MODIS surface albedo at local noon (Jin et al., 2003b) to the GC-Net estimated mean albedo from three hourly averages of shortwave downwelling and upwelling irradiance centered on local noon (11:00, 12:00 and 13:00) each day. Using the 3-hour (rather than 1-hour) average of irradiance around local noon to construct the in situ daily noon surface albedo smooths the albedo variation due to high frequency processes, such as wind blowing, fresh snow addition, and other disturbances. We also examine the GC-Net 24-hour average daily shortwave irradiance to discover its relationship with the 24-hour air temperature average. All irradiance data (SZA < 80°) are screened by a mandatory 10% outliers rule that albedo dropping beyond 10% against the mean of the two neighboring hourly data is replaced by the mean or just discarded for continuous non-credible data. MODIS cloud-cover value in the daily MODIS snow cover product (MCD10C1) is used to define a clear-sky day when daily cloud-cover fraction is less than 10%. All clear-sky daily values are converted into 16-day average in the same period as the MODIS albedo computation period. In situ albedo (daily and 16-day) is finally computed as the ratio of the mean upwelling shortwave irradiance to the mean downwelling irradiance. Observations with SZA larger than 75°are discarded in the comparison since at large zenith angle the ground solar irradiance measurements are degraded by cosine errors in clear sky (Albrecht & Cox, 1977). The MODIS atmospheric correction is not processed and the extrapolation from such a high SZA to a nadir illumination angle is also problematic for the MODIS BRDF model (Lucht, 1998).

Diurnal and seasonal cycle of shortwave radiation
Before examining surface albedo on climate timescales, it is useful to show the daily cycle of shortwave radiation at Summit, where the snow surface is more homogeneous than most of the other sites on Greenland (Fig. 3). On day 180 (June 28th) 2005, the SZA changes from 50°to 84°at Summit. Downwelling and upwelling shortwave radiation peaks at local noon. Snow albedo generally decreases with SZA, and reaches its minimal near local noon. This diurnal cycle with minimal surface albedo near noon is consistent with other reports including desert soil and grass land field observations (Jin et al., 2003b;Painter & Dozier, 2004;Wang et al., 2005;Liu et al., 2009) and with model predictions (Lucht, 1998). However, snow albedo does not always have the lowest value at local noon. As explained above in Fig. 1, atmospheric water vapor and ozone absorption shift the spectral intensity of surface insolation towards more reflective wavelengths, and many other factors, such as cloud, fresh snow addition, surface features, snow metamorphism and impurities, affect the SZA dependence of snow albedo. More importantly, the potential uncertainties of LI-COR 200SZ measured snow albedo may exceed the effect of SZA on snow albedo, especially on hourly and daily scales. Air temperature has similar shape as the cosine of the SZA with a two-hour phase shift to the right (not shown). The air temperature ranges from − 50°C to −5°C and is always below zero year-round at Summit.
The seasonal cycle (Fig. 4) shows that Summit is dark (SZA > 90°) from November through January. Summit receives attenuated illumination for about another three months (February and March 65°< SZA < 90°, and October 75°< SZA < 90°). During these 180 days, Summit receives less than 10% of its annual insolation. Summit has more cloud cover in the late summer and fall than the late winter and early summer.
Daily snow albedo at Summit exhibits large fluctuations in early spring and in late fall (Fig. 4). The seasonally symmetric portion of this variability is explained by the increasing importance of random noise in the decreasing (with increasing SZA) up-and down-welling fluxes whose ratio determines the albedo. This high ratio of noise to signal at high SZA may also explain the poor performance of the MODIS BRDF model. Hence the in situ albedo values for SZA > 65°must be used cautiously when evaluating satellite-retrieved albedos. From mid-April through September, snow albedo fluctuates little along with SZA   . 4. Seasonal cycle of daily mean SZA, cloud fraction, noontime albedo, relative insolation (fraction of the annual radiation) and its accumulated frequency at Summit (SMM). The right vertical scale is for the shortwave albedo and the shortwave downwelling insolation fraction. Daily radiation is the 24-hour average, and the daily mean is the 10-year average from 1997 to 2007 without 2001 on each day. Both albedo and frequency are finally computed using the daily mean radiation. The shortwave upwelling radiation frequency has similar pattern with the downwelling radiation frequency and is not shown on the figure. Mean cloud cover is the average of daily cloud-cover fraction from MODIS daily snow cover product (MOD10C1) from 2003 to 2007. The GC-Net in situ albedo data from the year 2001 were anomalously ∼ 10% higher than the 10-year mean and therefore are discarded from our analysis. and reaches its minimum in summer. As SZA increases by 21°from 49°( day 170) to 70°(day 258 or 86), the albedo increases from 0.82 to 0.86, similar to our simulations that the shortwave snow albedo increases from 0.76 to 0.8 when SZA increase from 55°to 80° (Fig. 1). Stroeve et al. (2005) also find that the GC-Net snow albedo values increase with SZA and latitude.

SZA dependence of and air temperature impacts on snow albedo
SZA is the most important environmental parameter in explaining seasonal variations of noontime surface albedo measured in the permanently covered, dry snow environment at Summit. The SZA cosine explains 43% of the springtime snow albedo variations from 1997 to 2007 (Fig. 5A). Note that the GC-Net in situ albedo data from the year 2001 were, for unknown reasons, anomalously~10% higher than the 10-year mean and therefore are discarded from our analysis. In contrast, daily air temperature does not have a significant relationship with daily snow albedo (Fig. 5B). This suggests that the surface air temperature-induced grain metamorphism (warming leads to fewer, larger grains) is insignificant in terms of its albedo effect (or is masked by other processes) at the daily scale on Summit (Flanner & Zender, 2006). The extremely low temperature at Summit (−50°C to −5°C) is less conducive to grain metamorphism than warmer temperatures (Taillandier et al, 2007). Table 2 further summarizes the correlations of snow albedo with cosine (SZA), air temperature and cloudiness during spring at Summit. Cosine of SZA has a significantly negative relationship with springtime snow albedo at Summit in eight of the ten years (that exclude 2001) from 1997 to 2007. In contrast, daily air temperature does not have significant correlations with snow albedo in seven years, and only has weakly significant positive correlation with snow albedo in three years. Cloudiness has only one year for significant correlation in five years. Theoretically, in springtime (days 86-172) both the increase of air temperature (which accelerates metamorphism) and the decrease of SZA contribute to the decrease of snow albedo (cf. Fig. 1). Processes that affect albedo but which we have not considered (due to lack of data) include meteorology (which determines initial snow grain size), surface processes (e.g., hoar formation, sastrugi), and transport of pollutants. In summary we find that the modest (4%) decline in surface albedo at Summit from spring to summer is best explained by the SZA effect, consistent with our model simulation (Fig. 1). The same behavior and explanation applies at the other stations.
In most years (1997, 1998, 1999, 2002 and 2004), negative seasonal-mean albedo anomalies were accompanied by positive seasonal-mean air temperature anomalies (Fig. 5C). In 2005 and 2007, positive seasonal-mean albedo anomalies were accompanied by negative seasonal-mean air temperature anomalies. This indicates that air temperature (far below 0°C) often does reduce snow albedo on longer timescales than daily at Summit. Black carbon and dust  aerosol could contribute to the negative albedo anomalies if sustained anomalously strong deposition occurred. While 1997 and 1998 were strong biomass burning years in Eurasia and North America ( Van der Werf et al., 2006;Randerson et al., 2006), we have no evidence that these aerosols made a difference at Summit (Flanner et al., 2007).

Comparison between GC-Net and MODIS albedo
Both GC-Net and MODIS-retrieved surface albedos are available for intercomparison at five permanently dry snow-covered stations in the five years from 2003 to 2007 (Fig. 6ABCEF). Only days with noontime SZA less than 70°are shown for reasons: (1) GC-Net in situ irradiance measurements are disturbed by noise at high SZA (Fig. 4); (2) the MODIS radiative transfer model and albedo retrieval algorithm do not work well at high SZA (Lucht, 1998;Schaaf et al., 2002;Wang et al., 2005). At South Dome (SDM), NGRIP(NGR) and Humboldt (HMG), MODIS albedos are always less than the GC-Net measurements by mean absolute differences of -0.04, − 0.12 and − 0.14, respectively when SZA is less than 60°(see Fig. 6AEF and Table 3). MODIS albedo most closely matches the GC-Net measurements at Saddle (SDL) with a mean absolute difference of −0.003 when SZA is less than 60°( Fig. 6B and  6C and Table 3). The discrepancy is relatively steady at other stations. This reminds us that the GC-Net instruments have 5% uncertainties for the downwelling and upwelling shortwave irradiance (Box & Steffen, 2000;LI-COR, 2005). The mean discrepancy increases from − 0.003 to − 0.14 for SZA < 60°with increasing station latitude due to the increased in situ albedo measurements and reduced MODIS albedo estimates. Similarly, the discrepancy dramatically increases at each station when SZA > 60°, mainly because MODIS albedo falls sharply as SZA increases above 55-60°( Fig. 6 and Tables 3 and 4). This is consistent with Stroeve et al. (2005) who find that MODIS snow albedo is reliable at low SZA, and that large errors (up to 0.2) may occur at high SZA. In contrast, we find systematic negative and larger bias of MODIS snow albedo compared to the GC-Net measurements at more northerly stations, instead of a systematic positive bias of 0.07 for the northerly stations in years from 2000 to 2003 for both main and backup retrievals (Stroeve et al., 2005), including Humboldt and NGRIP that show negative biases in 2003 through 2005 in this study.
Besides SZA and station latitude, the BRDF model retrieval quality also affects the MODIS albedo estimates. A retrieval quality flag (brdfq) of zero represents the highest quality and four the worst quality. The best quality retrievals (brdfq = 0) at Summit result in smaller (though usually still non-zero) albedo discrepancies (Fig. 6CD). The mean discrepancy is bigger at larger SZA (60-70°, Table 4) than smaller SZA (40-59°, Table 3) for all five stations. However, retrieval quality alone cannot explain the larger discrepancy because the mean discrepancy is still negative at Summit (from −0.03 to −0.07) and the other four stations when only considering the best retrievals (brdfq = 0). Thus SZA and the (related) retrieval quality both contribute to the larger discrepancy at large SZA.

Influence of SZA on MODIS albedo
GC-Net albedo measurements (Fig. 7A) and MODIS albedo estimates (Fig. 7B) at five stations in 2005 demonstrate the errant behavior of MODIS albedo at large SZAs. During the same period (each data point represents the mean clear-sky noon time albedo in the 16 days), contrary to the GC-Net measurements, MODIS albedos are systematically lower at higher latitude stations. The northernmost station, Humboldt, has the lowest MODIS albedo and yet has nearly the largest GC-Net albedo. When SZA < 55°, both MODIS and GC-Net albedos tend to increase with SZA as expected (cf. Fig. 1). However, MODIS albedo starts to dramatically decline when SZA exceeds 55-60°while GC-Net albedo still increases or remains relatively stable with SZA.
In order to demonstrate that the biases in MODIS albedo are systematic with SZA and extend beyond the 3 × 3 grid area surrounding the stations, we now examine the behavior of the mean MODIS albedo in ten-degree longitude (310-320E) by one-degree latitude swaths. Both MODIS white-sky broadband and narrowband (band4 = 545-565 nm) albedos increase with SZA for SZA < 55-60°i n 2005 (Fig. 8AB). Other years and bands, not shown, are similar. Of course, other factors besides SZA influence the snow albedo. For instance, as spring unfolds and the SZA decreases, air temperature, snow grain size and impurity content may all increase and compound the albedo decrease due to SZA (Flanner & Zender, 2006). As summer turns to fall and SZA increases again, fresh snow and more cloud (cf. Fig. 4) could also increase the snow albedo. No systematic SZA effect is discernible in either the broadband or narrowband MODIS albedos during summer (c. days 153-185) when SZA < 50°at all stations. This agrees with Petzold's (1977) empirical rule-of-thumb that snow albedo is virtually independent of SZA at SZA < 50°.
Broad and narrow band WSAs dramatically decrease with increasing SZA in late winter before day 113 and in fall after day 233 ( Fig. 8A and B). The MODIS-retrieved surface albedo in the northernmost zone (79°N) is always the darkest, while the southernmost zone (67°N) is always the brightest or second-brightest. The decline in broadband and narrow band WSAs is nearly coincident for each latitude zone. Hence the MODIS-retrieved albedos for large SZAs are self-consistent but are inconsistent with both theoretical estimates (Fig. 1) and in situ measurements (Figs. 3-4). In addition, both broad and narrow band WSAs are for the whole retrievals, including the "backup" and "main" algorithm retrievals. Similarly, systematic patterns of declining albedo at high SZA are also observed even for the "main" algorithm retrievals (not shown). Finally, Blacksky albedos behave similarly to, and have slightly smaller magnitude of decrease with SZA increase than their white-sky counterparts and therefore are not shown here.
It is instructive to examine the geographic distribution of the 16-day white-sky albedo (Fig. 9). The white areas are missing data. Coastal areas are snow free and have albedo < 0.5 (0.2-0.49). On day 65, when SZA > 75°, albedo is strongly correlated with latitude (and SZA). For latitudes north of 75°N, snow albedo is usually < 0.6. Such low albedos are not realistic; they must be discarded or adjusted for correct interpretation. The major reason for such distorted albedos is the MODIS BRDF model which has difficulty at high SZA (> 70°) because of the high noise-to-signal ratio. Further, the MODIS atmospheric correction is not applied for SZA > 75° ( Lucht, 1998). Similar albedo declines with latitude/SZA are evident in mid-spring (day 105) and early-fall (day 233). In summer (day 169 = June 18) when the SZA < 60°, snow albedo is much more uniform over Greenland (Fig. 9). The higher albedo at the southern crest may be due to fresher snow or surface hoar. The pattern of highest albedos coincident with smallest SZA is relatively insensitive to season and to absolute SZA.

Snow albedo dependence on SZA
Both theoretical model simulations and in situ measurements show that land surface albedo (snow, desert sand, and vegetated area) has a modest dependence on SZA, i.e., their albedo increases with SZA Dickinson, 1983;Jin et al., 2003b;Stroeve et al., 2005;Wang et al., 2005;Liu et al., 2009). For clear-sky conditions, the diurnal variation of snow (and vegetation) albedo generally shows a wide "U" shape since photons penetrate deeper into the surface and are more likely to be absorbed at smaller SZA, resulting in lower albedo (Kimes et al., 1987;Liu et al., 2009). For On plot D, the right vertical scale is for the mean BRDF retrieval quality code at 9 cells, 0 represents the best retrieval quality with >75% of full inversion, 2 represents 51-75% of full inversion, 3 represents 25-50% of full inversion, and 4 is the worst retrieval quality from the backup algorithm. In situ data in 2006 and 2007 are not available for both NGRIP and Humboldt stations. larger SZAs, photons interact with the snow grains over a longer path allowing more multiple scattering and less penetration through and absorption by snow Wang et al., 2005;Lucht, 1998). The theoretical dependence of snow albedo on SZA is <0.06 for a SZA increase from 0°to 90° (Fig. 1), and the GC-Net measurements show that the snow albedo change due to SZA and other factors is <0.04 for the observed SZA increase from 50°to 70°at Summit (Fig. 4).
The in situ GC-Net measurements show that the diurnal snow albedo is not always minimal at local noon, and that the annual minimum is not always at the summer solstice. There are at least three reasons why snow albedo does not always increase with SZA (Warren, 1982). Fig. 1 demonstrates the first two reasons. First, atmospheric water vapor and ozone absorption deplete near infrared (NIR) insolation which snow would otherwise likely absorb, and shift the spectral distribution toward visible wavelengths where snow is more reflective. Hence increased water vapor and ozone result in higher snow albedo. Second, clouds shift the spectral distribution of solar radiation by ice and liquid water like water vapor and change the effective incident angles by diffusing the incident solar radiation (Liu et al., 2009). Thus the broadband surface albedo under clear sky is lower than under cloudy sky (see Fig. 1). Third, and most important at explaining the spatial heterogeneity of snow albedo, are the surface features of snow (grain size distribution, black carbon, crusting, and sastrugi) that depend on local meteorology, roughness and topography. The model simulated albedo difference between 50 µm and 100 µm of effective snow grain radius is 0.02 (Flanner & Zender, 2006). The albedo reduction by black carbon (BC) in snow albedo is less than 0.01 for measured BC of 14.6 ng/g at Summit (Slater et al., 2002;Flanner et al., 2007). The sastrugi effect is up to 0.05 at the daily scale and 0.01 at the monthly scale in Antarctica (Pirazzini, 2004). The largest snow metamorphism, such as snow melting and refreezing, is excluded from our study sites although its effect on the snow albedo is up to 0.2 in the diurnal cycle around the freezing point (Pirazzini et al., 2006). Snow albedo variations caused by these properties likely exceed the relatively modest SZA effects (cf. Fig. 1) on small spatial and short temporal scales as shown by Hall et al. (2009b). Last, the potential uncertainties of LI-COR 200SZ measured snow albedo may also exceed the impacts of SZA on snow albedo, especially on the daily and hourly scales.

MODIS snow albedo bias at large SZA
The steep decline of MODIS snow albedo at large SZA (55°-75°) is physically unrealistic and is most likely a retrieval artifact (Figs. 6-9). MODIS-retrieved albedos at large SZA are self-consistent but are inconsistent with both theoretical estimates and in situ measurements. First, the in situ measurements of both the diurnal cycle ( Fig. 3) and seasonal cycle (Fig. 4) show that surface albedo remains constant or increases slightly, on average, as SZA increases. Second, theoretical predictions ( Fig. 1) of both the SZA effect and the spectral shift effect predict a slight increase of albedo with SZA. Third, the MODIS albedos are anomalously low, relative to in situ albedos, in both early spring and late fall (Figs. 6-9) for both "main" and "backup" algorithm retrievals. If melting surfaces contributed to the lowreported albedos (a possibility all but ruled out by our site selection and temperature criteria), this would likely occur in summer, not in late winter or late fall when air temperatures are far below 0°in Greenland. Stroeve et al. (2005Stroeve et al. ( , 2006 show that high-precision solar radiation measurements at Summit are systematically lower (and with a smaller RMSE uncertainty of 3.5%) than GC-Net albedo. However, such systematic biases only help explain the absolute discrepancy of MODIS albedo with GC-Net, and do not exclude the systematic and artificial decline of MODIS snow albedo at large SZA.
What causes the MODIS snow albedo decline as SZA > 55-60°? The SZA and (related) retrieval quality both appear to contribute to the decline at large SZA. A high SZA poses serious challenges not only to MODIS retrievals but also to in situ albedo measurements. As SZA increases, the surface solar radiation decreases, and the signal-tonoise ratio decreases, resulting in higher albedo change frequency and larger magnitude fluctuations. When SZA > 65°, the in situ snow albedos fluctuate at higher frequency and with larger magnitude compared to SZA < 65° (Fig. 4). Similarly, the MODIS sensors cannot effectively detect the reflected solar radiation at large SZA because of the reduced downwelling and upwelling solar radiation. In addition, the large MODIS view zenith angle ranges (0-55°), may contribute errors when converting the MODIS sensor retrieval from a large view zenith angle to a nadir view zenith angle, especially at a high SZA. In the legend, the first three letters are the station ID, and the last two numbers are their latitude (degrees North). SZA_HMG and SZA_SDM are the solar zenith angles at the northern-most station (Humboldt) and the southern-most station (South Dome), respectively, and other stations' SZAs are between them, proportional to their latitude. Minus sign means MCD43C3 albedos are lower than GC-Net albedos. Note: n/a means that the ground AWS data were not available.
Usually MODIS albedo starts to decline as SZA > 55°, though, in some years and stations (without any clear pattern), the decline may not commence until SZA exceeds 60°or even 65°. The SZA distortion of surface albedo for high SZA also exists at other regions in the world, such as Northern Eurasia and the Antarctic. These regions are not shown for brevity. However, a recognition and correction of the biases in the current MODIS snow albedo is necessary and important in order to make full use of the MODIS albedo product. Using the MODIS albedo to evaluate simulated albedo in snow covered regions with SZA > 55°could lead to erroneous conclusions. Oleson et al. (2003) conclude that the Community Land Model (CLM) overestimates visible (VIS) and near infrared albedo in North America, Eurasia north of about 55°, and Greenland. Zhou et al. (2003) conclude that the CLM overestimates winter MODIS VIS albedo by about 0.2 over Greenland and northern Canada. It is likely that in both cases the MODIS BRDF model underestimates the real snow albedo. Identifying the bias in MODIS snow albedo at large SZA highlights the current limitations of the MODIS snow albedo dataset, and should help reduce inappropriate conclusions drawn from data at large SZA. A logical next step, beyond the scope of the present study, would be to parameterize a large SZA correction to MODIS albedo. Such a correction would allow us to better harness the valuable MODIS albedo dataset for cryospheric monitoring and for climate model evaluation and parameterization.

Summary
In situ measurements of snow albedo at five stations along a northsouth transect in the perennially dry-snow covered interior of Green-land follow the theoretically expected dependence of snow albedo with SZA. GC-Net shortwave downwelling and upwelling solar radiation measurements from 1997 through 2007 exhibit the expected trend of modest surface brightening with increasing SZA on both diurnal and seasonal timescales. SZA explains up to 50% of seasonal albedo variability. The two other environmental factors explored, temperature and cloudiness, rarely play a statistically significant role in seasonal albedo variability at the five stations surveyed.
Compared to the 10-year record of these GC-Net measurements, a five-year record of MODIS satellite-retrieved snow albedo shows a systematic negative bias that increases with SZA for SZA larger than about 55°. The consistency and accuracy of the MODIS radiative transfer methods and BRDF models face serious challenges at high SZA, especially in polar regions. When SZA is smaller than 55°, MODIS albedos successfully capture the snow albedo dependence on SZA and show a relatively good agreement with GC-Net measurements. The discrepancy of MODIS with in situ albedo and with theory is caused mainly by two related factors, SZA and retrieval quality, that depend on both location and season. Although the in situ instruments suffer from significant uncertainties themselves, in aggregate they provide clear and compelling evidence for an artificial behavior in MODIS snow albedos at large SZA. Consequently, the current MODIS snow albedo estimates for SZA> 55°should be regarded skeptically. They are self-consistent and likely capture the morphology of spatial features correctly, but the accuracy of the snow albedos deteriorates for SZA > 55°and becomes physically unrealistic for SZA> 65°. Seasonal-to-interannual trends in surface brightness in Greenland, and in polar (i.e., large SZA) regions in general, and model simulations of these trends, should be evaluated in light of these limitations. Cover Fraction (SCF) is 100% except for latitude 67°in summer (SCF minimum is 98.8%). Albedo values before the 65th day and after the 281st day are not illustrated because of the high SZA (> 75°) and correspondingly limited solar incident radiation, resulting in unreliable albedo retrievals.