Constraining MODIS snow albedo at large solar zenith angles: Implications for surface energy budget in Greenland

An understanding of the surface albedo of high latitudes is crucial for climate change studies. MODIS albedo retrievals flagged as high quality compare well with in situ Greenland Climate Network (GC-Net) measurements but cover too little area to fully characterize Greenland's albedo in non-summer months. In contrast, poor quality MODIS retrievals provide adequate spatio-temporal coverage, but are not recommended for use at large solar zenith angles (SZAs) where they have a systematic low bias. We introduced an empirical adjustment to the poor quality data based on high quality reference albedos and constrained by GC-Net data and theory, and used the adjusted data to improve estimates and fill-in gaps of the year-round, Greenland-wide, albedo and surface energy budget. For observations made with SZAs between 55° and 75°, the mean differences (MODIS minus GC-Net) between our adjusted MODIS albedo and GC-Net measurements are -0.02 and -0.03 at Saddle and Summit, respectively, compared to -0.05 and -0.08 between the unadjusted MODIS albedo and GC-Net measurements. The adjusted MODIS snow albedos are usually between 0.75 and 0.87 over dry snow when SZA is larger than 55°, and they reduce unrealistic seasonal and meridional trends associated with MODIS retrievals at large SZA, defined as SZA > 55 ° and 70 ° , respectively, for low- and high-quality retrievals. The impact of the adjusted albedo on the surface energy budget, relative to the unadjusted albedo from all MODIS data, is least (-0.7±0.1W/m 2 ) in June, and greatest (-6.2±0.9 W/m 2 ) in September for the black-sky albedo (BSA). The mean annual absorbed solar radiation (ASR) reduction by the adjusted MODIS albedo in Greenland from 2003 to 2005 is 3.1±0.2 and 4.3±0.2 W/m 2 for BSA and white-sky albedo (WSA), respectively, about 8.0±0.5% and 10.8±0.4% of ASR based on the raw BSA and WSA. The ASR reduction by the adjusted blue-sky (actual) albedo is between 2.9 and 4.5 W/m 2 , enough to annually melt 27.1 to 41.7 cm snow water equivalent (SWE), or to sublimate 3.2 to 4.9 cm SWE. The ASR difference between the adjusted MODIS BSA and CERES albedo in March from 2003-2005 is only -0.1± 0.9 W/m 2 , much less than the difference (4.9±1.4 W/m 2 ) between the unadjusted MODIS BSA and CERES. The albedo adjustments exceed the likely direct anthropogenic radiative forcing experienced by Greenland due to greenhouse gases or aerosols. The proposed adjustment preserves most of the zonal and meridional structure of raw MODIS albedo, and extends its usefulness as a cryospheric climate record in times and regions of Greenland with large SZA.


Introduction
Snow covers most of Greenland year-round, and thus plays a pivotal role in determining the surface energy balance of Greenland which, by virtue of its area, ice-volume, and location near regions of North Atlantic deep-water formation, plays an important role in the climate system (Steffen and Box, 2001). The minimum (summer solstice) Solar Zenith Angles (SZAs) in southernmost and northernmost Greenland are 37° and 60°, respectively. Such high SZAs pose serious challenges not only to in situ solar radiation measurements (Augustine et al., 2000), but also to the consistency and accuracy of the MODerate resolution Imaging Spectroradiometer (MODIS) surface albedo retrieval algorithms applied to polar regions (Lucht, 1998;Liu et al., 2009). Due to its length (about 10 years) and relatively high spatio-temporal resolution (500 m pixels every eight days), the MODIS-retrieved surface reflectance has become a valuable climate data record for monitoring and evaluating the significance of snow albedo, and thus surface energy, changes in remote regions.
Theoretically snow surface albedo increases with SZA because the increased path over which obliquely-incident photons may interact with the snow grains allows more multiple scattering and less penetration through or absorption by snow. This results in a larger fraction of oblique solar radiation reflected by snow (Wiscombe and Warren, 1980;Warren, 1982;Lucht, 1998;Wang et al., 2005).
Several field studies show surface albedo increases with SZA over snow, desert, and vegetated surfaces (Warren and Wiscombe, 1980;Jin et al., 2003b;Wang et al., 2005;Liu et al., 2009). Although the MODIS albedo increases with the SZA increase at low SZA (Lucht, 1998;Jin et al, 2003b), Liu et al. snow albedo larger than 0.8 as SZA > 60°. In contrast, the MODIS-retrieved snow albedos in areas of Greenland known to be covered with dry snow dip as low as 0.6 for SZA > 70° for the magnitude inversions (quality flag Q>0) while theory and GC-Net measurements show that the snow albedo remains unchanged or increases slightly as SZA increases (Wang and Zender, 2010). It must be noted that Wang and Zender (2010) did not discriminate the retrievals by quality flag. More careful analysis shows that the biases they identified are attributable to poor quality retrievals for 55° < SZA < 70° (where the high quality retrievals perform well), and by all retrievals for SZA > 70°, which is beyond the recommended range of use of the albedo product (Stroeve et al., 2005). Liu et al. (2009) attribute the decline of the MODIS albedo with large SZA to the extrapolation algorithm and/or to the Bidirectional Reflectance Distribution Function (BRDF) model over snow surfaces. Wang and Zender (2010) concluded that the accuracy of MODIS albedos deteriorates for Q>0 data when SZA > 55° and for Q=0 data when SZA > 70° and that these albedos can be physically unrealistic for SZA > 70°.
Many studies intercompare GCM modeled albedo results with the combined (comprising both highand low-quality retrievals) MODIS dataset (Zhou et al., 2003;Oleson et al., 2003;Roesch 2006). These earlier uses of the combined MODIS dataset were more exploratory, and their comparisons over Greenland, Antarctica, and other high latitude snow-covered regions should be re-considered in light of subsequent studies that document low biases at large SZA (Stroeve et al., 2005;Liu et al., 2009;Wang and Zender, 2010).
Despite the problems noted above with retrieved snow albedos for SZA > 55-70°, an extensive literature search finds good agreements between the MODIS-estimated surface albedo and in situ observations. At most Surface Radiation Budget Network (SURFRAD) sites, the overall absolute accuracy of MODIS albedo is within 0.05 and shows an increasing negative bias and increasing root mean square error (RMSE) compared to ground observations as SZA increases beyond 70° (Liu et al., 2009). At the Gaize Automatic Weather Station on the western Tibetan Plateau with semi-desert or desert soil, the MODIS albedo does not show a distinctive bias with the ground-measured albedo (Wang et al., 2004). At the SURFRAD sites and Cloud and Radiation Testbed-Southern Great Plains (CART/SGP) sites, the MODIS surface albedo generally meets an absolute accuracy requirement of 0.02 with an RMSE less than 0.018 during April-September 2001(Jin at al., 2003b. In snow-free periods, MODIS albedo shows good agreements with worldwide Baseline Surface Radiation Network (BSRN) measurements (Rosech et al., 2004). At perennially snow-covered sites in Greenland, MODIS retrieves snow albedo with an RMSE of 0.04 for high-quality retrievals and 0.07 for all retrievals with SZA<70° in 2000-2003, relative to GC-Net measurements, which have an RMSE of 0.035 relative to the BSRN measurements at Summit (Stroeve et al., 2005). Moreover, MODIS snow albedos at large SZAs are self-consistent and are likely to capture the underlying spatial morphology, especially its zonal features.
Yet the best quality (Q=0) MODIS data cover less than 40% of Greenland on each day, and annually only 33% of all retrievals on days 41-297, 2005, when the majority (56%) of retrievals were quality level Q=2 data (11% are for data whose Q=1, 3, and 4). The best quality (Q=0) data congregate near smaller SZAs, i.e., towards southern Greenland and in less-cloudy months (e.g., May and August), when the best quality data is near 50% of all retrievals. The areal fractions of Greenland with the best quality data coverage in March, April, and May, are 13%, 30%, and 36%, respectively. Using data from all quality levels (Q=0-4) increases these monthly areal coverages to 80%, 93%, and 98% respectively. For Greenland north of 72° (i.e., north of Summit, which includes about 48% of Greenland's total area), the best quality data coverages in March, April, and May drop to 3%, 24%, and 35%, respectively, and utilizing lower quality data increases these to 75%, 91%, and 97%, respectively.
Since restricting our analysis to only those cells that are of best quality would eliminate too much data, we use all of the data, as in Oleson et al. (2003), to obtain a full spatial and temporal coverage in Greenland.
Our goal is to ensure that the community of scientists interested in surface processes in the cryosphere, and Greenland in particular, becomes aware of the importance of retrieval quality in assessing MODIS snow albedo and of the existence of the MODIS snow albedo bias at large SZAs (defined as SZA > 55° and 70° for low-and high-quality retrievals, respectively) and gains a good sense of its magnitude in Greenland (where we have data to evaluate it) and its implications there for surface processes such as sublimation or snow melt. The strategy we pursue to increase the usefulness of MODIS snow albedo to this community of researchers is to develop empirical adjustments that retain the accurate MODIS retrievals up to SZA = 55°, rely on high quality retrievals where possible for SZA between 55° and 70°, and mitigate the large SZA bias elsewhere. We first describe the patterns of the low-bias in MODIS snow albedo at large SZA in Greenland. Then we develop an empirical model to adjust the archived MODIS albedo at large SZA in order to improve its usefulness as a cryospheric climate record. Finally, we evaluate the corrected MODIS albedo using in situ measurements and Clouds and the Earth's Radiant Energy System (CERES) surface albedo, and demonstrate its implications for the surface energy budget in Greenland.

Study Area and GC-Net Albedo
The Greenland ice sheet is an ideal target to study snow and ice albedo from satellites and has been used as such in several studies ( Knap and Oerlemans, 1996;Stroeve et al., 1997Stroeve et al., , 2001Greuell and Oerlemans, 2005). First, Greenland has about 20 solar radiation monitoring stations that sample different snow zones. Second, the vast majority of Greenland is perennially snowcovered without the disturbance of vegetation. Third, outside the (increasing) melt-zone (Tedesco, between the in situ footprint and the satellite areal retrieval. Fourth, Greenland is far downwind from extensive human activities and dust sources. Concentrations of black carbon impurities, the most important absorbing impurity in present day Greenland, are less than 15 ppb and are thought to reduce broadband albedo by less than 1% (Flanner et al., 2007, Table 2). Fifth, Greenland's year-round high reflectivity in the midst of the dark North Atlantic Ocean plays an important role in the polar climate and in fresh water storage. Greenland's surface climatology and geographic features are documented by Steffen and Box (2001).
There are 21 GC-Net Automatic Weather Stations (AWS) in Greenland (Figure 1), which provide downwelling and upwelling shortwave irradiance measurement. Five stations (Humboldt GI-HMG, NGRIP-NGR, Summit-SMM, Saddle-SDL and South Dome-SDM) are selected in this study (Table 1). These five stations are located along the crest of the ice sheet and span Greenland from south to north at elevations close to or above 2000 m. Restricting our study to dry snow regions ensures that snow-melt does not contribute to the discrepancy between in situ and satellite-retrieved albedos.
The GC-Net shortwave solar downwelling and upwelling radiation are measured using a pair of horizontally leveled LI-COR 200SZ pyranometers in a narrow spectral range (0.4-1.1 µm) sampled at a 15-s interval and averaged over an hour. The LI-COR measurements have 5% uncertainties. The downwelling shortwave solar radiation value measured by a LI-COR 200SZ pyranometer is factorycalibrated to equal the spectral response from 0.28-2.8 µm measured by a more accurate Eppley Precision Spectral Pyranometer (PSP) (LI-COR, 2005;Stroeve et al., 2005). The PSP measures over 98% of the downwelling shortwave solar radiation.
The LI-COR pyranometers are calibrated against Eppley PSP measurements based on the spectral distribution of the downwelling solar radiation under clear sky. Thus systematic positive biases may exist in the GC-Net measured upwelling irradiances since the snow surface reflects over 90% of the visible solar radiation and depletes most radiation beyond 1.1 µm. Stroeve et al. (2005) compared the solar fluxes measured by LI-COR and Eppley pyranometer at Swiss Camp in Greenland, and found the upward and downward shortwave irradiance errors measured by LI-COR pyranometers did not exceed 2.7%. To compensate for the bias in the reflected LI-COR 200SZ pyranometer measured irradiances, they corrected the GC-Net albedo with a site-specific albedo offset 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 because the 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 correcting the GC-Net snow albedo, none of the downward or the 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. To avoid unexpected instrumental errors, we apply an outlier check for the hourly data such that albedo dropping more than 10% against the mean of the two neighboring hours when SZA < 80° is replaced by the mean or discarded as non-credible data.
In order to compare with the local noon MODIS albedo, the daily in situ albedo is derived from the average of shortwave downwelling and upwelling irradiance within a three-hour period centered on local noon (11:00, 12:00 and 13:00). The 16-day albedo is derived from the daily average of downwelling and upwelling solar radiation in clear skies for each 16-day period of MODIS albedo.
Clear sky is defined when the cloud cover fraction at the MODIS gridpoint coincident with the GC-Net site is less than 10% according to the daily snow cover product (MOD10C1), which contains the cloud cover fraction (Wang et al., 2008).

MODIS Albedo
The MODIS instruments aboard both NASA's Terra and Aqua satellites acquire daily images and provide global land surface albedos that include the daily unvalidated beta-test albedo product and the 16-day validated standard product (Schaaf et al., 2002;Stroeve et al., 2006). The MODIS standard 16-day surface albedo products are produced every 8 days within 16 days of acquisition. There are multi-angle observations in the 16-days. If the majority of observations are recorded as snow-covered, then the algorithm uses only snow-covered observations for the parameter retrievals; otherwise, the algorithm conservatively uses snow-free observations for parameter retrievals (Schaaf et al., 2002;Salomon et al., 2006). Meanwhile, if there are at least seven cloud-free observations during the 16 days, a full model inversion or the "main" algorithm is attempted; otherwise, a magnitude inversion or the "backup" algorithm is performed (Schaaf et al., 2002). MODIS albedo products provide directional hemispheric reflectance (black-sky albedo, BSA) and bihemispheric diffuse reflectance (white-sky albedo). Each BSA and WSA include seven narrow spectral bands (MODIS band 1-7) and three broadbands (0.3-0.7, 0.7-5.0, 0.3-5.0 µm). The broadband albedos are converted from the spectral reflectance via a narrow-to-broadband conversion factor (Liang et al., 1999;Stroeve et al., 2005). The shortwave broadband (0.3-5.0 µm) albedo (hereinafter referred to as MODIS albedo) has the best spectral match to most broadband instruments and is of the greatest interest to climate studies because it determines the net surface solar radiation flux. The actual albedo (also called blue-sky albedo) is a solar flux-weighted average of the intrinsic black-sky albedo and white-sky albedo, where the proportion of direct and diffuse solar radiation depends on the atmospheric conditions and SZA, usually reported at local noon for MODIS (Lewis and Barnsley, 1994;Schaaf et al., 2002). The MODIS Terra and Aqua combined albedo products have higher quality retrievals than Terra-only albedo products (Salomon et al., 2006).
The MODIS Terra and Aqua combined albedo has been available since Aqua MODIS operational retrievals began in June 2002. MODIS albedo products come in two grids and are distributed by the NASA Land Process Distributed Active Archive Center (LP DAAC) at https://lpdaac.usgs.gov/lpdaac/get_data/wist. One is a sinusoidal projection with 500 m and 1 km spatial resolution, and the other is the climate modeling grid (CMG) at 0.05 degree (about 5.6 km near the equator), which is aggregated from the 500 m product. Following Wang et al. (2005), this study uses the combined 0.05 degree MODIS albedo product (MCD43C3) from 2003 to 2007. MCD43C3 contains WSA and BSA for 7 spectral bands and three broadbands, snow cover fraction (SCF), local noon solar zenith angle (SZA), and the BRDF quality code. The quality flag has five values for the aggregated 0.05° product in MCD43C3. Q=0 represents the best quality and indicates that 75% or more of the 500 m pixels that compose the aggregated point were retrieved with the full inversion algorithm; Q=1 points are good quality, though less than Q=0 quality, and are also composed of 75% or more full inversions; Q= 2 means mixed quality, comprising fewer than 75% full inversions and fewer than 25% fill values; Q=3 suggests all magnitude inversions with fewer than 50% fill values; while aggregated points with Q=4 include more than 50% fill values. For near-homogeneous bright snow surfaces, both BSA and WSA are similar, and each closely represents the blue-sky albedo when SZA is less than 55° (Stroeve et al., 2005;Wang and Zender, 2010). This study reaches similar conclusions from independent analyses of BSA and WSA data, and occasionally, for brevity, we illustrate our results for BSA only.  , 2003b), any conclusion about the accuracy of the MODIS albedo must consider this and the uncertainties of the in situ shortwave downwelling and upwelling radiation measurements.

CERES albedo
Like MODIS, the Clouds and the Earth's Radiant Energy System (CERES) instrument is aboard both Terra and Aqua satellites. CERES measures broadband shortwave solar radiances at the Top Of Atmosphere (TOA), while MODIS measures narrow-band radiances in various shortwave bands (Jin et al., 2008). CERES generates global 1° gridded Monthly TOA/Surface Averages (SRBAVG) datasets. This study uses the CERES Terra FM2 Edition2D SRBAVG, which infuses observations from Geostationary Narrowband Radiances (GEO) and MODIS cloud data products and represents the most robust CERES TOA/surface monthly mean flux product (Wielicki, et al., 1996).

MODIS Albedo Bias
Before introducing the proposed correction to MODIS albedo at high SZA, it is instructive to summarize how albedo changes with large SZA in both temporal and spatial dimensions. At a given location, like Summit, the SZA is smallest at sidereal noon in the diurnal cycle, and in summer for the seasonal cycle. At a given day or time, SZA is smaller at a lower latitude. Wang and Zender (2010) document the theoretical and observed behavior of snow albedo (both GC-Net and MODIS) in dry snow regions of Greenland on the diurnal, seasonal, and inter-annual timescales. They show that temporally, over dry snow-covered regions of Greenland, the MODIS snow albedo (Q>0) decreases steeply as SZA exceeds 55-60°, and can reach values below 0.6 when SZA >75°. Our proposed albedo adjustment takes advantage of our new finding ( Figure 2) that these albedo decreases correlate positively and significantly with the cosine (SZA) in all years examined for both spring and fall at Summit. The SZA cosine explains 24-89% of the variance in BSA decline for large SZA at Summit. At vegetated (grassland with few trees) SURFRAD sites, the negative bias of MODIS albedo also increases as SZA increases beyond 70° (Liu et al., 2009). Other stations and years behave similarly to Summit for both BSA and WSA (not shown). Regressions and variance explained (R 2 values) are computed separately in spring and fall at each station and in each year in order to preserve any difference due to different seasonal properties of snow. In each season, there are only seven or eight data points available for SZA from 55° to 75° (recall that MODIS albedos are provided every eight days). Because of the limited data in each season, each datum has a strong influence on the final R 2 value. As a result, the regression coefficients vary greatly. Linear fits are better (higher R 2 ) in the fall than in the spring.  Figure 4c). Regionally, the low-bias in MODIS relative to CERES albedo increases up to 0.2 at latitudes > 75°N (Figure 4c). The lower surface albedo implies that the surface Absorbed Solar Radiation (ASR) is, on average, 5.18 W/m 2 greater for MODIS than CERES (Figure 4d). This ASR difference is about 23% of the ASR implied by the raw MODIS albedo, and 6% of the surface downwelling solar insolation. In spite of its potential uncertainty (overestimate at latitudes > 75°N) and coarser spatial resolution, CERES snow albedo increases with latitude and SZA (Figure 4a) as predicted by theory and as seen in in situ measurements (Wang and Zender, 2010). Table 2 demonstrate that MODIS snow albedo (constructed from all quality level retrievals) has a systematic negative bias at large SZA (>55-60°). This bias is closely related to SZA, and in particular for the magnitude inversions with Q>0 at large SZA. The Greenlandwide MODIS albedo is credible during June and July and with the best full inversions (Q=0) which appear unbiased for SZA < 70°. However, Q=0 retrievals comprise only 35% of the retrievals in June and July in Greenland, and congregate towards smaller SZAs, i.e., southern Greenland. The sparsity of high quality retrievals at large SZA means that all retrievals must be used to realize the full spatiotemporal coverage afforded by MODIS.

Proposed Albedo Correction
We propose a short-term empirical method to mitigate the MODIS albedo bias and to quantify its implications on the surface energy budget, from which the climate/cryosphere community would benefit. The theoretical dependence of spectral snow albedo on SZA, snow grain size and three singlescattering properties is well-known and has been expressed analytically in the delta-Eddington approximation (Wiscombe and Warren, 1980). The required inputs for that expression are unavailable in practice, so we start from an empirical surface albedo representation that depends on SZA and one free parameter (Dickinson,1983;Briegleb et al., 1986;Wang et al., 2005).
where A is the black-sky albedo, ө is SZA, ө 0 is the reference SZA (usually near 60°), and A(ө 0 ) is the albedo at the reference SZA which depends on date and location, and C is a surface-dependent empirical parameter of order 0.1.
Though equation (1) has been widely used in modeling land surface solar fluxes (Pinker and Laszlo et al., 1992) and land-atmosphere coupled models ( Briegleb et al., 1986;Kiehl et al. 1998, Hou et al., 2002, it cannot directly incorporate the MODIS retrieved albedo for SZA > ө 0 and we would like to retain the spatial morphology of MODIS retrieved albedos for large SZA. We consider five factors in modifying Equation (1) to fit our needs. First, no adjustment should be made where the raw MODIS snow albedo is presumed to be robust (Q=0 & SZA<70°, and Q>0 & SZA < 55°) , and these robust data should be used, where and when available, to adjust the non-robust data. Second, where adjustments are made, they should be be consistent with in situ observed dry snow albedo at large SZA. Third, the adjustment should be consistent with the SZA contribution to snow albedo at large SZA. Fourth, the adjustment should retain, when possible, spatial morphology of the raw MODIS albedo. Fifth, adjustments should not be made when non-snow or wet snow surfaces are indicated. To address these issues, we propose Equation (2) to adjust the raw MODIS snow albedo at large SZA: Here A(ө) is the raw MODIS retrieved albedo at ө SZA, and  A 0 is the reference albedo, taken to be the zonal mean albedo at the reference SZA ( 55°) for 100% snow-covered gridpoints. Thus,  A 0 varies with date and location. Since we set the reference SZA to 55° in Equation (2), we also changed the coefficient of C in Equation (1) to 1.74 to guarantee that the adjusted albedo ( A' ө ) increases with SZA. Correspondingly, C is changed to 0.15 for pure snow cover. This value constrains the total snow albedo increase caused by a 25° SZA increase from 55° to 80° to 0.04, in agreement with GC-Net observations and model simulations (Wang and Zender, 2010).
In Equation (2), the reference albedo,  A 0 , is the most critical parameter. It is computed in one of two ways. When the SZA at latitude 63°N (SZA 63N ) is less than 55° (e.g., from day 105 to day 241 in Greenland), then  A 0 is the zonal mean albedo at SZA=55°. When SZA 63N is larger than 55° (e.g., before day 105 and after day 241), then  A 0 is taken as the zonal mean albedo at latitude 63°N, and if this value is less than 0.8 then it is further adjusted by Equation (3) below. We choose the reference SZA = 55° because MODIS snow albedo is largest, and is still robust even for the magnitude inversions, at this SZA (Liu et al., 2009;Wang and Zender, 2010). This choice is also consistent with Petzold's (1977) empirical rule-of-thumb that snow albedo is virtually independent of SZA for SZA <50°.  A 0 is the average albedo of all 100% snow-covered gridpoints of the best (Q=0) retrievals whose noontime SZA = 55° (or, if SZA 63N > 55°, whose latitude is 63°N), and whose albedo is greater than 0.75 on the day of the adjustment. The albedo threshold of 0.75, which is derived from the MODIS best retrieval snow albedo in June and July, for contributing to  A 0 sets a minimum adjusted snow albedo when SZA = 55°. We choose latitude 63°N as the southernmost point to calculate the reference albedo when SZA 63N is larger than 55° because it is near South Dome (63.15°N, 44.81°W), the southernmost GC-Net station that is perennially snow-covered. However, the MODIS reference albedo at 63°N is often relatively too dark (< 0.8) because the large SZA (>65°) there produces an artificially low raw MODIS albedo that itself must be adjusted. Thus we adjust any  A  N 63  <0.8 before use in Equation (2) by applying Equation (3) which uses a ramp function to remap the reference albedo so that  A 0 > 0.8 when SZA 63N > 55°.
( 3) Equation (3) is applied only if SZA 63N > 55° (e.g., before day 105 and after day 241) and produces  A 0 consistent with GC-Net measurements at South Dome. In practice, Equation (3b&c) is only applied when SZA 63N > 65° because mean albedo for best retrievals at 63°N is larger than 0.8 when SZA 63N is less than 65°.  A 0 generally increases with SZA since  A 63N from MODIS decreases with SZA as SZA > 65°.
Equation (2) adjusts the original MODIS retrieved albedo A( ) ө by an offset scaled to the albedo difference between  A 0 and A( ) ө . In practice we apply equation (2) only at gridpoints that meet between 0.5 and 0.8, and (5) SCF=100% at the same grid but on day 161. The last three conditions help to exclude wet or melted snow by only allowing corrections at perennially snow-covered gridpoints.
The original MODIS snow albedo is not adjusted when any one of the above five conditions is not satisfied. This occurs mainly in the warm months in June and July. Condition (3) preserves the best retrieval (Q=0) MODIS snow albedo for SZA < 70°, and adjusts the albedo, using raw data of all qualities, when SZA ≥ 70°. Finally, we fill-in missing albedos with the mean adjusted albedo at the same latitude.
The maximum adjusted snow albedo produced by this algorithm is 0.87. A fresh dry snow albedo of 0.87 under clear skies is consistent with Stroeve et al. (2005). The lower bound on the adjusted snow albedo allowed by this algorithm is 0.75, which is also the lower bound of the reference albedo, and is consistent with the best retrieval snow albedo in Greenland in summer months.

Model Sensitivity Analysis
We now consider the sensitivity of the bias-adjustment procedure to choices made in its formulation. First we examine the behavior of different candidates for the reference albedo in Equation 3 ( Figure 5); then we examine the effects of different reference albedos on the adjusted albedo at Summit ( Figure 6); last, we demonstrate the sensitivity of the corrected albedos to different C values ( Figure 7).
Recall that the adjustment algorithm (2) employs as the reference albedo (  A 0 ) either the zonal mean albedo at SZA = 55°, or the (possibly adjusted) zonal mean albedo at latitude = 63°N. Figure 5 compares the temporal variation of the reference albedo (red curve) employed in 2005 to the unadjusted zonal mean albedo at latitude 63°N (green curve). The temporal variation of these albedos are quite similar, and are much the same for both WSA and BSA in 2005 (other years, not shown, are similar). The reference albedo employed is usually (especially during warm months) higher than the unadjusted mean albedo at 63°N. The advantage of a reference albedo that migrates in space and time (i.e., with SZA = 55°), over one fixed to a certain latitude (i.e., 63°N) is that the former puts more weight on the MODIS albedo retrievals located closer to those to be adjusted, and this better retains the spatio-temporal variations captured by the raw MODIS albedo.
The reference albedo is always calculated from from the best quality (Q=0) retrievals. Between 10 and 350 of these Q=0 gridpoints compose each  A 0 ( Figure 5). The number of Q=0 gridpoints used to calculate the reference albedo decreases in warm months when the latitude where SZA=55° migrates to northern Greenland, and in cold months when this latitude migrates to southern Greenland which is relatively narrow geographically.
The corrected albedo at Summit is affected by the choice of reference albedo ( Figure 6). In warm months when Summit's SZA < 55°(i.e., days 129-217), the WSA is close to the BSA and no adjustments are made even though the MODIS albedo is about 10% less than the GC-Net measurements. Raw WSA and BSA both reach maxima when SZA is around 55° on days 129 and 217/225, and then artificially decrease as SZA increases. Moreover, for SZA larger than 55°, WSA decreases faster with SZA than does BSA. Our algorithm, Equation (2), adjusts the artificial albedo decreases before day 129 and after day 217. The adjusted albedo is larger than 0.8 and remains constant or increases slightly with SZA. On days 105 to 121 and days 225 to 241, the adjustment uses the reference albedo at SZA=55°, and the adjusted albedo varies more than it would have had the unadjusted albedo at 63°N been used as the reference albedo. Before day 105 and after day 241, the adjustment algorithm uses the albedo at 63°N as the reference albedo. Equation (3b&c) remediates the artificial decrease of the reference albedo at 63°N and this causes the adjusted albedo increase slightly with SZA before day 89 and after day 273.
The C parameter in Equation (2) (2), increases along this transect by 0.03 for WSA and 0.02 for BSA as SZA increases from 65° (at 65°N) to 80° (at 80°N). The adjusted albedos increase slightly, rather than dropping sharply, for SZA > 71°. Raw albedos are only adjusted where they are less than the daily reference albedo (0.807 and 0.815) and, simultaneously, less than 0.80. to -0.05 and -0.08 between the original MODIS albedo and GC-Net observations. Here "mean difference" is the arithmetic mean of the difference between the MODIS albedo and the GC-Net measured albedo. As mentioned above, GC-Net overestimates snow albedo by about 0.03 relative to a limited collection of more precise measurements at Summit (Stroeve et al, 2005). Hence the adjusted MODIS snow albedo which is ~0.03 less than GC-Net measurements is consistent with those more precise measurements. The adjustment is a conservative estimate of snow brightness in Greenland in that it could be up to 0.04 brighter than the reference albedo without exceeding GC-Net measurements.

Adjusted MODIS Albedo
Nonetheless we err on the side of the minimal adjustments necessary to bring MODIS into agreement with other measurements and with theory at large SZA.
The geographic distribution of raw and adjusted albedos on two days in March (day 73) and April (day 105) 2005 demonstrates how the adjustment compresses the dynamic range of albedo, while retaining the observed structure, and imposing the theoretically predicted and in situ-observed brightening of snow with the increasing SZA ( Figure 9). The adjusted snow albedo also preserves the spatial continuity with the raw albedo for SZA less than 55° around the latitude of 65°N on day 105 ( Figure 9d). On these two days, the raw MODIS snow albedo generally decreases as latitude and SZA increase. The raw values are physically unrealistic (beneath 0.75) over large regions of Greenland.

Surface Energy Budget
The increased brightness of the adjusted MODIS snow albedo (BSA) reduces climatological mean (average of 2003 to 2005) monthly ASR by a minimum of 0.7±0.1 W/m 2 (in June) to a maximum of 6.2±0.9 W/m 2 (in September), where the uncertainties are indicated as plus/minus one δ, the standard deviation of the three years comprising the climatology (Table 2) In terms of surface processes in Greenland, the MODIS albedo and implied ASR biases have implications for other terms in the surface energy budget (SEB). MODIS albedos are often used to evaluate or constrain modeled albedos (e.g., Zhou et al., 2003;Oleson et al., 2003). Models that agree with a biased albedo from measurements must compensate with energetically equivalent and opposite biases in the non-constrained terms of the SEB. Simulations with the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM) (Flanner et al., 2007) suggest that, on average, the net longwave emission from Greenland remains relatively constant through the year and that the SEB fluxes most likely to compensate for seasonally varying ASR biases are those of latent heat (i.e., the net snow pack sublimation, snow melt) and sensible heat.
To illustrate the effects that the MODIS albedo and implied ASR biases would cause in a model, we show the snow melt and sublimation that would occur due to the MODIS ASR bias ( Figure   11 and Table 2) being fully converted to melt or sublimation. The MODIS albedo has the largest bias in months from November to February, yet the implied ASR bias and the consequent snow phase change is small because of the limited surface solar insolation in Greenland in winter. The strongest potential impacts of the ASR bias on snow phase change occur from March to May and from August to October for both WSA and BSA. The maximum annual potential snow melt due to this ASR bias is 40.2±1.5 cm and 28.9±1.8 cm snow water equivalent (SWE) for WSA and BSA, respectively. This change is probably closer to the potential snow sublimation between the BSA and WSA because the blue-sky snow albedo combines BSA and WSA, and because most of Greenland is below freezing when SZA is larger than 55°. The

Discussion and Summary
MODIS albedo retrieval algorithms generally work well for the full inversion data (Q=0) at low SZA, although the biases in comparison to in situ measurements are worse for snow than for other surface types (Jin et al., 2003b;Wang et al., 2004;Rosech et al., 2004;Stroeve et al., 2005;Liu et al., 2009). At large SZA, however, the negative bias of the shortwave MODIS albedo worsens with increasing SZA (Liu et al., 2009, Wang andZender, 2010). As a result, MODIS reports snow surfaces too dark by up to 30% in absolute albedo as compared to CERES although the CERES may report  (Zhou et al., 2003;Oleson, 2003;Rosech et al., 2005).
The difficulty of retrieving high quality surface albedos increases with solar zenith angles.
More recently, our analysis of MODIS albedo (MCD43C3) in the dry snow covered regions of Greenland (Wang and Zender, 2010) documents a physically unrealistic snow albedo decline for SZA > 70° for the best quality (Q=0) data, and for SZA > 55° for lower quality (Q > 0) retrievals, also called magnitude inversion retrievals. For these reasons, the MODIS surface albedo science team considers the best retrieval data for SZA > 70° to be suspect, and recommends avoiding use of all magnitude inversion retrievals (Q>0) and of the best retrieval (Q=0) data for SZA > 70°. However, the dearth of in situ measurements means that satellite-based estimates are necessary to characterize surface albedo and its changes in remote polar regions like Greenland.
The strategy of this study is to enhance the MODIS poor quality data by using model simulation and in situ data from the GC-Net to characterize snow albedos at large SZAs in Greenland and hereby remediate the implied bias in the surface energy budget. The adjustment is based on a semi-empirical parameterization of the surface albedo dependence on SZA (Dickinson,1983;Briegleb et al., 1985;Wang et al., 2005). The adjustment hinges on a reference albedo defined, to ensure robustness, as the mean snow albedo either at SZA=55° or at latitude 63°N (in Southern Greenland). The dynamic snow reference albedo relies on nearby gridpoints at the same SZA, and preserves the original temporal variation and spatial continuity with albedo retrievals for SZA < 55°. The various thresholds in Equation (3)  adjustments in regions influenced by wet snow, vegetation, or exposed rock. However, it is very challenging to determine which MODIS gridpoints have wet or dry snow or melt ponds using the MODIS snow albedo value alone. The influences of some small (several meters) and shallow melt ponds may be negligible for the MODIS (MCD43C3) grid size of 5.6 km x 0.8-2.8 km in Greenland.
Here we rely on a threshold albedo of 0.5 to exclude melt ponds, wet snow, or exposed rocks. In addition, snow surface temperature can indicate whether snow is wet or dry. Hall et al. (2009) use -1°C land surface temperature to distinguish "melt" snow from dry snow. The noontime SZA is also anticorrelated with maximum diurnal temperature. When SZA is larger than 55°, the 2-m air temperature in Greenland is usually far below 0°C according to GC-Net measurements. Moreover, the snow surface air temperature is usually much lower than the 2-m air temperature. Thus, the albedo (> 0.5) and SZA (> 55°) thresholds both help ensure that our adjustment is applied only to dry snow gridpoints.
The adjusted MODIS snow albedo behaves in a manner consistent with model simulations, GC-Net measurements, and CERES surface albedos. It has a smaller magnitude of spatial variations than the raw MODIS albedo and, to a lesser extent, than the CERES surface albedo. Compared to GC-Net measurements, the corrected MODIS snow albedo is more realistic than the raw MODIS snow albedo and, at latitudes north of 75°N, than the CERES albedo. The raw MODIS albedo derived ASR bias of 2.9-4.5 W/m 2 exceeds the estimated present day direct climate forcings in Greenland due to anthropogenic agents such as CO 2 (1.4 W/m 2 ) and black carbon (BC) aerosol deposition (0.6 W/m 2 ) (Koch and Hansen, 2005;Flanner et al., 2007& 2009). We intend one consequence of our improved estimates of Greenland surface albedo from MODIS retrievals to be an increased signal-to-noise ratio albedo is best attempted in regions of smaller MODIS albedo uncertainty, i.e., snowy regions of relatively low SZA, and therefore strong surface insolation, such as the Tibetan Plateau, Colorado Plateau, Rocky Mountains, and western European Alps (Flanner et al., 2007;Painter et al., 2007;Ming et al., 2009).
The combined data (comprising low-and high-quality retrievals) MODIS large SZA albedo and implied ASR biases have implications for the other terms in the Surface Energy Budget (SEB) in Greenland. The snow melt and sublimation that would occur due to the MODIS ASR bias being fully converted to melt or sublimation are greatest from March to May and from August to October for both WSA and BSA. The actual effect of the ASR bias on snow phase change is probably closer to the potential snow sublimation between the BSA and WSA because the blue-sky snow albedo combines BSA and WSA, and because most of Greenland is below freezing. Overestimating snowpack ASR can also be expected to trigger positive snow-albedo feedbacks such as the snow cover-albedo feedback and snow albedo reductions due to grain growth and atmospheric conditions (Flanner and Zender, 2006). Simulations (Flanner et al., 2007) with CAM estimate that about 34% of the annual surface net solar radiation in Greenland converts to latent heat. Thus, we expect that models constrained to agree with the raw MODIS ASR would overestimate snow phase changes by about one third of the maximum snow melt or sublimation shown in Figures 11 and 12.
Our empirical adjustment to MODIS albedo is based on snow optical properties and in situ measurements over purely snow-covered regions on Greenland. Further examination utilizing in situ data is necessary to test whether systematic adjustments would improve MODIS snow albedo in other regions with large SZA such as Antarctica, and in partially vegetated regions like North America and Eurasia. Other regions may require different formulations for an adjustment algorithm due to regional variations in snow condition and structure (e.g., Sturm et al., 1995).          31   739  740  741  742  743  744  745  746  747  748  749  750  751  752  753  754  755  756  757  758  759  760  761  762  763   764  765  766  767  768  769  770  771  772  773  774  775  776  777  778  779  780 zenith angle (SZA) in degrees. The albedo adjustments are applied before day 97 and after day 249 at Saddle, and before day 129 and after day 217 at Summit. The GC-Net measured snow albedos shown have a mean positive bias of ~0.03 compared to more accurate Kipp and Zonen CM 21 pyranometer measurements (not shown). Figure 9. Raw (a, b, raw) and adjusted (c, d, new) MCD43C3 shortwave black-sky albedo in Greenland on days 73 (March) and 105 (April), 2005. The SZAs at latitude 65° and 75° on day 73 are 65° and 75°, respectively, and on day 105 are 55° and 65°, respectively. Differences were constructed from gridpoints with both MODIS and CERES albedo > 0.5 to exclude non-snow covered areas. Figure 11. Seasonal cycle of the potential maximum snow melt and sublimation in snow water equivalent (SWE, cm) that would be produced by converting the MODIS absorbed solar radiation bias (raw-new) in Greenland to melt or sublimation, respectively. The vertical whiskers on the SWE_melt curve show the range of ± one standard deviation (δ) uncertainty for each month computed from the three years (2003)(2004)(2005). Surface solar insolation is obtained from CERES estimates. The right handside vertical axis is for albedo.  32   782  783  784  785  786  787  788  789  790  791  792  793  794  795  796  797  798  799  800  801  802  803  804  805  806  807  808   Table 1. Greenland Climate Network (GC-Net) Automatic Weather Stations (AWS) used in this study (Steffen and Box, 2001).   36   856  857  858  859  860  861  862  863  864  865  866  867  868  869  870  871  872  873  874  875  876  877  878  879  880  881  882  883  884  885  886 Figure 3. Correlation coefficient R (blue columns) between cos(SZA) and the zonal mean MODIS black-sky snow albedo (all retrievals) in Greenland from 50°W-40°W averaged into 1° latitude bands from 65°N to 80°N. Filled blue columns represent significant correlation at 0.95 confidence level. Red columns are the albedo difference (alb_65N minus alb_80N) between zonal mean albedo at latitude 65°N and at 80°N. The two green curves show the noontime SZA (using the right hand vertical axis) at latitudes 65°N and 80°N.  . Seasonal cycle of the potential maximum snow melt and sublimation in snow water equivalent (SWE, cm) that would be produced by converting the MODIS absorbed solar radiation bias (raw-new) in Greenland to melt or sublimation, respectively. The vertical whiskers on the SWE_melt curve show the range of ± one standard deviation (δ) uncertainty for each month computed from the three years (2003)(2004)(2005). Surface solar insolation is obtained from CERES estimates. The right handside vertical axis is for albedo.