2.2.1. VJB Method

The directional effects correction is done using the VJB method, which is based on the Ross Thick Li-Sparse Reciprocal (RTLSR) BRDF model [25,26]. This model expresses the surface reflectance *ρ* as:

$$\rho(\theta\_{\text{v}}, \theta\_{\text{s}}, \phi, t) = k\_{\text{iso}}(t) \left[ 1 + \mathbb{R} F\_{\text{geo}}(\theta\_{\text{v}}, \theta\_{\text{s}}, \phi) + \mathbb{V} F\_{\text{val}}(\theta\_{\text{v}}, \theta\_{\text{s}}, \phi) \right] \tag{2}$$

where, (*Fgeo*, *Fvol*) are the geometric and volumetric scattering components that provide the basic BRDF shapes for characterizing the heterogeneous scattering of the soil-vegetation system; (*kiso*, *kgeo*, *kvol*) are the weight coefficients of the isotropic, geometric and volumetric kernel functions; **V** = *kvol*/*kiso*, **R** = *kgeo*/*kiso*.

It assumes that the difference between consecutive observations is due principally to directional effects, with small variations of the isotropic weight coefficient kiso.

$$\frac{\rho(t\_i)}{1 + \mathsf{V} F\_{\text{vol}}^i + \mathsf{R} F\_{\text{geo}}^i} \approx \frac{\rho(t\_{i+1})}{1 + \mathsf{V} F\_{\text{vol}}^{i+1} + \mathsf{R} F\_{\text{geo}}^{i+1}}.\tag{3}$$

The values of V and R are the unknowns and can be derived using N observations from a certain pixel by minimizing the merit function:

$$M = \sum\_{i=1}^{N-1} \frac{\left(\rho\_{i+1} \left[1 + VF\_1^i + RF\_2^i\right] - \rho\_i \left[1 + VF\_1^{i+1} + RF\_2^{i+1}\right]\right)^2}{day^{i+1} - day^i + 1}.\tag{4}$$

The minimization is done by the classic derivation:

$$\frac{dM}{dV} = \frac{dM}{dR} = 0.\tag{5}$$

This leads to the system of equations:

$$
\begin{pmatrix}
\begin{array}{cc}
\sum^{N-1}\Delta^{i}\rho F\_{1}\Delta^{i}\rho F\_{1} & \sum^{N-1}\Delta^{i}\rho F\_{1}\Delta^{i}\rho F\_{2} \\
\sum^{N-1}\Delta^{i}\rho F\_{1}\Delta^{i}\rho F\_{2} & \sum^{N-1}\Delta^{i}\rho F\_{2}\Delta^{i}\rho F\_{2}
\end{array} \\
\begin{array}{c}
\sum^{N-1}\Delta^{i}\rho F\_{1}\Delta^{i}\rho F\_{2}\Delta^{i}\rho F\_{2}
\end{array}
\end{pmatrix} \otimes \begin{pmatrix} V \\ R \end{pmatrix} = \begin{pmatrix}
\end{pmatrix} \tag{6}
$$

where:

$$
\Delta^i d = d a y\_{i+1} - d a y\_i + 1$$

$$
\Delta^i \rho = (\rho\_{i+1} - \rho\_i) / \sqrt{\Delta^i d}$$

$$
\Delta^i \rho F\_{1,2} = \left(\rho\_{i+1} F\_{1,2}^i - \rho\_{i+1} F\_{1,2}^{i+1}\right) / \sqrt{\Delta^i d} \,. \tag{7}$$

Through the inversion of Equation (6), we can now obtain a V and R parameter for every pixel. To perform the instantaneous directional correction for every observation, V and R parameters are needed for every date, so the inversion of the V and R parameters is done for five different NDVI populations. A linear regression is performed for V and R as a function of the NDVI (Equation (8)). This retrieves a slope and intercept for every pixel which allow the computation of V and R parameters for a certain date given the surface's NDVI value.

$$V = V\_0 + V\_1 \* N DVI$$

$$R = R\_0 + R\_1 \* N DVI.\tag{8}$$

 These V and R can now be used to calculate the normalized surface reflectance (*ρ<sup>N</sup>*) at *θs* = 45◦, and nadir observation using Equation (9) [18]:

$$\rho^N(45,0,0) = \rho(\theta\_s, \theta\_{v\prime}, \phi) \ast \frac{1 + VF\_1(45,0,0) + RF\_2(45,0,0)}{1 + VF\_1(\theta\_{s\prime}, \theta\_{v\prime}, \phi) + RF\_2(\theta\_{s\prime}, \theta\_{v\prime}, \phi)}\tag{9}$$

The NDVI can now be obtained using the normalized reflectances:

$$NDVI^N = \frac{\rho\_{NIR}^N - \rho\_{rad}^N}{\rho\_{NIR}^N + \rho\_{rad}^N}.\tag{10}$$

## 2.2.2. BRDF Parameters Relationship

To minimize the propagation of the atmospheric errors into the BRDF correction, we analyze the V and R parameters for bands 1 (red) and 2 (NIR). The goal is to find a physical relationship between the BRDF parameters of both bands to avoid using the noisy red AVHRR band. Therefore, we can derive the BRDF parameters of the NIR band and then estimate the red band based on these parameters. To build this physical relationship we use MODIS data since we need data with the least error possible.

To do this, we extract the band 1, band 2 and NDVI time series of one pixel. We then sort them in ascending values of NDVI. We divide these values into five groups, using the 20th, 40th, 60th and 80th percentiles as group edge values. Each of these groups is defined as an NDVI population with its own average. We then extract the V and R values of each one using their band 1, band 2 and NDVI. This means that for every pixel and every band, we obtain 5 different V and R values. When using all the BELMANIP2 sites we obtain a total of 445 × 5 = 2225 points for each band. Finally, we do simple regressions between the obtained parameters for the different bands and with the NDVI itself.
