2.1.

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.

Table 1

The list of MODIS albedo production periods (day of year) used to make the monthly albedo climatology (day-of-year 361 comes from previous year).

2.2.

MODIS Land Cover Product

The yearly MODIS IGBP maps¹⁹ 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.

2.3.

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.²⁰ 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.

2.4.

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 pixel-specific 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.²¹ 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 pixel-scale 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.

3.1.

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. 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).

Table 2

List of the selected multiscale layers and their statistics over all land cover types globally [“Min. req. pixels” represents the minimum required number of samples at the corresponding spatial resolution in order to generate statistics. “Percent valid” represents the percentage of cells that have a sufficient number of valid values and are included in the look-up map (LUM)].

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.

3.2.

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 ($\sim 1\mathrm{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 $\sim 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., $\mathrm{stdev}\_\mathrm{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, $\sim 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.

3.3.

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):

The standard deviation ($\sigma$) of the modeled albedo can be computed using Eq. (2):

3.4.

Computing Albedo from BRDF LUM

Different albedo values may be obtained by assuming purely direct illumination (directional-hemispherical 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 black-sky 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 high-latitude 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.

4.1.

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 snow-free 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.

Fig. 3

An illustration of albedo of mixed forests from albedo LUMs (shortwave broadband white-sky albedo) in August under snow-free conditions (a) and in March under snow-covered conditions (b) at 0.05 deg or coarse resolution (pixel resolution varies for regions). Standard deviations of two albedo LUMs for mixed forest are shown in (c) and (d), respectively. The gray shows the regions that use global average values since mixed forests were not available in a wide range (10 deg).

4.2.

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.

Fig. 4

An illustration of BRDF parameters for mixed forests from BRDF LUMs (shortwave broadband) in March under snow-covered conditions. The window size for the statistics of isotropic (a), volumetric (b), and geometric (c) parameters is shown in (d). The gray in (a) to (c) and the red in (d) represent global statistics.

4.3.

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 white-sky albedo for shortwave broadband (0.3 to $5.0\mu \mathrm{m}$) using the 0.05 deg albedo LUM under snow-free [Fig. 5(a)] and snow-covered conditions [Fig. 5(d)]. Comparison of the snow-free estimate to the retrieved snow-free MODIS albedo [Fig. 5(c)] from the MODIS CMG product (0.05 deg) of the July composite (median value from day 177, 185, 193, and 201) indicates that the differences between the reconstructed albedo [Fig. 5(a)] and actual albedo [Fig. 5(c)] are very small. In addition, the reconstructed albedo map gives a more complete coverage globally. 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).

Fig. 5

Reconstructed white-sky albedo (top) and standard deviation (middle) for shortwave broadband using LUMs and the 2005 MODIS International Geosphere-Biosphere Programme (IGBP) land cover map under snow-free (left panel, July 2005) and snow-covered (right panel, March 2005) conditions. The MODIS albedos under snow-free and snow-covered conditions are shown in the bottom line. The modeled albedo under snow-covered conditions assumes 100% snow cover globally.

Table 3

Mean bias differences (MBD, reconstructed/modeled minus MODIS), mean absolute differences (MAD), and coefficient of determination (R2) between LUM reconstructed albedo and the MODIS climate modeling grid albedo product (white-sky albedo, shortwave broadband, 0.05 deg) in 2005 under snow-free and snow-covered conditions.

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 $\mathrm{MBD}=-0.002$, $\mathrm{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–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.

4.4.

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.

Fig. 6

Modeled white-sky albedo for shortwave broadband at 1 deg resolution using albedo LUM (1 deg) and MODIS IGBP and land-use harmonization IGBP map under complete snow-free [(a) and (b)] and snow-covered [(c) and (d)] conditions in March.

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 $\pm 0.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.

Table 4

Globally averaged monthly and annual mean white-sky albedo in shortwave broadband for land derived based on the MODIS International Geosphere-Biosphere Programme (IGBP) and LUH IGBP map and the albedo LUM at 1 deg resolution under snow-free and snow-covered conditions.

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.²¹ 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.²¹