*2.2. Methods*

#### 2.2.1. Current HLS BRDF Normalization

The HLS product (version 1.4) is normalized for BRDF effects using the method proposed by Roy [5]. This method uses fixed BRDF coefficients for each spectral band (i.e., a constant BRDF shape), derived from a large number of pixels in the MODIS 500 m BRDF product (MCD43) that are globally and temporally distributed.

$$
\rho^N(\theta\_s^{out}, 0, 0; \lambda) = \rho(\theta\_s, \theta\_{\upsilon}, \phi; \lambda) \* c(\lambda) \tag{1}
$$

$$x(\lambda) = \frac{k\_0(\lambda) + k\_1(\lambda)F\_1(\theta\_s^{out}, 0, 0) + k\_2(\lambda)F\_2(\theta\_s^{out}, 0, 0)}{k\_0(\lambda) + k\_1(\lambda)F\_1(\theta\_s, \theta\_v, \phi) + k\_2(\lambda)F\_2(\theta\_s, \theta\_v, \phi)} \tag{2}$$

where *θs* is the sun zenith angle, *θv* is the view zenith angle, *φ* is the relative azimuth angle, *F*1 is the volume scattering kernel, based on the Rossthick function derived by Roujean [18], and *F*2 is the geometric kernel, based on the Li-sparse model [19] but considering the reciprocal form given by Lucht [20] and corrected for the Hot-Spot process proposed by Maignan [21]. *θout s* is the sun zenith angle of the normalized data, which is set to constant value depending on the tile central latitude [1]. The BRDF coefficients of the model (*k*0, *k*1, and *k*2) are provided in [1,5].

#### 2.2.2. Proposed BRDF Normalization Method

For our proposed method, we first subset a 9 × 9 km area centered over the albedometer tower locations. All images are then masked for clouds and cloud shadow effects using the product's cloud mask. Next, we run an unsupervised classification over every scene using the Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) developed by Ball, G.H. and Hall, D.J. (1965). This classification is used to unmix the BRDF parameters retrieved from MODIS (at 1 km spatial resolution) to HLS spatial resolution (30 m). We run the ISODATA algorithm for maximum three iterations with the maximum number of clusters setup at 15.

Following the Vermote et al. (2009) notations, the surface reflectance (*ρ*) is written as:

$$\rho(\theta\_{s}, \theta\_{v}, \phi) = k\_{0} \left[ 1 + \frac{k\_{1}}{k\_{0}} F\_{1}(\theta\_{s}, \theta\_{v}, \phi) + \frac{k\_{2}}{k\_{0}} F\_{2}(\theta\_{s}, \theta\_{v}, \phi) \right] \tag{3}$$

where *F*1 and *F*2 are fixed functions of the observation geometry, and *k*0, *k*1, and *k*2 are free parameters. From this notation, we define *V* as *k*1/*k*0 and *R* for *k*2/*k*0. In order to invert the MODIS BRDF

parameters (*V* and *R*) we use the VJB method [12,13]. This method assumes that the difference between the successive observations in time is mainly attributed to directional effects while the variation of *k*0 is supposed small. Additionally, it assumes that *R* and *V* are represented by a linear function of the NDVI.

$$V = V\_0 + V\_1 \text{ NDVI} \tag{4}$$

$$R = R\_0 + R\_1 \text{ NDVI} \tag{5}$$

We unmix the BRDF parameters to the HLS 30 m spatial resolution, by assuming that the surface reflectance of a MODIS pixel can be written as a weighted sum of *n* Landsat classes:

$$\begin{split} k\_0^{\text{Lkm}\_x} + k\_1^{\text{Lkm}\_x} F\_1 & \quad (\theta\_\nu, \theta\_\nu, \mathcal{Q}) + k\_2^{\text{Lkm}\_x} F\_2(\theta\_\nu, \theta\_\nu, \mathcal{Q}) \\ &= A\_x \Big( k\_0^{\text{C}\_1} + k\_1^{\text{C}\_1} F\_1(\theta\_\nu, \theta\_\nu, \mathcal{Q}) + k\_2^{\text{C}\_1} F\_2(\theta\_\nu, \theta\_\nu, \mathcal{Q}) \Big) \\ &= B\_x \Big( k\_0^{\text{C}\_2} + k\_1^{\text{C}\_2} F\_1(\theta\_\nu, \theta\_\nu, \mathcal{Q}) + k\_2^{\text{C}\_2} F\_2(\theta\_\nu, \theta\_\nu, \mathcal{Q}) \Big) + \cdots \\ &= N\_x \Big( k\_0^{\text{C}\_n} + k\_1^{\text{C}\_n} F\_1(\theta\_\nu, \theta\_\nu, \mathcal{Q}) + k\_2^{\text{C}\_n} F\_2(\theta\_\nu, \theta\_\nu, \mathcal{Q}) \Big) \end{split} \tag{6}$$

where *Ax*, *Bx* ... , *Nx* represent proportions of each the *n* classes within the *x* MODIS pixel. BRDF parameters for each class, which are unknowns in Equation (6), can be derived through matrix inversion. We only invert the *k*1 and *k*2 parameters since they describe the directional shape, while *k*0 describes the reflectance magnitude. Considering the HLS surface reflectance and the classification image, we derive *kHLS* 0as shown in Equation (7).

$$k\_0^{HLS} = \rho^{HLS}(\theta\_{\rm s}, \theta\_{\rm v}, \phi) - k\_1^{HLS} F\_1(\theta\_{\rm s}, \theta\_{\rm v}, \phi) - k\_2^{HLS} F\_2(\theta\_{\rm s}, \theta\_{\rm v}, \phi) \tag{7}$$

Figure 2 shows the red band BRDF parameters inverted using this method on an HLS image on June 20th of 2017 over the Yuma desert in Arizona (US).

**Figure 2.** Red band (**a**) R and (**b**) V parameters applying inverting Equation (6) over an HLS image (tile 11SQS) on June 20th of 2017 in Yuma, Arizona (US).

The BRDF inversion is applied to each image. However, this may generate noise in the BRDF inversion caused mainly by clouds. The presence of clouds reduces the number of cloud-free pixels considered in the BRDF inversion, limiting the number of observations available for inversion (note that just images with a minimum of 30% cloud free pixels can be processed to have enough information to invert the model). Additionally, the inaccuracy of the cloud mask [22] may introduce noise in the BRDF unmixing of a given class.

Based on the individual BRDF unmixing of each image in the HLS database (starting from 2013 for Landsat 8, 2015 for Sentinel 2A and 2017 for Sentinel 2B), our goal is to reproduce the BRDF parameters *V*0, *V*1, *R*0, and *R*1 in Equations (4) and (5) at HLS spatial resolution. Next, we apply a linear regression for the *V* and *R* parameters versus the NDVI to derive the *V*0, *V*1, *R*0, and *R*1 parameters at HLS pixel level.

Finally, with the BRDF parameters at HLS resolution, we apply Equation (8) to derive the normalized surface reflectance (*ρ<sup>N</sup>*) at 45 degrees of solar zenith angle and nadir observation (Vermote et al., 2009):

$$\rho^N(45,0,0) = \rho(\theta\_{\delta\prime}\theta\_{\upsilon\prime}\phi) \ast \frac{1 + VF\_1(45,0,0) + RF\_2(45,0,0)}{1 + VF\_1(\theta\_{\delta\prime}\theta\_{\upsilon\prime}\phi) + RF\_2(\theta\_{\delta\prime}\theta\_{\upsilon\prime}\phi)}\tag{8}$$

#### 2.2.3. Temporal Evaluation of Homogeneous Sites

We test the performance of the BRDF normalization by applying the proposed method to two homogeneous sites: forest and desert. To evaluate the correction and assuming that these surfaces remain stationary, we estimate the coefficient of variation (CV), for the reflectance time series. The CV is defined as the ratio of the standard deviation to the mean (Equation (9)).

$$\text{CV} = \frac{stddev(\rho\_{0\prime}\rho\_{1\prime}\dots\rho\_{n})}{average(\rho\_{0\prime}\rho\_{1\prime}\dots\rho\_{n})} \tag{9}$$

where *ρi* is the surface reflectance of day *i*.

#### 2.2.4. Spatial Evaluation of an Equatorial Region

We evaluate the spatial impact of the view zenith angle variation in a dense Amazonian forest close to the equator. To do this, we extract a transect consisting of seven circular homogeneous regions of 800 pixels next to the Urubu river, which cover different values of the VZA. The condition for homogeneity is that the standard deviation is lower than 5% in any band. From these we can estimate the average and standard deviation for each band. The results from using the directional reflectance, the HLS BRDF normalized, and the proposed BRDF normalization are then compared. A *p*-value for each case was calculated which provides a test of the null hypothesis that a coefficient for each regression is equal to zero, i.e., the regressor has no effect on the dependent variable (VZA in our case). A small *p*-value indicates that the null hypothesis can be rejected, meaning that there is a statistically significant relationship between the regressor and dependent variable. We used this statistic index to assess how significant the reflectance dependency with the VZA was.
