*3.3. Surface Albedo Validation*

We also evaluate the BRDF normalization by using the BRDF parameters to derive the surface albedo, which can be validated with field measurements. Figure 8 shows the validation of the surface albedo over (a) the SURFRAD network, (b) the OzFlux network, and (c) both networks. Each marker in the plot represents an individual Landsat 8 or Sentinel 2 overpass. OzFlux validation shows lower errors (0.015–0.016) than SURFRAD (0.020). Sentinel 2 albedo validation in SURFRAD shows an overestimation for albedos higher than 0.2. However, this is observed in a total of 5 dates in Desert Rock and may be caused by an inaccurate atmospheric correction during those dates. Nevertheless, the overestimation is included within the RMSE of the sites. Combining both networks, the surface albedo validation shows a RMSE of 0.019 for Landsat 8 and 0.018 for Sentinel 2. Figure 8d shows the result of applying the narrow to broadband equation to the directional reflectance. In this case, the errors increase to 0.028 for Landsat 8 and 0.030 for Sentinel 2. Note that an equivalent evaluation cannot be included for the current HLS BRDF normalization method because in Roy [5] it is explicitly stated that "this approach is not applicable for generation of Landsat surface albedo, which requires a full characterization of the surface BRDF".

**Figure 8.** Broadband blue sky surface albedo validation of all the (**a**) SURFRAD, (**b**) OzFlux sites, and (**c**) combining both sites considered from 2013 to 2017. (**d**) The broadband directional surface reflectance comparison with surface albedo measurements.

#### **4. Discussion and Conclusions**

Quality of atmospheric and BRDF correction is essential for any application of the HLS dataset. In this work, we present a new method to improve the BRDF normalization of the HLS data. It is based on the Franch [11] method, which was originally designed and evaluated to derive the surface albedo of Landsat TM and ETM+ data. In this manuscript we improve it, apply it to HLS data, explore its feasibility to correct the BRDF effects of the dataset, and validate it over two homogeneous sites and several land cover types of the SURFRAD and OzFlux networks.

The proposed improvement is based on reproducing the BRDF parameters' NDVI dependency [12,13] at moderate spatial resolution (HLS). In the operational context, this means that we create an HLS pixel-based database of the BRDF parameters (V0, V1, R0, and R1) that is regularly updated to account for any land cover change. In this way, the effect of outliers and poor-quality pixels on the BRDF

inversion is minimized, resulting in a more stable and robust BRDF model. We acknowledge that the assumption of a BRDF model being a function of the NDVI has limitations. For example, on sparse forests where the NDVI is not a good descriptor of canopy structure [24]. However, this simple model might be used for a rough correction of BRDF effects in reflectance time series. Although a full inversion of the BRDF model will give better results, some applications, such as real time processing, may want to trade accuracy for simplicity [25].

When applying the proposed method, the results show a decrease in the surface reflectance timeseries CV of 4% and a decrease of the correlation coefficient with the SZA for the forest site, and little to no dependency on the desert site (which is barely affected by angular effects). In contrast, the current HLS BRDF normalization algorithm under-corrects the BRDF effects on the forests site, and increases the CV on the desert site. It is known that vegetation in the Amazon has a phenological cycle as a result of small variations in temperature and precipitation throughout the year [26]. However, directional effects on surface reflectance can be larger than the differences in reflectance between floristically different vegetation types [27,28].

The evaluation of the spatial variability of the view zenith angle showed a clear dependency in the backscatter direction, leading to higher directional reflectance for larger view angles. The current BRDF method reduces this dependency, but still shows an under-correction of the signal in the dense forest area analyzed. This under correction was already reported by Roy [5]. Finally, the proposed method shows no dependency with the view zenith angle.

The evaluation against surface albedo field measurements show an RMSE error of 0.019 for Landsat 8 and 0.018 for Sentinel-2. The OzFlux network error is similar to Franch [11] (RMSE of 0.016), while SURFRAD sites show slightly larger errors. Note that Franch [11] evaluated the method over five SURFRAD sites from 2003 to 2006. The higher error in the SURFRAD network compared to OzFlux might be caused by higher errors in the atmospheric correction algorithm since the atmosphere over Australia exhibit lower aerosol content than over the US. Compared to previous studies that provide albedo validation, the proposed method shows lower errors than Shuai [6] (RMSE 0.024). These results are further evidence of the good performance of the surface reflectance product of Landsat 8 and Sentinel 2.

**Author Contributions:** Conceptualization and methodology, B.F. and E.V.; data curation: S.S., J.J. and A.S.-A.; supervision: J.-C.R. and J.M.; validation B.F.; writing: B.F., J.L.V.-N. and A.S.-A.

**Funding:** This work was supported by the NASA gran<sup>t</sup> "Support for the HLS (Harmonized Landsat-Sentinel-2) Project" (no. NNX16AN88G).

**Acknowledgments:** This work was supported by the NASA gran<sup>t</sup> "Support for the HLS (Harmonized Landsat-Sentinel-2) Project" (no. NNX16AN88G).

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