**1. Introduction**

The Harmonized Landsat/Sentinel-2 (HLS) project [1] provides a surface reflectance product that combines observations from USGS/NASA's Landsat 8 and ESA's Sentinel-2 satellites at moderate spatial resolution (30 m). The main goal is to provide a unique dataset based on both satellites' data to improve the revisit time to 3–5 days depending on the latitude [1]. Along with a common

atmospheric correction algorithm [2], geometric resampling to 30 m spatial resolution and geographic registration [1], the product is also corrected for BRDF effects and band pass adjustment.

The purpose of the band pass adjustment is to correct discrepancies between data coming from different sensors due differences their Relative Spectral Response (RSR) functions. In the case of the near infrared (NIR) spectral band, the spectral difference between Sentinel-2 and Landsat-8 is ~0.03%, due to the almost identical spectral functions of the analogous bands. In the red band, however, these differences are of ~3% [3] and need to be corrected. To correct the differences in the red band, Villaescusa-Nadal [3] employed a Spectral Band Adjustment Factor (SBAF) exponential model dependent on the NDVI which reduced the discrepancy between the sensors by ~55%.

Regarding the BRDF correction, Landsat and Sentinel-2 sampling characteristics provide nearly constant observation geometry and low illumination variation within each scene. However, when a surface reflectance time series combining measurements from both sensors is created, variations due to differences in view geometry between them arise. In extreme cases, the differences in the view zenith angle for a ground target can reach 20◦ from adjacent orbits only a few days apart. Additionally, variations in the seasonal illumination also impact the surface reflectance value. Gao [4] concluded that for Landsat-like narrow swath sensors, the major BRDF effect arises from the day of year effect, and can cause variations of 0.04–0.06 reflectance compared to mid-summer observations. Such angular effects can be corrected using a Bidirectional Reflectance Distribution Function (BRDF) model. However, the narrow angular sampling of moderate resolution sensors such as Landsat and Sentinel 2 complicates the BRDF parameters retrieval.

The official HLS product (version 1.4) provides a BRDF normalized product based on the BRDF normalization proposed by Roy [5]. In this product, the view angle is set to nadir and the solar zenith angle is fixed through time but varies for each tile based on the latitude. However, this correction is not optimal because the considered coefficients are constant and do not depend on the surface type or condition. Previous studies have worked towards BRDF estimation at moderate resolution for correcting the BRDF effects or for generating an albedo product of medium spatial resolution sensors. Shuai [6] used the MODIS Bidirectional Reflectance Distribution Function (BRDF) MDC43 product at 500 m [7] to generate a Landsat surface albedo product that later was implemented back to the 80s to the "pre-MODIS era" [8]. Their method achieves a Root Mean Square Error (RMSE) generally lower than 0.03 (12%) when compared with in situ data. However, this method does not provide a BRDF product at Landsat scale. Flood [9] used constant BRDF parameters to adjust surface reflectance from Landsat (TM and ETM+) and SPOT 5 images to a standard angular configuration. Gao [4] proposed a method to correct BRDF effects on Landsat and the advanced wide field sensor (AWiFS) based on detecting homogeneous areas with USDA cropland data layer (CDL) map to build a BRDF Look-Up Map (LUM) based on MODIS BRDF data. Van Doninck [10] assessed the directional effects on Landsat TM/ETM+ imagery over Amazonian forests and evaluated five different normalization techniques, including MODIS-based and Landsat-based BRDF. They concluded that the two MODIS-based methods analyzed only removed half of the angular effect in the infrared bands while the best results were obtained parametrizing the BRDF model using pairs of Landsat images. Franch [11] presented an alternative algorithm to calculate the Landsat (TM and ETM+) surface albedo based on the BRDF parameters estimated from the MODerate Resolution Imaging Spectroradiometer (MODIS) Climate Modeling Grid (CMG) surface reflectance product (M{O,Y}D09) using the VJB method [12,13]. The validation of the results against field measurements showed an RMSE of 0.015 (7%). Franch [14] tested this method to correct the BRDF effects of the HLS data over the Amazon forest in Peru. Their results show a good performance of the model when applied to HLS data with an error of 0.01 in the surface albedo retrieval.

In this paper, we describe the adaptation of the method of Franch [11,14] for HLS BRDF normalization, apply the proposed method to the HLS product from 2013 to 2017 and validate its applicability to Landsat 8 and Sentinel 2 data, test it on three homogeneous sites to assess its capability to remove BRDF artifacts, and finally evaluate it through the comparison to ground albedo measurements over SURFRAD and OzFlux networks.

## **2. Materials and Methods**
