Multiscale climatological albedo look-up maps derived from moderate resolution imaging spectroradiometer BRDF/albedo products

Abstract Surface albedo determines radiative forcing and is a key parameter for driving Earth’s climate. Better characterization of surface albedo for individual land cover types can reduce the uncertainty in estimating changes to Earth’s radiation balance due to land cover change. This paper presents albedo look-up maps (LUMs) using a multiscale hierarchical approach based on moderate resolution imaging spectroradiometer (MODIS) bidirectional reflectance distribution function (BRDF)/albedo products and Landsat imagery. Ten years (2001 to 2011) of MODIS BRDF/albedo products were used to generate global albedo climatology. Albedo LUMs of land cover classes defined by the International Geosphere-Biosphere Programme (IGBP) at multiple spatial resolutions were generated. The albedo LUMs included monthly statistics of white-sky (diffuse) and black-sky (direct) albedo for each IGBP class for visible, near-infrared, and shortwave broadband under both snow-free and snow-covered conditions. The albedo LUMs were assessed by using the annual MODIS IGBP land cover map and the projected land use scenarios from the Intergovernmental Panel on Climate Change land-use harmonization project. The comparisons between the reconstructed albedo and the MODIS albedo data product show good agreement. The LUMs provide high temporal and spatial resolution global albedo statistics without gaps for investigating albedo variations under different land cover scenarios and could be used for land surface modeling.

Multiscale climatological albedo look-up maps derived from moderate resolution imaging spectroradiometer BRDF/albedo products

Introduction
Land cover changes affect the global climate through both biochemical and physical mechanisms. Biochemical mechanisms include the release of carbon and nitrogen compounds through burning and decomposition. Physical mechanisms include changes in surface albedo, evapotranspiration, and the release of aerosols through burning. The relative strengths of these factors depend on both the type and location of the land cover change or disturbance. 1 Historical land use has long been recognized as a potential driver for global climate change, but the uncertainty of the net radiative forcing arising from land use change is quite high. 2 A better and more accurate description of global albedo, both present and historical, is required to reduce the uncertainty in land use radiative forcing. Moderate resolution imaging spectroradiometer (MODIS) bidirectional reflectance distribution function (BRDF)/albedo products have been available since 2000. 3 They provide necessary information for evaluating the spatial and temporal variability of surface albedos. The variability of albedos depends on the land cover type, geographic location, season, and snow cover as revealed from both field and space measurements. [4][5][6][7][8] Most climate and ecosystem models currently rely on the land cover types present in a map grid cell to compute or deduce land surface albedo. 6,9 Each land cover type is associated with seasonal albedo values through a look-up table (LUT). Therefore, the accuracy of the land cover albedo LUT directly affects the modeling results. Different land surface models may use different land cover classifications and need different albedo LUTs. For example, Carrer et al. built a global vegetation and bare soil albedo map using the 10-year MODIS data products. 10 In the past, we have examined the inter-and intraannual variability of albedo for 17 International Geosphere-Biosphere Programme (IGBP) classes under snow-free and snow-covered conditions. 7 The work is based on the MODIS coarse resolution (0.05 deg) climate modeling grid (CMG) albedo product before 2004 and albedo statistics are summarized in 10 deg latitude bands. 7 However, these albedo LUTs are either too coarse in spatial resolution or too abstract in surface types for evaluating albedo variations that are caused by various land cover changes at the local scale. In addition, the land cover examples derived from MODIS at 0.05 deg resolution are often not pure, but are actually mixed pixels. This is particularly problematic for human-dominated environments, which may be heterogeneous at fine scales.
To evaluate albedo variations caused by land cover changes, a high-quality look-up map (LUM) at a fine spatial and temporal resolution is required. This kind of LUM should include albedo values for each land cover type at a fine spatial resolution. For modeling past and future land cover conditions, it is also desirable to have albedo values for all land cover types available globally. However, albedo data for certain land cover types do not always exist, depending on present land cover distributions, for example, the albedo for evergreen broadleaf forests may be available at tropical regions and not available at high-latitude regions. The availability of MODIS data products is also restricted by cloud contamination and data acquisition strategies. In addition, seasonal LUMs are required to capture albedo values at the different vegetation growth stages. Snow coverage brings the biggest change for the surface albedo and should be considered separately in the LUMs. Uncertainties of albedo from LUMs are needed for better characterization of the albedo changes. These requirements have not been fully considered in past work. 7 This paper presents an improved approach for building a multiscale hierarchical BRDF and albedo LUMs using a decade of MODIS BRDF/albedo products and the Global Land Survey (GLS) collection of Landsat imagery. Under the multiscale hierarchical data structure, BRDF and albedo LUMs were organized in multiple layers. Each LUM layer was generated at a specific spatial resolution. The principal improvements associated with this product are as follows: • Spatially explicit (geographic) characterization of albedo for each IGBP type (i.e., grassland in northern Eurasia can have a different albedo from grassland in South America). • Monthly temporal resolution derived from MODIS climatology.
• Pure examples of each IGBP type, based on Landsat spectral assessment. • Spatially hierarchical representation of albedo, such that the finest possible resolution is captured for all locations and a value for each IGBP type is available for every location on the globe.
The climatological albedo LUM results were applied to yearly MODIS land cover maps and compared to the MODIS albedo product to assess the albedo LUM. To extend the application of the albedo LUM beyond the MODIS era, the current land use map from the Intergovernmental Panel on Climate Change (IPCC) land-use harmonization (LUH) project 11,12 was used to estimate albedo and was evaluated with the MODIS albedo product. The agreement of MODIS albedo and the LUH land use map derived albedo during the MODIS era will support our ongoing investigation of albedo variations for the past and the projection of potential albedo changes for the future under the land use scenarios constructed from the LUH project. The main objective of this paper is to build a high-quality albedo abstraction at high temporal and spatial resolutions for evaluating climate forcing changes under different land use scenarios. We start by describing the required data products and the approach and then present BRDF and albedo LUM results, followed by concluding remarks. Our results rely on the high-quality retrievals of MODIS BRDF/albedo products. The MODIS BRDF/albedo products were extensively validated previously 8,13-18 and will not be discussed in this paper.

Data Preparation
The albedo and BRDF LUMs were built based on a decade of MODIS data products, including the BRDF/albedo products and the land cover product. The required data products for this paper are described below.

MODIS BRDF/Albedo Products
The global gap-filled Terra and Aqua combined MODIS BRDF/albedo products (MCD43GF, geographic projection, 30 arcsec pixel resolution) from 2001 to 2011 were used to build albedo climatology. The 16-day MCD43GF products (eight-day overlapping) were grouped to build monthly composites. Each month was composited of four production periods ( Table 1). The median albedo from four periods during 2001 to 2011 (one month has four production periods * 11 years' inputs in total) were selected to make the monthly climatology albedo composites. Multidate MODIS data products use the first day-of-year (DOY) to name the production period. For the 16-day MODIS BRDF/albedo product, the production period covers dates from day 1 to day 15. In Table 1, the DOY 361 for January composite comes from the previous year, which covers the period from about December 27 to January 11. The monthly climatology albedo represents the less-noisy median values from a decade of high-quality BRDF/albedo retrievals at a fine spatial resolution.

MODIS Land Cover Product
The yearly MODIS IGBP maps 19 from 2001 to 2011 were used to generate a multiyear composite IGBP map by selecting the majority classes during the period. The yearly MODIS IGBP map is a major land cover product from MODIS land cover products. The 17 IGBP surface types include water, evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest, mixed forests, closed shrubland, open shrubland, woody savannas, savannas, grasslands, permanent wetlands, croplands, urban and built-up, cropland/natural vegetation mosaic, snow and ice, and barren or sparsely vegetated. Our composite IGBP map represents the most representative IGBP class during the period. Both the composite IGBP map and the climatology BRDF/albedo at a 30 arcsec pixel resolution were used to create albedo and BRDF statistics. The yearly MODIS IGBP map was used to compute and examine the reconstructed albedo from the albedo LUMs.

Landsat Data
The GLS collections of Landsat images from 2000 and 2005 were used to measure the subpixel homogeneity of MODIS data products. The GLS collection includes a global dataset of the 30-m resolution Landsat imagery (over 8700 scenes from GLS 2000 and over 9700 scenes from GLS 2005) to support land cover and land use change studies. 20 Most Landsat images were acquired during the growing season, which is ideal for detecting surface homogeneity for the growing season. However, they may limit some MODIS homogenous samples during the nongrowing season in the vegetated areas.

LUH Land Cover Map
Land cover maps from the LUH project were used to evaluate MODIS albedo LUMs. The LUH project generates a harmonized set of land use scenarios that smoothly connect historical reconstructions of land use with future projections. The LUH product describes spatially explicit fractional land use patterns and annual transitions for both historical (1500 to 2005) and future (2005 to 2100) scenarios at a 0.5 deg resolution within five classes: croplands, urban, pasture, primary and secondary lands, and water/ice. 11,12 These five classes were converted to different proportions of IGBP 17 land cover classes based on a pixelspecific map intersection at 1 deg with the 30 arcsec 1992 to 1993 AVHRR IGBP product of Eidenshink and Faudeen (1994) for the period of overlap, as described in Ghimire et al. 21 LUH cropland and urban classes were directly mapped to the same IGBP classes, while pasture was mapped to shrublands, savannas, and/or grasslands, and primary and secondary lands were mapped to forests, shrublands, savannas, and/or grasslands, each based on pixelscale proportions. This mapping rule was assumed to be temporally invariant and was applied across time. The resulting product provides spatially explicit dataset on land cover change with detailed land cover classes for evaluating albedo changes. In this paper, we use MODIS albedo LUMs and the LUH-derived IGBP map to compute land surface albedo and then compare the synthetic albedo map to that derived from using the MODIS IGBP map at a 1-deg resolution.

Approach
The base (finest) spatial resolution for the multiscale LUMs in this paper is defined as 0.05 deg (the same as the MODIS CMG products). Ideally, the BRDF and albedo LUMs require statistics (mean and standard deviation) for all IGBP classes at the base resolution. However, given the current geographic distribution of land cover classes, in most cases, there are not enough MODIS samples to generate robust statistics for every land cover class for all 0.05 deg cells around the globe. In order to build an LUM that includes the statistics for every land cover class, a multiscale approach was developed. Multiple layers of albedo statistics at different spatial resolutions for different IGBP classes were generated. This data structure allowed the generation of a complete LUM by combining multiple LUM layers from different spatial resolutions such that each value for a different land cover type can be appropriately described. Figure 1 illustrates the data flow chart for generating and evaluating albedo LUMs. The gap-filled MODIS albedo products (MCD43GF) from 2001 to 2011 were used to make monthly climatology albedo composites. The GLS collections of Landsat images were used to generate the map of pure MODIS pixels through a MODIS homogeneity test. The gap-free albedo LUMs at multiple spatial resolutions were generated from climatology albedo composites and the pure MODIS pixels based on multiscale hierarchical albedo statistics. The albedo LUMs were used to compute albedo based on the MODIS IGBP land cover map and the LUH land use map. The resulting albedos were compared to the standard MODIS albedo products. The application of the modeled albedos for 1500 to 2100 based on the LUH land use maps are presented and analyzed separately in Ghimire et al. 21 The major steps in Fig. 1 are discussed in the following subsections.

Multiscale Hierarchical Layers
The multiscale LUM data structure used hierarchical layers to build a complete LUM (Fig. 2). For each LUM layer, albedo and BRDF statistics were generated for a given land cover class if there were enough samples. Table 2 lists the selected layers we generated and the minimum required numbers of pixels for generating the statistics over all IGBP classes. At the base resolution (0.05 deg CMG resolution), the minimum required numbers of pixels for an IGBP class was defined as 2 (of the total of 36 pixels in a 0.05 deg cell) so that more spatial detail could be reserved. At a coarser resolution, the required numbers of pixels were raised (>20) so the statistics are more robust. The sample size of 20 allows finding an estimate of the population standard deviation within 30% of the true value with a confidence level of 95% for a normal population. In the work, we included spatial resolutions at 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7.5, 10, and 360 deg. The percentages of valid statistics for all IGBP classes (land surface) increase as the cell spatial resolution decreases. These percentages were computed based on the total numbers of global pixels including both water and ocean, thus the percentages (mostly land surface) are relatively low. For the coarsest global layer (360 deg, corresponding to a single, global LUT), statistics are available for every IGBP class under snow-free conditions, but some IGBP classes (e.g., evergreen broadleaf forest) still lack sufficient samples under snow-covered conditions. Fig. 1 Data processing flow chart for generating and assessing albedo look-up maps (LUMs) (rectangles represent data; italic texts describe processing; gray rectangles show input data products). Fig. 2 The data structure for LUM layers at different spatial resolutions. When the BRDF/albedo value for a pixel is not available (blank small rectangles), a substitute value from the upper (coarser resolution) layer will be used if it is available (represented by solid small rectangles).
The LUM layers have many gaps (missing values) as shown in Table 2. The number of missing values decreases as the layer spatial resolution decreases. When the layer resolution reaches the global scale (360 deg cell), statistics for every land cover type are normally available except under some extreme cases. Based on the LUM layers, we can build a complete LUM at a given spatial resolution in the multiscale layers. Starting from the given layer, if statistics for a pixel are available from the starting layer, the LUM will use them. Otherwise, the algorithm will search coarser layers until the statistics for the land cover type become available (Fig. 2). In the worst case, the algorithm may end up with the global statistics (same as global LUT). In this way, we can preserve available spatial details and also maintain a complete LUM for all land cover classes.

MODIS Homogeneity Test
In the albedo LUM, a basic assumption is that each IGBP land cover class should associate with an albedo value. Therefore, we need to ensure the albedo statistics are generated from these pure MODIS pixels with a dominate land cover class. However, at the 30 arcsec (∼1 km) pixel resolution, a MODIS pixel may represent the mixture of multiple IGBP classes, and the IGBP class assigned to this pixel only represents the majority class. The albedo/BRDF statistics computed from the mixed pixels could be biased. To limit the use of the mixed pixels, we used Landsat imagery to quantify class homogeneity at the MODIS resolution, thus to derive statistics from relatively homogeneous samples. The 2000 and 2005 GLS Landsat datasets were used to generate indicators of homogeneity. Each Landsat scene was first reprojected and aggregated from 30 m Landsat resolution to 30 arcsec MODIS resolution. There are ∼1000 Landsat pixels included in each 30 arcsec MODIS pixel. The global map of the standard deviation of NDVI at a 30 arcsec pixel resolution was then computed. We estimated the homogeneity of a MODIS pixel using the standard deviation of NDVI (std_ndvi) from Landsat pixels. A smaller standard deviation of NDVI means more uniform Landsat NDVI values exist in the MODIS pixel, thus it is more likely to be homogeneous. After extensive tests and visual interpretation, we selected the threshold of 0.12 (i.e., stdev ndvi < 0.12) as an indicator to choose the most homogeneous MODIS pixels. This simple threshold minimizes obviously mixed pixels (e.g., shorelines and the boundaries of different land cover types) yet still keeps enough pixels for the LUM generation. Note that the threshold is an empirical indicator for filtering out highly heterogeneous pixels in our work and can be adjusted to include or exclude more pixel samples. Only MODIS pixels that pass the homogeneity test are used to compute mean values and the standard deviation. In this paper, ∼17.5% of MODIS land pixels were detected as mixed pixels and were excluded from the LUM computing. Since the GLS Landsat images were mostly acquired during the growing season, we may exclude some MODIS pure samples during the nongrowing season in vegetated areas in the LUM generation.

Computing Albedo Using Albedo LUMs
Using the albedo LUM and land cover map at the designed spatial resolution, we can compute the albedo (a) for each pixel using Eq. (1): where i represents the IGBP land cover class and n is the total number of classes within the pixel; a i is the albedo for class i from the albedo LUM; and p i is the percentage area of class i from the land cover map within the pixel. In the albedo LUM, all a i were derived from the pure and homogeneous class that passed the homogeneous test. Since albedos at a fine resolution can be linearly aggregated to a coarse resolution, Eq. (1) is scalable for different spatial resolutions, which means energy can be conserved at different spatial resolutions. The standard deviation (σ) of the modeled albedo can be computed using Eq. (2): where σ i is the standard deviation of the albedo from the albedo LUM.

Computing Albedo from BRDF LUM
Different albedo values may be obtained by assuming purely direct illumination (directionalhemispherical or black sky), purely diffuse illumination (bihemispherical or white sky), or actual conditions based on cloud cover and atmospheric scattering (actual or blue sky). 3,8 The MODIS albedo product provides the black-sky (direct) albedo at local solar noon. The statistics of blacksky albedo at a very coarse resolution could vary depending on the pixels used due to the variation of solar zenith angles at local solar noon. To compute the black-sky albedo at a different solar zenith angle, BRDF parameters are needed. In addition, the full expression of the actual albedo needs to consider the effects of anisotropic diffuse illumination and multiple scattering. Recent studies show that the isotropic diffuse illumination assumption underestimates the albedo relative to the more realistic full expression albedo for the snow-covered conditions at the highlatitude areas. 14,22 This is mainly due to the longer path length of light in the atmosphere that increases multiple scattering and the proportion of diffuse radiance at large solar zenith angles. 23,24 Here, we have also built snow BRDF LUMs using the full expression of the actual albedo. The procedures to compute the actual albedo using BRDF parameters have been documented in much research. 3,14,22 The snow-covered BRDF LUMs were built in a slightly different way from the albedo LUMs. To keep the BRDF properties from individual pixels at the original spatial resolution, a nearest-neighbor search approach was adopted and combined with the multiscale hierarchical approach. The BRDF parameters from the 30 arcsec MCD43D CMG BRDF products for each month and for each pixel were first averaged from a decade of high-quality BRDF retrievals. The BRDF parameters are preserved if there are valid values from the 30 arcsec pixel. Otherwise, the algorithm will search BRDF parameters for each land cover type within a large window, first at 10 arcmin and then increasing to 20 and 30 arcmin if valid BRDF values are not available. The multiscale hierarchical approach is activated if there are still no valid BRDF parameters that can be located. At that point, the BRDF parameters are filled by the multiscale LUM layers at 1, 5, 10, 20 deg and the single, global layer.

Albedo LUM
Based on the multiscale albedo layers, we can build a complete LUM at a given spatial resolution in the multiscale layers. Figure 3 demonstrates an example of the albedo LUM for mixed forests under snow-free [ Fig. 3(a), August] and snow-covered [ Fig. 3(b), March] conditions. The standard deviations for two albedo LUMs are shown in Figs. 3(c) and 3(d), respectively. Under snowfree conditions, for the regions (e.g., Boreal forest) where the albedo for mixed forests existed, the spatial variations of the albedo were reserved. For the regions (e.g., tropical forest) where mixed forests were rare, the albedo statistics were obtained from a coarse resolution layer and the images show large blocks. For the regions (e.g., ocean shown in gray) where mixed forests were not available in a wide range (10 deg), the global average (0.130 for mixed forests) and standard deviation (0.016) were used. Our albedo LUM for mixed forests in Fig. 3 shows spatial variations of albedo that will be otherwise missed if one LUT value of 0.130 is used. In March under snow-covered conditions, the albedo LUM for mixed forest cannot provide spatial variations for most regions due to the lack of MODIS observations under snow-covered conditions. A global average (0.282) and standard deviation (0.072) were used for these regions (gray). The standard deviation map reveals the variations of MODIS albedo samples for different regions.

BRDF LUM
The BRDF LUMs were built under snow-covered conditions. Figure 4 shows an example of three BRDF parameters (isotropic, volumetric, and geometric) for mixed forest in March. The BRDF parameters were preserved in the North high-latitude band since the original 30 arcsec BRDF parameters existed. Otherwise, the statistics of the same land cover type from a large window were used. The large blocks in Fig. 4 are the values produced from these large windows. The window sizes are shown in Fig. 4(d). The BRDF parameters were filled by single global statistics [gray in Figs. 4(a) to 4(c) and red in Fig. 4(d)] if they are not available in a 20 deg BRDF LUM layer. In this example, the gray area in Figs. 4(a) to 4(c) represents the global values of 0.317, 0.077, and 0.034 for isotropic, volumetric, and geometric parameters, respectively.

Reconstructed Albedo Using Yearly MODIS IGBP Product
The global surface albedo can be computed from the BRDF/albedo LUMs described above combined with a global IGBP land cover map using Eq. (1). To test the realism of our albedo LUMs, we used the 2005 MODIS IGBP land cover map (MCD12Q1, 500 m) to compute the albedo. The fractions of each IGBP class in every 0.05 deg grid cell were generated. Figure 5 shows the     Table 3 lists the comparisons between the reconstructed albedo from LUMs and the albedo from MODIS albedo products for all months in 2005. The mean bias difference (MBD, reconstructed minus MODIS) and mean absolute difference (MAD) are less than −0.002 and 0.01. The standard deviation of the reconstructed albedo is computed using Eq. (2). In Fig. 5(b), the standard deviations for >97% of pixels are <0.02. A higher standard deviation is observed for surface types with a high albedo (e.g., desert).
Similar to the snow-free albedo, snow-covered albedo can be computed using the corresponding albedo LUMs and the IGBP land cover map. Figure 5(d) shows the modeled white-sky albedo for shortwave broadband under snow-covered conditions for March 2005. In computing, we assume snow presents globally for every pixel. The standard deviation of the modeled snow albedo [ Fig. 5(e)] shows a higher uncertainty (red) in the regions that lack the snow albedo data, like deserts and tropics. The monthly median value composite from the MODIS CMG albedo product shows limited areas covered by snow and had the highest BRDF/albedo retrieval quality in March 2005 [ Fig. 5(f)]. In contrast, the LUM modeled albedo provides an estimation of snow albedo for most regions except for tropical forests that lacked snow-covered albedo examples for the LUM. Table 3 shows MBD, MAD, and R 2 between the reconstructed albedo and MODIS albedo product. The reconstructed albedo agrees with the MODIS albedo monthly composite very well under the snow-free condition (annual MBD ¼ −0.002, MAD ¼ 0.009, and R 2 ¼ 0.969). Since the IGBP classes are generally stable during the period from 2001 to 2011 at the MODIS spatial resolution, the MODIS albedo during this period can be reconstructed with the climatological albedo LUMs. Note that the capability of reconstruction of albedo is simply a verification of the albedo LUMs. Albedo validations and comparisons were performed and have been discussed in previous literature. 8,[13][14][15][16][17][18] The snow-covered albedo shows larger differences between the LUM modeled and the MODIS albedo. The MBD and MAD under a snow-covered condition are within the range of (−0.014, 0.005) and (0.014, 0.037), respectively. This is mainly due to annual variations of snow coverage and conditions. The statistics for snow-covered albedo are affected by many factors, such as snow ages, snow depth, snow grain size, the geometrical structure of snow-free surface, etc. The snow-covered albedo LUMs have not included these details. A decade of MODIS data products may not be long enough for generating climate data record, especially under snow-covered conditions. In order to compensate for this, a maximum albedo layer from all months under the snow-covered condition is also produced and saved as a separate LUM.

Modeled Albedo Using Land Use Harmonization Map
The percentages of IGBP types from the LUH map and the albedo LUMs from this work are used to compute the albedo under snow-free and snow-covered conditions. Figure 6 shows the modeled diffuse (white-sky) albedo for shortwave broadband in March using LUH and MODIS IGBP fraction maps at 1 deg spatial resolution. The MODIS IGBP fraction map was derived based on the MODIS era climatology IGBP map. The LUH IGBP fraction maps during this period are very stable. Figure 6 shows the mean albedo values in March for the MODIS era. Both snow-free and snow-covered conditions are assumed globally in March. Although the LUH maps were cross-walked from 5 LUH classes to 17 IGBP classes, the spatial variabilities of the modeled albedo are very consistent between the LUH and MODIS modeled albedos under both conditions. For comparison purposes, Fig. 6 only demonstrates the overlapped areas of the modeled albedo from MODIS and LUH. Table 4 lists the mean albedo modeled from LUH and MODIS IGBP maps for each month and all the months combined. The averaged values are weighted using the actual area of each 1 deg cell, so they are the equal area averages. Under a snow-free condition, the annual mean albedos are very close (0.170 from MODIS and 0.169 from LUH) and the annual R 2 is 0.915. For each individual month, the mean bias difference is within AE0.002. Under snow-covered conditions, the annual mean albedos are also similar (0.536 from MODIS and 0.538 from LUH). The monthly mean albedo and the differences between MODIS and LUH show some seasonal variations. The differences between the two modeled snow-covered albedo maps are larger from June to September than the other months. This is due to the lack of high-quality snow albedo retrievals in the Northern Hemisphere during summer. However, the annual differences of the two modeled albedo maps are small (0.002 in absolute albedo value or <0.4% in relative value) as compared to the generally high values of snow-covered albedo.
The modeled albedo based on albedo LUMs from this paper and LUH land use map from the IPCC LUH project have been generated from 1500 to 2005. Global albedo change and radiative cooling from anthropogenic land cover change were investigated. 21 The study found that the mean annual global albedo was increased due to land cover transitions from natural vegetation to agriculture during 1700 to 2005. 21

Conclusions
This paper builds monthly climatology BRDF and albedo LUMs from a decade of MODIS data products and Landsat imagery using a multiscale hierarchical design. The hierarchical approach preserves the spatial resolution from MODIS data products and fills gaps that are missing from the fine spatial resolution data. The LUMs from this paper provide an abstraction of BRDF/ albedo for 17 IGBP classes at an appropriate spatial resolution for each location. Since the LUMs include albedo and BRDF statistics for each IGBP class for all possible locations, they can be used to generate albedo maps for historic or projected (future) land cover that may be very different from current conditions. This feature allows the investigation of the albedo change trends under different land use and land cover scenarios.
Initial results show that the LUMs can be used to reconstruct MODIS albedo using MODIS IGBP map under snow-free and snow-covered conditions. The application of the LUMs to a land use map from the LUH project reveals close results in comparison to the reconstructed albedo from MODIS land cover. There is no bias on the annual albedo under snow-free conditions and <1% relatively under snow-covered conditions, which are affected by different conditions. The capability of extending MODIS IGBP-based albedo LUMs to the LUH land use map during the MODIS era would allow the investigation of albedo variations for the past and projection of potential albedo changes for the future if the albedo properties of IGBP classes are similar.
The BRDF/albedo LUMs provide spatially and temporally complete statistics of albedo and BRDF for each IGBP class at a global scale, which could be used in land surface modeling studies. The BRDF/albedo LUMs in this paper were built based on the IGBP classification. Tao He is a research associate in the Department of Geographical Sciences at the University of Maryland, College Park. His areas of interest include surface anisotropy and albedo modeling, data fusion of satellite products, and long-term regional and global surface radiation budget analysis.
Zhuosen Wang is a NASA postdoctoral fellow at NASA Goddard Space Flight Center, Greenbelt. His area of interests include monitoring the high latitude snow melt and energy balance, modeling and evaluation of land surface albedo, monitoring vegetation phenology and climate change, estimating canopy structure using optical remote sensing and Lidar data, and simulating carbon and energy cycle and fluxes using a regional land surface model in coastal areas.
Bardan Ghimire is a postdoctoral research fellow in the Climate Sciences Department at Lawrence Berkeley National Laboratory. His research expertise involves understanding how ecosystem processes respond to climate, nutrients (e.g., nitrogen and phosphorus) and landcover change including large-scale land surface feedbacks to climate through changes in carbon, water and energy cycles.
Yanmin Shuai is a research scientist at the Earth Resources Technology Inc., working on NASA GSFC projects. Her research interests are focused on the development of routine direct broadcast anisotropic and radiative products for MODIS as well as approaches for Landsat albedo, the detection and monitoring of phenological events over agriculture and forested regions, and the radiative evolution of terrestrial ecosystems disturbed by fires, harvesting, and insect epidemics.

Jeffrey
Masek is a research scientist in the Biospheric Sciences Branch at the NASA Goddard Space Flight Center. His research interests include mapping landcover change in temperate environments, application of advanced computing to remote sensing, and satellite remote sensing techniques.
Crystal Schaaf is a professor in the Environmental, Earth and Ocean Sciences Department at the University of Massachusetts, Boston. Her current interests include modeling reflectance