An Algorithm for the Retrieval of High Temporal-Spatial Resolution Shortwave Albedo from Landsat-8 Surface Reﬂectance and MODIS BRDF

: Variations in surface physicochemical properties and spatial structures can prominently transform surface albedo which conversely inﬂuence surface energy balances and global climate, making it crucial to continuously monitor and quantify surface dynamics at ﬁne scales. Here, we made two improvements to propose an algorithm for the simultaneous retrieval of 30-m Landsat albedo, based on the coupling of Landsat-8 and MODIS BRDF. First, two kinds of prior knowledge were added to disaggregate BRDF, including the Anisotropic Flat Index (AFX) and the Albedo-to-Nadir reﬂectance ratio (AN ratio), from MODIS scales into Landsat scales. Second, a simpliﬁed data fusion method was used to simulate albedo for the same, subsequent, or antecedent dates. Finally, we validated the reliability and correlations of the algorithm at six sites of the Surface Radiation (SURFRAD) budget network and intercompared the results with another algorithm called the ‘concurrent approach’. The results showed that the proposed algorithm had favorable usability and robustness, with a root mean square error (RMSE) of 0.015 (8%) and a mean bias of − 0.005; while the concurrent approach had a RMSE of 0.026 (14%) and a mean bias of − 0.018. The results emphasized that the proposed algorithm has captured subtle changes in albedo over a 16-day period.


Introduction
Surface albedo is the ratio of the reflected light to the total incident energy of a surface in the hemispherical space [1]. It plays an essential role in global climate change, which affects the Earth's energy budget [2][3][4]. Variations in the physicochemical properties and spatial structure of local objects in turn transform surface albedo [5][6][7][8]. At present, surface albedo has been widely used in modelling for radiation and energy balances, numerical weather prediction, atmospheric circulation, and land surface processes [9][10][11][12][13]. Studies have shown that uncertainties in albedo can cause the results of climate and phenology simulations to deviate greatly from real-world conditions [5,[14][15][16]. To minimize any adverse effects on calculations, surface albedo products are usually required to meet the standards for climate modelling, with an uncertainty of less than 5% [17]. Therefore, there is an urgent need to derive a method for obtaining high-accuracy surface albedo measurements, for monitoring global climate and vegetation phenology.    (Table 1) used in this study were acquired from high-resolution satellite data on https://www.google.com/earth. The red markers and black circles represent the location of the measurement towers and the 126-m diameter area they covered, respectively.  (Table 1) used in this study were acquired from high-resolution satellite data on https://www.google.com/earth. The red markers and black circles represent the location of the measurement towers and the 126-m diameter area they covered, respectively.
The albedo value measured in situ at Landsat local time is the average of the fieldmeasured values of the first half hour and the next half hour of Landsat's transit time. As instantaneous conditions such as clouds or birds cause instantaneous disruptions to the accuracy of pyranometer measurements, noise is included in the analysis of albedo records. Therefore, the measured albedo and the ratio of diffuse light to sky light per minute with unusually large standard deviations are excluded [46]. Additionally, while all in situ field measurements with clear days and no clouds were included, field measurements with cloud cover greater than 30% and diffuse-to-total irradiance ratios greater than 30% at local solar noon (LSN) were excluded. The measurement towers cover the same area as a Remote Sens. 2021, 13, 4150 5 of 24 3 × 3 Landsat grid. For comparisons of the Landsat albedo to the measured in situ values, we also used the cosine-factor-based upscaling method to aggregate the Landsat albedos from 30-m resolutions to the measurement tower area of coverage [37,63].

Satellite Data
In this study, we used the latest version of the daily concurrent MODIS MCD43 V006 product and Landsat-8 OLI surface reflectance product to retrieve albedo for the period between 2013 and 2016. The Landsat-8 data were obtained from the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov), including 85 cloudless scenes over nine different paths/rows ( Table 1). The accuracy of Landsat-8 OLI surface reflectance depends on atmospheric correction; thus, we used the Land Surface Reflectance Code (LaSRC), which is an accurate atmospheric correction algorithm specially designed by USGS for Landsat-8 data [64], for atmospheric correction. The concurrent MODIS BRDF/albedo data (Collection 6) were downloaded from the Land Processes Distributed Active Archive Center of the US National Aeronautics and Space Administration (NASA) and USGS partnership (https://lpdaac.usgs.gov). The MCD43 BRDF V006 product, which has been validated in the third stage [65], uses high-quality cloudless MODIS surface reflectance data within 16-days to retrieve surface anisotropy information through a kernel-driven model [20,61,[66][67][68]. It includes surface BRDF parameters (MCD43A1), BRDF parameter inversion quality (MCD43A2), albedo (MCD43A3), and Nadir BRDF-Adjusted Reflectance (NBAR) (MCD43A4), and encourages users to use frequency-band-specific quality control (QA) to obtain the highest quality complete inversion results [45,61,69]. Therefore, the high quality BRDFs marked 0 were used. Both MCD43 BRDF and Landsat-8 OLI data are available from https://earthdata.nasa.gov.

Theoretical Basis to Retrieve Albedo
The albedo inherent to natural surfaces depend on the anisotropic properties of the land itself, the BRDF, which describes how surface reflectivity changes with the Solar Zenith Angle (SZA) and viewing angle [70]. Currently, the recognized kernel-driven model consists of an empirical, linear combination of three scattering properties, that can accurately characterize the BRDF of surface patches [18,71,72]: Among them, R is the Bidirectional Reflectance Factor (BRF); λ is the given MODIS spectral band; θ v is the observation zenith angle; θ s is the SZA; ϕ is the relative azimuth; θ v and ϕ are set to 0 to simplify the calculation when calculating the BRF for Landsat viewing. F iso , f vol and f geo are the three spectral constants of BRDF, representing isotropic scattering, volume scattering and geometric optical scattering, respectively. K vol is the volume scattering kernel called Ross-Thick [71,72]. K geo is the geometric optical scattering kernel called Li-Sparse-Reciprocal [18,73]. Once the surface anisotropy is clearly represented by the BRDF parameters, the Black-Sky Albedo (BSA, R λ ) is obtained by integrating the viewing angles over the entire hemisphere at any given SZA. The diffuse light is further integrated at all illumination angles to obtain the White-Sky Albedo (WSA, R λ ).
where H is the integral value of the kernel function, please refer to Lucht et al. [18] for specifics. At the same time, if the ratio of diffuse light radiation to the total radiation is given under certain light conditions (S θ s ), combined with the BSA and WSA, the actual (or blue-sky) albedo (R λ (θ s )) for a certain time such as could be obtained at the surface area by satellites for that SZA [18,74,75].
At present, it is hard to directly retrieve finer-resolution albedo with the kernel-driven model, due to the difficulty in obtaining finer-resolution BRDFs. Fortunately, the concurrent approach proposed by Shuai et al. [37] provides an idea for obtaining finer-resolution BRDF. The concurrent approach made two assumptions that the land surface remains invariable over a 16-day period and the AN ratio are similar in the homogeneous regions covered by both MODIS and Landsat pixels (Equation (5)). Then, after unsupervised classification, it finds the representative pixels (i.e., pure pixel) under the spatial resolution of MODIS, and calculates the black-sky and white-sky AN ratio, a λ and a λ (Equation (6)).
where m denotes the spatial resolution of MODIS 500-m, and l represents the spatial resolution of Landsat 30-m. Ω l is the geometry for the Landsat radiance observation. Then, assuming the land surface remains invariable over a 16-day period, the spectral AN ratio of these pure pixels are used to calculate the spectral albedo for Landsat resolution through Equations (7) and (8), according to the frequency band correspondence between Landsat and MODIS sensors (as shown in Table 2).

BRDF Inversion by Prior AN Ratio and AFX
Although previous studies of surface anisotropy [33,76,77] have discovered an association between characteristic BRDF shapes and various homogeneous surfaces, it is unavoidable that the BRDF shapes retrieved from satellite data are frequently inexact because veracity is muted by heterogeneous surfaces and complicated vegetation types and fragmented structures within the same pixel [37]. The lack of classification of pure pixels in the entire Landsat scene has further hampered the acquisition of exact BRDF shapes [39,45,46]. To address this limitation, we propose to assume that, the near-infrared (NIR) band prior knowledge can represent 500-m information, which is roughly the same as 30-m information within.
In the search for suitable prior knowledge that satisfies this hypothesis, the anisotropic information and ground scattering properties were considered. The spectral AN ratio, which is related to anisotropic information and the SZA, has demonstrated the capacity to capture BRDFs under specific solar illumination conditions. In particular, it is the key factor for computing spectral albedo. In this study, the AN ratio was used as prior knowledge for the determination of BRDF. Figure 2A,B shows the spectral AN ratio of the Landsat red band and NIR band, using the spectral conversion coefficients provided by He et al. [78]. The histograms of AN ratios in the red band and NIR band were obtained through unsupervised classification ( Figure 2C,D). It is worth noting that the AN ratio shown in Figure 2A,B is slightly different from the AN ratio under real-world conditions, because the 500-m grids smoothed the 30-m signal. The red band AN ratio is spatially fragmented (Figure 2A), while the NIR band AN ratio is relatively smoother within 500 m ( Figure 2B). At the same time, the NIR band AN ratio of each type ( Figure 2D) is more concentrated than that of the red band ( Figure 2C), which indicates that the NIR AN ratio is more suitable to satisfy the hypothesis. Therefore, the NIR band AN ratio was used as one of the links that connect different scales. The BRDF shapes of the ground features depend not only on the spectral AN ratio but also on the scattering characteristics. Therefore, it is necessary to introduce further a priori knowledge to simultaneously describe the scattering properties of ground objects. The anisotropic flat index (AFX) is a vegetation index that represents a nonlinear response to changes in surface structure [79]. Further classification of AFX can produce BRDF archetypes, which can be effectively used to classify BRDFs based on ground scattering properties [50,[79][80][81]. A study by Jiao et al. [49] showed that the range of variation of AFX in the NIR band is lower than that in the red band, indicating that the NIR band AFX is more stable than the red band AFX. Therefore, the NIR AFX was also used to retrieve the BRDFs of pure pixels, as follows: The NIR AFX depends on two parameters: geometric optics and volume scattering. The BRDF shapes of the ground features depend not only on the spectral AN ratio but also on the scattering characteristics. Therefore, it is necessary to introduce further a priori knowledge to simultaneously describe the scattering properties of ground objects. The anisotropic flat index (AFX) is a vegetation index that represents a nonlinear response to changes in surface structure [79]. Further classification of AFX can produce BRDF archetypes, which can be effectively used to classify BRDFs based on ground scattering properties [50,[79][80][81]. A study by Jiao et al. [49] showed that the range of variation of AFX in the NIR band is lower than that in the red band, indicating that the NIR band AFX is more stable than the red band AFX. Therefore, the NIR AFX was also used to retrieve the BRDFs of pure pixels, as follows: The NIR AFX depends on two parameters: geometric optics and volume scattering. When the volume scattering effect is greater than the geometric optical effect, AFX > 1; when the geometric optical effect is larger than the volume scattering effect, AFX < 1; when the two effects are very close, AFX ≈ 1. AFX reveals the relationship between the surface scattering properties and the BRDF archetypes [49,79].
In summary, we used the NIR AFX and NIR AN ratio as prior knowledge for BRDF inversion, fully utilizing all the high-quality BRDF data. Additionally, the BRDF spatial distribution of each type was more suitable for obtaining the scattering properties and angular characteristics of each location without the use of other auxiliary data. Refer to Section 3.2. for specific operations.

Narrowband-to-Broadband Conversion
The spectral albedo cannot satisfy the research of energy balance on the land surface, it is necessary to further convert the narrowband spectral albedo into broadband shortwave (SW) albedo [57]. Narrowband-to-broadband (NTB) conversion coefficients can be generated by simulating the spectral responses of various ground objects in the digital spectrum library, based on the sensor's extensive radiative transmission [57,82]. This method has been verified to be sufficient for estimating shortwave albedo from spectral albedo [30,57,63,82]. The surface albedo is generated through the conversion coefficients. In this study, we use the conversion coefficients for Landsat-8 OLI developed by Wang et al. [30] under the condition of no snow.

Landsat Albedo Modulation
The assumption that the surface reflectance remains unchanged within 16-day periods is reasonable when the surface changes slowly during the dormant season [37]. However, during the growing season or the period of large-scale changes in land cover (such as fire, snow cover, and melt), this assumption is broken because land cover types change rapidly. Hence, we adopted a multi-temporal fusion method of MODIS-Landsat data based on the MODIS backup algorithm [52], assuming that the MODIS modulation term C represents the change in reflectance of the Landsat wavelength with similar spectra at the Landsat scale. The prediction of Landsat reflectance at t 2 from the Landsat observation data at t 1 is defined as: As the MCD43A1 product provides shortwave BRDF parameters, we simplify the modulation term C. Similarly, we suppose that the MODIS short-wave (SW) albedo change is similar to the short-wave albedo change at the Landsat scale [42]. C SW represents the multiplicative modulation term of the short-wave albedo, which is described as: Then, we can combine Equations (4), (10), and (11), based on the Landsat surface reflectance data on the i-th day and the MODIS data on the j-th day (j ranges from 8 days before i to 8 days after i); finally, retrieve the actual (blue sky) albedo of Landsat on the j-th day: Thus, the calculation procedure can be simplified, based on the existing data, avoiding any impact arising from the disparity between the spectral response functions for the two sensors [42,78].

. Data Preparation
The data processing in this article is shown in the flowchart ( Figure 3). In order to combine MODIS BRDF and Landsat surface reflectance, we need to establish the connection between them. First, a cluster-based unsupervised classification (Iterative Self-Organizing [ISO] data analysis) is used to classify the Landsat-8 data after atmospheric correction (the number of categories is set to [10][11][12][13][14][15]. Then, Landsat-8 coordinate system is converted from the default WGS-84 projection to the SIN projection to associate with MODIS anisotropy information. Next, the percentage of each type of 30-m patch throughout the MODIS grid is calculated. The higher this percentage is, the more applicable the BRDF of coarseresolution pixels is to the finer 30-m pixels. All the 500-m pixels, whose category percentage is higher than the given threshold (60%), are marked as pure pixels [40,83], that is, the representative pixels. At the same time, the percentage is exported as a quality control dataset. Lastly, the highest-quality MODIS BRDF parameters are extracted from the concurrent MCD43A1 BRDF products (the quality is determined by the MCD43A2 QA marks). For those categories that do not have pixels to meet the threshold conditions, the pixels in the top 15% of the percentage in those categories are extracted, and the extracted pixels of those categories are marked as low quality (Table 3).

BRDF Inversion
In the regions with similar spectral types of ground objects, BRDF parameters of all high-quality pure pixels can be extracted by the quality control QA (MCD43A2).  Table 3. The outcomes of BRDF inversion and the corresponding quality labels for each case.

Cases
Representation

BRDF Inversion
In the regions with similar spectral types of ground objects, BRDF parameters of all high-quality pure pixels can be extracted by the quality control QA (MCD43A2). However, the value of BRDF parameters is still dispersed ( Figure 4A). After the band quality control QA in the MCD43A2, we find the high-quality NIR parameters of the BRDF in the MCD43A1 and then calculate the NIR AN ratio (Equation (6)). The WSA after quality control in MCD43A3 and the isotropic BRDF parameters in MCD43A1 are utilized to obtain the white-sky AFX (Equation (9)) of each 500-m patch. To determine BRDF shapes, the consideration is required to get appropriate threshold ranges of the prior AN ratio and AFX. For AFX, we refer to the threshold ranges of six BRDF archetypal parameters by Jiao et al. [49], as shown in Table 4 However, the value of BRDF parameters is still dispersed ( Figure 4A). After the band quality control QA in the MCD43A2, we find the high-quality NIR parameters of the BRDF in the MCD43A1 and then calculate the NIR AN ratio (Equation (6)). The WSA after quality control in MCD43A3 and the isotropic BRDF parameters in MCD43A1 are utilized to obtain the white-sky AFX (Equation (9)) of each 500-m patch. To determine BRDF shapes, the consideration is required to get appropriate threshold ranges of the prior AN ratio and AFX. For AFX, we refer to the threshold ranges of six BRDF archetypal parameters by Jiao et al. [49], as shown in Table 4.   To seek the exact properties that each pixel has, we build a look-up table (LUT) for  To seek the exact properties that each pixel has, we build a look-up table (LUT) for the AN ratio and AFX of each category [63]. Then, each type of pure pixel that satisfies both the AN ratio and the AFX threshold range is a spatial representative pure pixel. The average value of selected pure pixels is taken as the reliable BRDF parameters, and their values ( Figure 4B) are more centralized than before ( Figure 4A). The values of AFX and AN ratio obtained by the new algorithm (dashed line) are no longer fixed to those corresponding to the concurrent approach (solid line), but change with different spatial positions ( Figure 4C). Note that when the inversion fails, the BRDF average value of the successful inversion within 16 days is used to fill, and the pixels of these categories are marked as low quality. All quality marks are output with the BRDF inversion results, see Table 3 for details.

Shortwave Albedo Generation
In the Landsat scene, each 30-m patch has been assigned a daily change BRDF within a 500-m grid. This provides a basis for the estimation of 30-m daily surface reflectance and shortwave albedo under any zenith angle. After obtaining the BRDF parameters of all plots, these BRDF parameters are used to calculate the black-sky and white-sky albedo of each band under the Landsat viewing angle and solar illumination. Then, the spectral AN ratio of all bands is obtained through Equation (6). Landsat data provide the surface reflectance of each 30-m patch, which can be correlated with the spectral AN ratio of each patch to gain the spectral albedo through Equations (7) and (8). Then, shortwave albedo is generated through narrowband-to-broadband conversion, based on extensive radiant transmission simulation. The modulation term C SW can reflect the changes in the albedo of the date that Landsat has not observed. Hence, we use the BRDF parameters of shortwave albedo in MCD43A1 to calculate the actual albedo and gain the multiplication modulation term C SW of daily albedo. Then, C SW is used to transform the albedo that has been retrieved for the day through Equation (12). Finally, the Landsat albedo is generated.

Validation Results at SURFRAD
The results show the time series comparison of Landsat albedo obtained by the concurrent approach and the new algorithm at six sites ( Figure 5A-F). Figure 6 shows the validation results of the concurrent approach product (hereafter 'CAP') and the new albedo product (hereafter 'NP') at each site. In the extended Landsat coverage area, the ground truth was well described by in situ measurements. Landsat albedo retrieval results showed seasonal characteristics, with the highest value during the dormant period and the lowest value when entering a complete growth period.
During growth and dormancy periods, the DRA was uniformly covered by deserts and sparse vegetation, which occupied most regions covered by the Landsat data ( Figure 5A). Similarly, around the FPK and TBL sites ( Figure 5C,E), there was relatively homogeneous land cover, which enabled the retrieval of albedo from satellite data, and the collection of a large amount of effective and coherent anisotropy information (i.e., pure pixels). At three representative sites, the two Landsat albedo retrievals were in good agreement with the field measurement, with positive and negative deviations evenly distributed. The slopes of representative sites were almost within the 5% uncertainty range ( Figure 6A,C,E). Therefore, both albedo products showed high consistencies with field measurements at the representative sites. However, with regard to the non-representative sites of BON, PSU, and GWN, CAP satellite retrieval results poorly represented seasonal characteristics, translating to large deviations during dormant periods and small deviations during growing periods ( Figure 5B,D,F). The uncertainties of non-representative sites were higher than 5%, but the slopes of NP were closer to 1 than the slopes of CAP ( Figure 6B,D,F). Table 5 summarizes the evaluation of the accuracy of the two albedo products relative to in situ measurements at six sites. After excluding cloud-contaminated albedo values, the blue-sky albedo values from NP and CAP against in situ were used for statistics. Uncertainty (root mean square error [RMSE]), accuracy (median deviation), offset (mean bias) and linear correlation (slope and correlation coefficient [R 2 ]) for satellite retrievals were included in the statistical processing [65]. The RMSE, median deviation, and mean bias of CAP at representative sites were 0.015, −0.014, and −0.012, respectively, but showed poor results at non-representative sites, with RMSE of 0.032, median deviation of −0.032, and mean bias of −0.030 (Table 5). However, NP showed a slight improvement at the three representative sites, with RMSE, median deviation and mean bias of 0.011, 0.004 and 0.002, respectively. For the three non-representative sites of satellite retrievals, the RMSE dropped from 0.032 to 0.020, median deviation decreased from −0.032 to −0.019 and the mean bias decreased from −0.030 to −0.015. Therefore, NP had less uncertainty and better accuracy than CAP. During growth and dormancy periods, the DRA was uniformly covered by deserts and sparse vegetation, which occupied most regions covered by the Landsat data (Figure  Table 5 summarizes the evaluation of the accuracy of the two albedo products relative to in situ measurements at six sites. After excluding cloud-contaminated albedo values, the blue-sky albedo values from NP and CAP against in situ were used for statistics.  Although the albedo retrieval for each representative site was in good agreement with the field measurements, the R 2 for some sites were relatively abnormal ( Table 5). The R 2 was relatively poor for representative site, such as DRA (R 2 = 0.39). Conversely, the R 2 may have been higher for each non-representative site, such as GWN (R 2 = 0.83). It is worth noting that R 2 would be very low even if the RMSE is small, if there is a weak linear relationship between satellite retrieval and ground data, or if the amount of data analyzed is small [37,39]. To reduce the impact, the data were analyzed using two groups: representative and non-representative sites, as shown in Table 5. The results of CAP and NP had high correlation coefficients at representative sites, 0.96 and 0.94, and lower correlation coefficients at non-representative sites, 0.77 and 0.78, with overall correlation coefficients of 0.90 and 0.91. Therefore, the correlations between the two products and the ground measurements were similar. Figure 7 shows the compared albedo retrieval results from the two algorithms for all sites and all ground measurements; after aggregation, the total number of samples was 85. At all sites, the slopes from the new algorithm ( Figure 7B) were closer to 1 than those of concurrent approach ( Figure 7A). The overall RMSE decreased from 14% to 8%, the median deviation dropped from −0.028 to −0.011 and the mean bias decreased from −0.018 to −0.005. Therefore, the new algorithm has corrected the underestimation and improved the accuracy of albedo retrieval.  Table 5 for the slopes of the fitted lines. Figure 8 shows the spatial patterns of the albedo in the 15 × 15 km range of the GWN site on 9 April 2014. The albedos from these three products visually had similar spatial  Table 5 for the slopes of the fitted lines. Figure 8 shows the spatial patterns of the albedo in the 15 × 15 km range of the GWN site on 9 April 2014. The albedos from these three products visually had similar spatial distributions in overall patterns, while the MODIS albedo smoothed tiny detail about the surface. Therefore, the two kinds of Landsat retrievals have the potential to finely describe tiny-scale natural surface detail in albedo. Figure 7. Results of albedo retrieval verified on all sites by concurrent approach (A), compared with results of albedo retrieval by the new algorithm (B). Dotted gray lines indicate 5% uncertainty of all sites, solid lines are one-to-one lines. The slopes of the fitted lines were obtained by the least square method, with intercepts of all the fitted lines set to 0. Please refer to the Table 5 for the slopes of the fitted lines. Figure 8 shows the spatial patterns of the albedo in the 15 × 15 km range of the GWN site on 9 April 2014. The albedos from these three products visually had similar spatial distributions in overall patterns, while the MODIS albedo smoothed tiny detail about the surface. Therefore, the two kinds of Landsat retrievals have the potential to finely describe tiny-scale natural surface detail in albedo.  Figure 9 shows the histograms of the three albedo products in Figure 8. Here, the standard deviation (Stdev) was used to quantify the magnitude of spatial heterogeneity between the MODIS product, CAP, and NP. Unsurprisingly, the standard deviations of both CAP and NP were significantly higher than that of MODIS albedo. The standard deviation of the MODIS product was 0.009, while that of CAP and NP were 0.026 and 0.027, respectively, indicating that the two finer albedo products had similar spatial heterogeneity.

Comparison of Three Albedo Products
To test the spatial characteristics of the three albedo products, a random transect line was graphed in Figure 8, from which the data were picked up for spatial characteristics analysis. The results implied that the Landsat albedo retrievals had a better spatial  Figure 9 shows the histograms of the three albedo products in Figure 8. Here, the standard deviation (Stdev) was used to quantify the magnitude of spatial heterogeneity between the MODIS product, CAP, and NP. Unsurprisingly, the standard deviations of both CAP and NP were significantly higher than that of MODIS albedo. The standard deviation of the MODIS product was 0.009, while that of CAP and NP were 0.026 and 0.027, respectively, indicating that the two finer albedo products had similar spatial heterogeneity.
Remote Sens. 2021, 13, x FOR PEER REVIEW 18 expression capability, because they were more dramatic and greater fluctuation than MODIS product ( Figure 10A). Figure 10B further qualitatively analyzes the spatial c acteristics of three products. Obviously, the whiskers of NP albedo had greater range other products, which further indicated that NP had better spatial expression capabi Figure 9. Histograms of the three black-sky albedo products shown in Figure 8. To test the spatial characteristics of the three albedo products, a random transect line was graphed in Figure 8, from which the data were picked up for spatial characteristics analysis. The results implied that the Landsat albedo retrievals had a better spatial expression capability, because they were more dramatic and greater fluctuation than the MODIS product ( Figure 10A). Figure 10B further qualitatively analyzes the spatial characteristics of three products. Obviously, the whiskers of NP albedo had greater range than other products, which further indicated that NP had better spatial expression capability. Figure 9. Histograms of the three black-sky albedo products shown in Figure 8. To test the spatial characteristics of the three albedo products, a random transect line was graphed in Figure 8, from which the data were picked up for spatial characteristics analysis. The results implied that the Landsat albedo retrievals had a better spatial expression capability, because they were more dramatic and greater fluctuation than the MODIS product ( Figure 10A). Figure 10B further qualitatively analyzes the spatial characteristics of three products. Obviously, the whiskers of NP albedo had greater range than other products, which further indicated that NP had better spatial expression capability.

Capture Daily Surface Dynamics
Although a finer resolution of surface albedo has been derived to describe fragmented details, the 16-day Landsat reflectance makes it difficult to capture surface dynamics, especially in the growing and senescent seasons, which are of vital importance in monitoring the seasonal characteristics of climate change. The new algorithm in this study generated a synthetic 30-m high temporal-spatial resolution albedo to solve this problem. In the growing season of GWN, we differentiated between the inversion results from 9 April 2014 and that from adjacent dates to identify daily changes in black-sky albedo in the 15 × 15 km range ( Figure 11). As the new algorithm modulated albedo values, the significant variations in albedo (difference > 0.01) were captured (Figure 11a-l; Figure 11c-p).

Capture Daily Surface Dynamics
Although a finer resolution of surface albedo has been derived to describe fragmented details, the 16-day Landsat reflectance makes it difficult to capture surface dynamics, especially in the growing and senescent seasons, which are of vital importance in monitoring the seasonal characteristics of climate change. The new algorithm in this study generated a synthetic 30-m high temporal-spatial resolution albedo to solve this problem. In the growing season of GWN, we differentiated between the inversion results from 9 April 2014 and that from adjacent dates to identify daily changes in black-sky albedo in the 15 × 15 km range ( Figure 11). As the new algorithm modulated albedo values, the significant variations in albedo (difference > 0.01) were captured (Figure 11a-l; Figure 11c-p). We calculated the daily mean albedo of the two products and the spatial correlation between two products' variation in the eight days as shown in Table 6. The results indicated that NP not only reflected the daily albedo change such as MODIS with R 2 higher than 0.7, but also described the range of change more finely. Table 6. Statistical results of the mean albedo of the two products, and the correlation between two products' variation in the eight days adjacent to 9 April 2014 at GWN (Figure 11).  Table 6. Statistical results of the mean albedo of the two products, and the correlation between two products' variation in the eight days adjacent to 9 April 2014 at GWN (Figure 11).

Performance Evaluation and Seasonal Deviation Analysis
High temporal-spatial resolution albedo is helpful for understanding the relationship between surface energy balances, vegetation cover dynamics, and ecosystem disturbance. In this study, NP had a lower uncertainty (RMSE of 8%) than CAP (RMSE of 14%), which was approaching the requirement of the World Meteorological Organization (5%) [17]. Therefore, NP had favorable usability and robustness to produce high spatial-temporal resolution albedo. Similar to previous retrieval results using the concurrent approach [37,39], CAP was underestimated at all sites, compared with field measurements ( Table 5). The deviation is mainly due to the underestimation of dormant albedo at non-representative sites. The reason could be that the single-type composition of these sites still contained a variety of soil conditions (such as soil moisture, soil hardness, and planting types) [39,45]. In addition, the fragmented vegetation at these sites is dominated by human disturbance events, because the GWN site is covered with pasture, and both BON and PSU are agricultural research fields. This means that albedo decreases rapidly after crops are planted and increases rapidly after the crops are harvested. For these non-representative areas, the single-type anisotropy information collected by the concurrent approach contains a variety of other vegetation information, which hinders the accurate determination of BRDFs [45]. Therefore, in these complex areas where land cover changes faster or slower than the average rate of land cover in single type, seasonal deviations are likely to occur.
Fortunately, because the new algorithm incorporates the addition of prior NIR AFX [79] and NIR AN ratio [37] in spatial differences, these complicated areas could be further distinguished into combinations of various homogeneous surfaces, so that the seasonal deviations were taken into account in an improved way. The retrieval errors have been limited to each 500-m grid by the use of prior knowledge, instead of accumulating with changes in spatial position. Therefore, the new method has shown favorable retrieval capabilities at all sites, especially for those that are non-representative.

Spatial Heterogeneity
Regarding the spatial heterogeneity of three albedo products, Landsat albedo products were significantly higher than that of MODIS ( Figure 9). Especially, NP had better spatial expression capability ( Figure 10B). However, there was a consistent mismatch of mean albedo among the MODIS albedo and NP in the extended heterogeneous area of GWN (Table 6). According to Cescatti et al. [84], MODIS albedo retrievals for non-forested land underestimate the field measurements of the entire seasonal cycle. This mismatch is likely due to the extreme fragmentation of these landscapes [84].

Quality Control Suggestions and Limitations
Surface albedo is sensitive to atmospheric conditions, which the synthetic albedo products cannot control the quality of albedo at simulated illuminations. Although MODIS quality markers cannot describe the atmospheric conditions under the specific simulated illumination of Landsat, data from multiple samples taken at different times gives a better picture of the day's weather. Additionally, the reliability of these variations depends on MODIS albedo quality markers ( Figure 11C), because the modulation term C SW at MODIS scale is limited by the quality of satellite acquisitions [20]. Therefore, we suggested to use both MODIS quality markers and BRDF inversion markers for the 30-m albedo quality control to obtain the high-quality albedo.
In agricultural areas of GWN, spatial variations in albedo may be caused by crop growth, human disturbance, or small-scale microclimates [40]. As a result, these finer resolution albedo images not only captured details of surface heterogeneity, but also reflected the small-scale changes in albedo caused by variations in surface conditions, which is helpful for monitoring continuous small-scale changes in albedo (e.g., vegetation phenology, fire, snow, and melt). Further study is needed to extend this study to different land cover over a longer time-series, to accurately obtain the seasonal characteristics of heterogeneous ground.

Other Sources of Errors and Potential Improvements
There are potential random and systematic errors in data processing, mainly caused by MODIS/Landsat scale conversions and narrow-to-broad band conversion. Here, we focused on the errors of the MODIS/Landsat scale conversions. First, it retrieves favorable results for Landsat albedo, using NIR AFX and NIR AN ratio as a connection between Landsat and MODIS at different scales, to convert heterogeneous surfaces into appropriate pure pixel values. Although using these two types of prior knowledge can most effectively produce appropriate information for ground objects, outputs may be incorrect for some land surfaces, because the properties of some 30-m surfaces may not fall within the 500-m resolution of information provided by high-quality BRDF. Therefore, further research on these specific land surfaces is needed, to find a more reliable retrieval method that suits them. Second, although the simplified data fusion method can bypass the need to convert multi-band spectral response functions and capture variations in albedo, the modulated results might be unreliable because the 500-m smooth modulation term, C, may not match fragmented Landsat patches. Hence, a more advanced method is required to reasonably disaggregate the 500-m modulation term, C, to a 30-m resolution.
In addition, we note that another potential source of errors in albedo retrieval may be the algorithm itself. With regard to the studies on validating albedo at ground observation stations, the albedo retrieved by the kernel-driven model has generally lower values than those from in situ site measurements [62,75]. Liu et al. [60] pointed out that with increases in SZA, the negative deviation and RMSE of the albedo retrieval results from the kernel-driven model increase. This means that the uncertainties in albedo retrieval might be related to the SZA through the kernel-driven model. In this study, because Landsat's observation zenith angle and relative azimuth angle are both set to 0 to simplify the calculation, the value of the volume scattering kernel is always zero, which would result in erroneous BRF calculations for some SZAs via the kernel-driven model. Therefore, using the kernel-driven model to calculate the AN ratio might produce inaccuracies for some SZAs, causing systematic errors in the retrieval results. Further study is needed to validate the reliability of AN ratio.
Practical solutions can reduce procedural difficulties and errors. In the Landsat 180 × 180 km scene range, using prior knowledge could reduce the number of pure pixels that meet requirements or even lead to retrieval failure. To increase the possibility of pure pixels being identified, Li et al. [46] suggested that multiple adjacent images collected on the same day should be mosaicked to increase the chance of finding these pure pixels. It is also possible to couple VIIRS BRDF with Landsat surface reflectance, and then select for qualified VIIRS pure pixels to improve the quality of albedo retrievals [46]. Additionally, terrain effects can cause distortion of the reflection shape and hemispherical distribution, thereby affecting surface reflectivity [85]. In future research, a three-dimensional kerneldriven model could be used to reduce the impact of terrain effects on albedo retrieval.

Conclusions
This study proposed a new method for estimating Landsat albedo by coupling MODIS BRDF and Landsat-8 surface reflectance. These efforts have been verified at six SURFRAD sites, and the main conclusions are the following:

•
The results demonstrated that the new method had favorable usability and robustness, with less uncertainty (RMSE of 0.015 [8%]), less offset (mean bias of −0.005), and better accuracy (median deviation of −0.011) than CAP. In particular, the RMSE for non-representative sites was significantly improved, mainly by corrections for seasonal deviations. • The new method can limit the seasonal deviations and capture subtle changes in surface albedo of an extended heterogeneous surface. As a result, the new method expands the capacity to retrieve albedo for complex heterogeneous surfaces, because the retrieval errors have been limited to each 500-m grid by prior BRDF knowledge, instead of accumulating with changes in spatial position.

•
The new Landsat albedo product can accurately, finely, and continuously reveal more dynamic surface information. In addition, this simple operability could enable users to continuously and accurately retrieve albedo products with high spatial and temporal resolution in the absence of other auxiliary data (such as topography, land cover types, or disturbance of nature). Therefore, the new method is quite practical and thus very attractive.
Future work will be devoted to the additional verification of Landsat albedo products for different types of surface areas over a longer period, and consider the impact of topographical effects.

Data Availability Statement:
The MODIS BRDF/Albedo data (Collection 6) were downloaded from the LPDAAC of the US NASA and USGS partnership (https://lpdaac.usgs.gov). The Landsat data were obtained from the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) (https://espa.cr.usgs.gov). Thanks to NOAA SUR-FARD for the field albedo measurements (https://gml.noaa.gov/grad/surfrad/index.html). Images of the six ground sites used in this study were acquired from high-resolution satellite data on https://www.google.com/earth.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this article: AFX Anisotropic Flat Index AN ratio Albedo-to-Nadir reflectance ratio AVIRIS