# Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints

Jonathan León-Tavares, Jose Gómez-Dans, Jean-Louis Roujean, Véronique Bruniquel

► To cite this version:

Jonathan León-Tavares, Jose Gómez-Dans, Jean-Louis Roujean, Véronique Bruniquel. Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints. Remote Sensing of Environment, 2024, 302, pp.113967. 10.1016/j.rse.2023.113967. hal-04599891

HAL Id: hal-04599891

https://hal.science/hal-04599891

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints

Jonathan León-Tavares, Jose Gómez-Dans, Jean-Louis Roujean, Véronique Bruniquel

► To cite this version:

Jonathan León-Tavares, Jose Gómez-Dans, Jean-Louis Roujean, Véronique Bruniquel. Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints. Remote Sensing of Environment, 2024, 302, pp.113967. 10.1016/j.rse.2023.113967. hal-04599891

HAL Id: hal-04599891

https://hal.science/hal-04599891v1

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints

Jonathan León-Tavares<sup>a, *,</sup> , Jose Gómez-Dans<sup>b, c, f</sup> , Jean-Louis Roujean<sup>d</sup> , Véronique Bruniquel<sup>e</sup>

<figure><img src="image_3.png" /></figure>

## ABSTRACT

Sentinel-3 supports governmental policies and European programmes such as the Copernicus Global Land Service with daily monitoring of the land surface by measuring the Earth's surface with its multi-band wide-swath visible and near-infrared radiometers. To fully exploit Sentinel-3 data sets, performing a correction of directional effects caused by the land surface reflectance anisotropy properties is mandatory. This correction requires first an estimation of the bidirectional reflectance distribution function (BRDF) to be properly achieved. The quality and robustness of global estimates of key land surface variables derived from space-borne sensors depend clearly on an accurate assessment of the spectral BRDF.

Here, we present an algorithm aiming at harnessing the Sentinel-3's wide imaging swath to retrieve the land surface BRDF, approximating it via a kernel-driven semi-empirical model which can then be used to normalise the observations to a common Sun-sensor geometry. Our algorithm, named ReBeLS (Regularised BRDF inversion for Land Surface), uses a temporally regularised BRDF inversion approach and, together with prior knowledge of the land surface BRDF obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors, assimilates Sentinel-3 surface reflectance observations. We focus on the results derived from observations provided by the Sentinel-3 Ocean and Land Colour Instrument (OLCI) sensors; however, the algorithm can also be applied to reflectance observations acquired by the Sea and Land Surface Temperature Radiometer (SLSTR). The retrieved BRDF model parameters reproduce within uncertainties the directional signature as seen by the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor. The generated Sentinel-3 normalised surface reflectance products are a significant advance compared with time series of directional surface reflectance and maximum value composite images, showing a statistically significant (p<0.05) correlation (r > 0.9) with VIIRS nadir BRDF-adjusted surface reflectance. Furthermore, the ReBeLS Sentinel-3 BRDF products offer enough flexibility to allow users to compute spectral surface albedos or to normalise surface reflectance to any Sun-sensor configuration.

## 1. Introduction

The Sentinel-3 (S3) constellation sensors' large imaging swaths – Ocean and Land Colour Instrument (OLCI): 1270 km<fnref n="1" />; Sea and Land Surface Temperature Radiometer (SLSTR): 1400 km [nadir] and 740 km [oblique]<fnref n="2" /> – secure a daily coverage of the Earth at the expense of a different geometry of observations between satellites overpasses. It is well known that the radiance of the Earth's land surface is highly anisotropic (Breon and Maignan, 2017). Consequently, measured radiances from Sentinel-3 are impacted by the sensor acquisition geometry and the Sun's spectral irradiance direction. It follows that those radiances issued from the different parts of the Sentinel-3 wide imaging swaths cannot be directly compared without performing an anisotropy

effects correction.

To fully exploit the high revisit period provided by the Sentinel-3 constellation, knowledge of the surface's bidirectional reflectance distribution function (BRDF) is required. The BRDF describes the radiation scattering in one direction given the direction of illumination. By approximating the bidirectional reflectance factor (BRF) (Nicodemus et al., 1977; Schaepman-Strub et al., 2006) derived from the observed radiances of an arbitrary surface with a BRDF model, the directional effects induced by the land surface anisotropy can be modelled and its effects removed from the radiometric signal. The latter procedure is commonly known as BRDF correction, which leads to an angular normalisation or BRDF-adjustment. This can be achieved by finding a set of model parameters (BRDF descriptors) obtained via a BRDF model

<footnote><sup>a</sup> Vlaamse Instelling voor Technologisch Onderzoek (VITO), Boeretang 200, 2400 Mol, Belgium</footnote>

<footnote><sup>b</sup> Department of Geography, King's College London, Strand, London WC2B 4BG, UK</footnote>

<footnote><sup>c</sup> NERC National Centre for Earth Observation (NCEO), UK</footnote>

<footnote><sup>d</sup> CESBIO—Centre d'Etudes Spatiales de la BIOsphère, CESBIO UMR 5126, 31401 Toulouse, France</footnote>

<footnote><sup>e</sup> ACRI-ST, 260 Route du Pin Montard, 06904 Sophia Antipolis, France</footnote>

<footnote><sup>f</sup> Leverhulme Centre for Wildfires, Environment and Society, United Kingdom</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: jonathan.leontavares@vito.be (J. León-Tavares).</footnote>

<footnote n="1">¹ https://sentinel.esa.int/web/sentinel/user-guides/sentinel-3-olci/coverage</footnote>

<footnote n="2">² https://sentinel.esa.int/web/sentinel/user-guides/Sentinel-3-slsr/coverage</footnote>
<figure><img src="image_2.png" /></figure>

# Retrieving land surface reflectance anisotropy with Sentinel-3 observations and prior BRDF model constraints

Jonathan León-Tavares<sup>a, *,</sup> , Jose Gómez-Dans<sup>b, c, f</sup> , Jean-Louis Roujean<sup>d</sup> , Véronique Bruniquel<sup>e</sup>

<figure><img src="image_3.png" /></figure>

## ABSTRACT

Sentinel-3 supports governmental policies and European programmes such as the Copernicus Global Land Service with daily monitoring of the land surface by measuring the Earth's surface with its multi-band wide-swath visible and near-infrared radiometers. To fully exploit Sentinel-3 data sets, performing a correction of directional effects caused by the land surface reflectance anisotropy properties is mandatory. This correction requires first an estimation of the bidirectional reflectance distribution function (BRDF) to be properly achieved. The quality and robustness of global estimates of key land surface variables derived from space-borne sensors depend clearly on an accurate assessment of the spectral BRDF.

Here, we present an algorithm aiming at harnessing the Sentinel-3's wide imaging swath to retrieve the land surface BRDF, approximating it via a kernel-driven semi-empirical model which can then be used to normalise the observations to a common Sun-sensor geometry. Our algorithm, named ReBeLS (Regularised BRDF inversion for Land Surface), uses a temporally regularised BRDF inversion approach and, together with prior knowledge of the land surface BRDF obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors, assimilates Sentinel-3 surface reflectance observations. We focus on the results derived from observations provided by the Sentinel-3 Ocean and Land Colour Instrument (OLCI) sensors; however, the algorithm can also be applied to reflectance observations acquired by the Sea and Land Surface Temperature Radiometer (SLSTR). The retrieved BRDF model parameters reproduce within uncertainties the directional signature as seen by the Visible Infrared Imaging Radiometer Suite (VIIRS) sensor. The generated Sentinel-3 normalised surface reflectance products are a significant advance compared with time series of directional surface reflectance and maximum value composite images, showing a statistically significant (p<0.05) correlation (r > 0.9) with VIIRS nadir BRDF-adjusted surface reflectance. Furthermore, the ReBeLS Sentinel-3 BRDF products offer enough flexibility to allow users to compute spectral surface albedos or to normalise surface reflectance to any Sun-sensor configuration.

## 1. Introduction

The Sentinel-3 (S3) constellation sensors' large imaging swaths – Ocean and Land Colour Instrument (OLCI): 1270 km<fnref n="1" />; Sea and Land Surface Temperature Radiometer (SLSTR): 1400 km [nadir] and 740 km [oblique]<fnref n="2" /> – secure a daily coverage of the Earth at the expense of a different geometry of observations between satellites overpasses. It is well known that the radiance of the Earth's land surface is highly anisotropic (Breon and Maignan, 2017). Consequently, measured radiances from Sentinel-3 are impacted by the sensor acquisition geometry and the Sun's spectral irradiance direction. It follows that those radiances issued from the different parts of the Sentinel-3 wide imaging swaths cannot be directly compared without performing an anisotropy

effects correction.

To fully exploit the high revisit period provided by the Sentinel-3 constellation, knowledge of the surface's bidirectional reflectance distribution function (BRDF) is required. The BRDF describes the radiation scattering in one direction given the direction of illumination. By approximating the bidirectional reflectance factor (BRF) (Nicodemus et al., 1977; Schaepman-Strub et al., 2006) derived from the observed radiances of an arbitrary surface with a BRDF model, the directional effects induced by the land surface anisotropy can be modelled and its effects removed from the radiometric signal. The latter procedure is commonly known as BRDF correction, which leads to an angular normalisation or BRDF-adjustment. This can be achieved by finding a set of model parameters (BRDF descriptors) obtained via a BRDF model

<footnote><sup>a</sup> Vlaamse Instelling voor Technologisch Onderzoek (VITO), Boeretang 200, 2400 Mol, Belgium</footnote>

<footnote><sup>b</sup> Department of Geography, King's College London, Strand, London WC2B 4BG, UK</footnote>

<footnote><sup>c</sup> NERC National Centre for Earth Observation (NCEO), UK</footnote>

<footnote><sup>d</sup> CESBIO—Centre d'Etudes Spatiales de la BIOsphère, CESBIO UMR 5126, 31401 Toulouse, France</footnote>

<footnote><sup>e</sup> ACRI-ST, 260 Route du Pin Montard, 06904 Sophia Antipolis, France</footnote>

<footnote><sup>f</sup> Leverhulme Centre for Wildfires, Environment and Society, United Kingdom</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: jonathan.leontavares@vito.be (J. León-Tavares).</footnote>

<footnote n="1">¹ https://sentinel.esa.int/web/sentinel/user-guides/sentinel-3-olci/coverage</footnote>

<footnote n="2">² https://sentinel.esa.int/web/sentinel/user-guides/Sentinel-3-slsr/coverage</footnote>
inversion (Roujean, 2017).

Because of their good performance and reasonable computational costs, the category of semi-empirical kernel-driven linear BRDF models has a long history in operational chains of production of BRDF and surface albedo products (Carrer et al., 2021; Carrer et al., 2018; Geiger et al., 2008; Schaaf et al., 2002). Inverting kernel-driven BRDF model is often an ill-posed problem (Combal et al., 2003; Kimes et al., 2000) due to typically insufficient angular sampling and the directional distribution of the available sampling (Gao et al., 2002; Lucht and Lewis, 2000). To improve the ill-posed nature of the problem, additional constraints are often deployed, such as the assumption of surface reflectance stability over a relatively short period of time (Schaaf et al., 2002), or its slowly varying evolution (Quaife and Lewis, 2010).

The Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43 product (Schaaf et al., 2002) has received wide usage from the community. The product accumulates all the cloud-free observations in a 16-day time frame to form a data set serving to invert the BRDF model and retrieve the descriptor parameters. For MODIS Collection 6 (Wang et al., 2018), a sliding exponential temporal weighting is used to provide a daily estimate of the BRDF parameters. This approach is supplemented by a backup “magnitude inversion” algorithm (Strugnell et al., 2001) when there are not enough observations within the 16-day window, or outliers due to cloud, cloud shadow or otherwise leading to a poor retrieval. A similar approach is used for the VNP43 product derived from the Visible Infrared Imaging Radiometer Suite (VIIRS) sensors (Liu et al., 2017). According to Wang et al. (2018), over the European continent, in the MCD43 C6 product, only 35% of the retrievals are produced with the full algorithm, whereas 38% and 27% of the retrievals being respectively backup and “no retrieval”. The backup algorithm can provide reliable results, but caution must be exercised. While both MCD43 and VNP43 have a quality assurance (QA) layer, neither of the two products provides a per-pixel uncertainty estimates when such information in derived products has been recently highlighted in the recent ESA Earth observation calibration and validation strategies (Niro et al., 2021). The use of a 16-day window, based on Privette et al. (1997) is inherent to the MODIS revisiting strategy and results in large gaps within a European context.

The Copernicus Global Land Service, Copernicus Climate Change Service, and study 2 of the Copernicus 4 Core Product Algorithm Studies (COPA) project, require the development of a state-of-the-art algorithm (Gastellu-Etchegorry et al., 2022) to retrieve the land surface BRDF from Sentinel-3 constellation's BRF measurements. The goal was to provide the long-term continuity of the BRDF-corrected NDVI (León-Tavares et al., 2021) and surface albedo (Carrer et al., 2021) products as a follow-on of the PROBA-V mission (Dierckx et al., 2014), decommissioned in 2021, with Sentinel-3 (Sánchez-Zapero et al., 2023). To satisfy users requirements of constructing a dynamic database of global BRDF products, the ReBeLS (Regularised BRDF inversion for Land Surface) algorithm has been developed. This approach builds on previous BRDF and albedo products, such as MCD43 or GlobAlbedo, and exploits both the simplicity and the suitability of linear BRDF kernel models. ReBeLS uses a prior estimate of land surface reflectance anisotropy to constrain the parameter inversion, as in GlobAlbedo (Muller et al., 2013), and imposes an additional constraint on the smooth temporal evolution of the BRDF parameters (Quaife and Lewis, 2010). Phrasing the problem in a Bayesian context allows for an effective and traceable treatment of per-pixel uncertainties. The methodology implemented in ReBeLS does not require a backup algorithm and provides an estimate of the land surface BRDF even in periods of persistent cloud coverage, the quality of which can be quantified by its associated uncertainty. Unlike the MCD43 and VNP43 methodologies, the BRDF inversion problem in our approach is stated so that no fixed accumulation period is needed.

Phrasing the inversion problem in Bayesian terms provides an elegant way of coherently adding constraints expressed as priors and propagating uncertainties. One such way of adding extra information is to add an a priori estimate of the model parameters derived from

independent observations, or from a climatology (Muller et al., 2013). An additional useful constrain is the already mentioned smooth variation of the model parameters over time (Quaife and Lewis, 2010; Wang et al., 2018).

This article is structured as follows: Section 2 describes the datasets used to feed and evaluate the algorithm. Section 3 describes the ReBeLS methodology. It is important to notice that the ReBeLS BRDF product is a set of spectral BRDF descriptors that serve as a weight to the kernels of the semi-empirical kernel-driven BRDF model. As such, they cannot be compared directly with in-situ measurements nor to BRDF descriptors retrievals obtained with another sensor. However, the ReBeLS BRDF descriptors can be evaluated by their ability to reproduce the land surface anisotropy observed by another sensor. The results of the latter validation approach are described in Section 4. Finally, our conclusions and the algorithm's outlook are presented in section 5. Hereafter, we indistinctly refer to the BRF as surface reflectance throughout the manuscript, and all mentions of MODIS and VIIRS products concern version V006 and V001, respectively.

## 2. Data

### 2.1. Sentinel-3 observations

The input dataset to the BRDF retrieval methodology presented in this article is the atmospherically corrected surface reflectance and acquisition geometries provided by the Copernicus Global Land Service. Although this Sentinel-3 derived land surface reflectance dataset is not yet in the public domain, its preprocessing stages are fully described in its ATBD.<fnref n="3" /> To summarise the latter, the Sentinel-3 (A and B) OLCI and SLSTR (nadir) L1b data products are coregistered and resampled to a common spatial resolution using the open-source Level 1C SYN tool.<fnref n="4" /> Pixel identification and reprojection processing are performed on the resulting product, referred to as the L1C SYN product.

Further, an atmospheric correction based on SMAC (Rahman and Dedieu, 1994) – and auxiliary meteorological data from the MERRA-2 database<fnref n="5" /> – is applied to the OLCI and SLSTR L1C SYN TOA radiances to obtain Top of Canopy (TOC) surface reflectances. Full details on the atmospheric correction can be found in the Copernicus Global Land Service ATBD of the atmospheric correction applied to Sentinel-3 data.

The Sentinel-3 TOC surface reflectance products are defined on a global 333 m WGS 84 grid divided into 10° × 10° tiles that follow the PROBA-V grid tiling system. These tiles are distributed as NetCDF files. TOC reflectances and uncertainties (associated with atmospheric correction procedure) are provided in each file for 20 Sentinel-3 spectral bands (15 OLCI bands, 5 SLSTR bands). Moreover, four quality information flags and geometry layers for OLCI and SLSTR-nadir are also included in the NetCDF files. For full details on these datasets, we refer the reader to the aforementioned ATBD.

### 2.2. BRDF prior database (MCD43P)

Numerical inversion of kernel-driven semi-empirical models is typically used to approximate the land surface BRDF (for a review, see Roujean, 2017). The linear least-squares method is commonly used to solve the BRDF inversion problem. Sometimes, the solution may prove unstable and physically unacceptable (e.g. negative values associated with BRDF descriptors and BRDF-adjusted surface reflectances). This occurs when the BRDF inversion is undetermined due to insufficient angular sampling, the presence of outliers or other sources of uncertainty in the surface reflectance observations, which leads to an ill-posed

<footnote n="3">³ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_S3-AC-V1_11_30.pdf</footnote>

<footnote n="4">⁴ https://github.com/bcdev/11c-syn-tool</footnote>

<footnote n="5">⁵ https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/</footnote>
inversion (Roujean, 2017).

Because of their good performance and reasonable computational costs, the category of semi-empirical kernel-driven linear BRDF models has a long history in operational chains of production of BRDF and surface albedo products (Carrer et al., 2021; Carrer et al., 2018; Geiger et al., 2008; Schaaf et al., 2002). Inverting kernel-driven BRDF model is often an ill-posed problem (Combal et al., 2003; Kimes et al., 2000) due to typically insufficient angular sampling and the directional distribution of the available sampling (Gao et al., 2002; Lucht and Lewis, 2000). To improve the ill-posed nature of the problem, additional constraints are often deployed, such as the assumption of surface reflectance stability over a relatively short period of time (Schaaf et al., 2002), or its slowly varying evolution (Quaife and Lewis, 2010).

The Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43 product (Schaaf et al., 2002) has received wide usage from the community. The product accumulates all the cloud-free observations in a 16-day time frame to form a data set serving to invert the BRDF model and retrieve the descriptor parameters. For MODIS Collection 6 (Wang et al., 2018), a sliding exponential temporal weighting is used to provide a daily estimate of the BRDF parameters. This approach is supplemented by a backup “magnitude inversion” algorithm (Strugnell et al., 2001) when there are not enough observations within the 16-day window, or outliers due to cloud, cloud shadow or otherwise leading to a poor retrieval. A similar approach is used for the VNP43 product derived from the Visible Infrared Imaging Radiometer Suite (VIIRS) sensors (Liu et al., 2017). According to Wang et al. (2018), over the European continent, in the MCD43 C6 product, only 35% of the retrievals are produced with the full algorithm, whereas 38% and 27% of the retrievals being respectively backup and “no retrieval”. The backup algorithm can provide reliable results, but caution must be exercised. While both MCD43 and VNP43 have a quality assurance (QA) layer, neither of the two products provides a per-pixel uncertainty estimates when such information in derived products has been recently highlighted in the recent ESA Earth observation calibration and validation strategies (Niro et al., 2021). The use of a 16-day window, based on Privette et al. (1997) is inherent to the MODIS revisiting strategy and results in large gaps within a European context.

The Copernicus Global Land Service, Copernicus Climate Change Service, and study 2 of the Copernicus 4 Core Product Algorithm Studies (COPA) project, require the development of a state-of-the-art algorithm (Gastellu-Etchegorry et al., 2022) to retrieve the land surface BRDF from Sentinel-3 constellation's BRF measurements. The goal was to provide the long-term continuity of the BRDF-corrected NDVI (León-Tavares et al., 2021) and surface albedo (Carrer et al., 2021) products as a follow-on of the PROBA-V mission (Dierckx et al., 2014), decommissioned in 2021, with Sentinel-3 (Sánchez-Zapero et al., 2023). To satisfy users requirements of constructing a dynamic database of global BRDF products, the ReBeLS (Regularised BRDF inversion for Land Surface) algorithm has been developed. This approach builds on previous BRDF and albedo products, such as MCD43 or GlobAlbedo, and exploits both the simplicity and the suitability of linear BRDF kernel models. ReBeLS uses a prior estimate of land surface reflectance anisotropy to constrain the parameter inversion, as in GlobAlbedo (Muller et al., 2013), and imposes an additional constraint on the smooth temporal evolution of the BRDF parameters (Quaife and Lewis, 2010). Phrasing the problem in a Bayesian context allows for an effective and traceable treatment of per-pixel uncertainties. The methodology implemented in ReBeLS does not require a backup algorithm and provides an estimate of the land surface BRDF even in periods of persistent cloud coverage, the quality of which can be quantified by its associated uncertainty. Unlike the MCD43 and VNP43 methodologies, the BRDF inversion problem in our approach is stated so that no fixed accumulation period is needed.

Phrasing the inversion problem in Bayesian terms provides an elegant way of coherently adding constraints expressed as priors and propagating uncertainties. One such way of adding extra information is to add an a priori estimate of the model parameters derived from

independent observations, or from a climatology (Muller et al., 2013). An additional useful constrain is the already mentioned smooth variation of the model parameters over time (Quaife and Lewis, 2010; Wang et al., 2018).

This article is structured as follows: Section 2 describes the datasets used to feed and evaluate the algorithm. Section 3 describes the ReBeLS methodology. It is important to notice that the ReBeLS BRDF product is a set of spectral BRDF descriptors that serve as a weight to the kernels of the semi-empirical kernel-driven BRDF model. As such, they cannot be compared directly with in-situ measurements nor to BRDF descriptors retrievals obtained with another sensor. However, the ReBeLS BRDF descriptors can be evaluated by their ability to reproduce the land surface anisotropy observed by another sensor. The results of the latter validation approach are described in Section 4. Finally, our conclusions and the algorithm's outlook are presented in section 5. Hereafter, we indistinctly refer to the BRF as surface reflectance throughout the manuscript, and all mentions of MODIS and VIIRS products concern version V006 and V001, respectively.

## 2. Data

### 2.1. Sentinel-3 observations

The input dataset to the BRDF retrieval methodology presented in this article is the atmospherically corrected surface reflectance and acquisition geometries provided by the Copernicus Global Land Service. Although this Sentinel-3 derived land surface reflectance dataset is not yet in the public domain, its preprocessing stages are fully described in its ATBD.<fnref n="3" /> To summarise the latter, the Sentinel-3 (A and B) OLCI and SLSTR (nadir) L1b data products are coregistered and resampled to a common spatial resolution using the open-source Level 1C SYN tool.<fnref n="4" /> Pixel identification and reprojection processing are performed on the resulting product, referred to as the L1C SYN product.

Further, an atmospheric correction based on SMAC (Rahman and Dedieu, 1994) – and auxiliary meteorological data from the MERRA-2 database<fnref n="5" /> – is applied to the OLCI and SLSTR L1C SYN TOA radiances to obtain Top of Canopy (TOC) surface reflectances. Full details on the atmospheric correction can be found in the Copernicus Global Land Service ATBD of the atmospheric correction applied to Sentinel-3 data.

The Sentinel-3 TOC surface reflectance products are defined on a global 333 m WGS 84 grid divided into 10° × 10° tiles that follow the PROBA-V grid tiling system. These tiles are distributed as NetCDF files. TOC reflectances and uncertainties (associated with atmospheric correction procedure) are provided in each file for 20 Sentinel-3 spectral bands (15 OLCI bands, 5 SLSTR bands). Moreover, four quality information flags and geometry layers for OLCI and SLSTR-nadir are also included in the NetCDF files. For full details on these datasets, we refer the reader to the aforementioned ATBD.

### 2.2. BRDF prior database (MCD43P)

Numerical inversion of kernel-driven semi-empirical models is typically used to approximate the land surface BRDF (for a review, see Roujean, 2017). The linear least-squares method is commonly used to solve the BRDF inversion problem. Sometimes, the solution may prove unstable and physically unacceptable (e.g. negative values associated with BRDF descriptors and BRDF-adjusted surface reflectances). This occurs when the BRDF inversion is undetermined due to insufficient angular sampling, the presence of outliers or other sources of uncertainty in the surface reflectance observations, which leads to an ill-posed

<footnote n="3">³ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_S3-AC-V1_11_30.pdf</footnote>

<footnote n="4">⁴ https://github.com/bcdev/l1c-syn-tool</footnote>

<footnote n="5">⁵ https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/</footnote>
problem (Combal et al., 2003; Kimes et al., 2000). The ReBeLS algorithm uses a climatology of BRDF descriptors (k<sub>iso</sub> - Isotropic kernel weight; k<sub>vol</sub> - Volumetric kernel weight; k<sub>geo</sub> - Geometric kernel weight) built from MODIS MCD43A1<sup>6</sup> and MCD43A2<sup>7</sup> products version V006 (Schaaf et al., 2011; Schaaf et al., 2002; Wang et al., 2018) to constrain the results of the BRDF inversion, as in GlobAlbedo (Muller et al., 2013). The MCD43P dataset is here used as a priori parameter distribution that emphasises solutions that differ from climatology but that may yield marginally better fits to the observations.

As a first step in assembling the MCD43P prior database, MCD43A1 and MCD43A2 products for each Monday between January 1, 2013, and December 31, 2019, are downloaded from the MOTA LP DAAC server<sup>8</sup>. Then, we converted MCD43A1 and MCD43A2 tiles from MODIS (sinusoidal, 500 m) to Sentinel-3 OLCI (Plate-Carre, 333 m) grid projection and pixel size using GDAL (GDAL/OGR contributors, 2022). Afterwards, good-quality BRDF descriptors are screened by decoding the BRDF_Albedo_Band_Quality (for all bands) layer in the MCD43A2 product. The BRDF descriptors' mean and variance are computed as in GlobAlbedo<sup>9</sup> (Lewis et al., 2012b; Muller et al., 2011, 2013), taking into account statistical weights associated with the qualitative quality assurance flag of the MCD43A1 products. Lastly, linear interpolation is applied to gap-fill the MCD43P climatology. It should be emphasised that this final step is applied to all tiles, although predominantly affecting those at the highest latitudes. Hereafter, the MCD43P dataset is referred throughout the article as the dynamic prior. The term “dynamic” in this context means that the prior BRDF descriptors are not constant throughout the year and it should not be interpreted as regularly updated. Thus, once the dynamic prior has been assembled, the operational status of MODIS (or any other sensor used to assemble the dynamic prior) does not pose any risk to ReBeLS processing chains.

## 2.3. PROBA-V

Time series of atmospherically corrected surface reflectance for the Aerosol Robotic NETwork (AERONET)<sup>10</sup> sample were extracted from the daily synthesis of 333 m PROBA-V TOC product (Dierckx et al., 2014) from July 2018 to June 2019. The quality of the PROBA-V atmospherically corrected surface reflectances was screened using the Status Map (SM) layer included in the products and following the recommendations given in the PROBA-V product manual (Wolters et al., 2018). We only keep pixels when SM denotes a clear measurement associated with land and good radiometry quality in the NIR channel.

## 2.4. VIIRS

The VIIRS VNP09GA V001<sup>11</sup> product provides atmospherically corrected surface reflectances derived from the VIIRS sensor onboard the Suomi National Polar-Orbiting (Suomi NPP) platform. We consider the VNP09GA v001 product as a reference dataset to cross-validate the BRDF model retrieved from Sentinel-3 OLCI surface reflectance measurements - see section 4.1.1. Since the Sentinel-3 BRDF-descriptors are retrieved for each Monday of the period ranging from July 2018 to June 2019, we downloaded VNP09GA tiles comprising the spatial extent of Europe for each Monday of the period mentioned above from the LP DAAC server.<sup>12</sup> Further, using GDAL (GDAL/OGR contributors, 2022), we converted the VNP09GA tiles from sinusoidal to PlateCaree projection and resampled its original spatial resolution of 500 m to 333 m to

match Sentinel-3 OLCI data. Following the coding of the VNP09GA v001 quality flags presented in Table 10 of the VNP09 ATBD (see Roger et al., 2016), we rejected pixels if the quality flag QF1 indicates confident cloudy or night. In addition, the radiometric quality in spectral bands I1 and I2 should be categorised in the quality flag QF6 as good for pixels to be considered in the assembled reference dataset.

The VIIRS VNP43IA4 v001 nadir BRDF-adjusted reflectance (NBAR) product<sup>13</sup> derived from VIIRS surface reflectance observations is used in our analyses to explore the agreement of our Sentinel-3 BRDF-adjusted reflectance estimates with a reference BRDF-adjusted reflectance product. We applied the same methodology as described above to transform the original VNP43IA4 data into the spatial grid and pixel size of the Sentinel-3 data. By interpreting the quality layer in the VNP43IA4 product, we consider only those pixels for which the BRDF_Albedo_Band_Mandatory_Quality_I2 indicates full BRDF inversions.

# 3. Methodology

## 3.1. BRDF model inversion

The bidirectional surface reflectance (ρ) measured with the Sentinel-3 optical sensors can be approximated as the linear combination of three angular functions (BRDF kernels - f<sub>iso</sub>, f<sub>vol</sub>, f<sub>geo</sub>) weighted by the BRDF descriptors (k<sub>iso</sub>, k<sub>vol</sub>, k<sub>geo</sub>),

$$ \rho(\theta_s, \theta_v, \phi_r, \lambda) = k_{iso}(\lambda)f_{iso}(\theta_{SZA}, \theta_{VZA}, \phi_r) + k_{vol}(\lambda)f_{vol}(\theta_{SZA}, \theta_{VZA}, \phi_r) + k_{geo}(\lambda)f_{geo}(\theta_{SZA}, \theta_{VZA}, \phi_r) \quad (1) $$

where λ represents wavelength, θ<sub>SZA</sub> and θ<sub>VZA</sub> are the Sun and viewing (sensor) zenith angles, respectively. The relative azimuth angle (φ<sub>r</sub> = φ<sub>VAA</sub> - φ<sub>SAA</sub>) is the difference between viewing azimuth angle (φ<sub>VAA</sub>) and Sun azimuth angle (φ<sub>SAA</sub>). The BRDF kernels are approximations made from physically-based models (i.e. radiative transfer theory) that aim to retain some physical interpretation while allowing efficient numerical inversion. f<sub>iso</sub> is an isotropic function that accounts for bidirectional reflectance while observing at nadir viewing and the overhead Sun (f<sub>iso</sub> = 1). f<sub>vol</sub> is a function that replicates the surface's volume scattering characteristics (Ross, 1981; Roujean et al., 1992a) and f<sub>geo</sub> is a geometric function that accounts for the impact of shadows and the geometrical structure of surface protuberances (Li and Strahler, 1985).

MODIS has been observing the Earth since 2000, providing the remote sensing community with a large variety of Earth observation products. Among those, the MCD43 suite of products arises as the longest temporal record of land surface BRDF and albedo estimates at sub-km spatial resolution. The MCD43 algorithm assumes that a linear combination of the RossThickLiSparseReciprocal (commonly known as “MODIS kernels”) angular functions can approximate the land surface reflectance. Full expressions for these BRDF kernels can be found in the literature (Muller et al., 2013; Wanner et al., 1995). Despite the limitations of the MODIS family of kernels to accurately reproduce the hot-spot effect (Maignan et al., 2004), it becomes natural to exploit the long-term time series of MCD43 BRDF descriptors to calculate a climatology that could serve as prior information on the land surface BRDF. For this reason, ReBeLS adopts the same family of BRDF kernels used in the MCD43 product to approximate the land surface reflectance.

The wavelength-dependent BRDF descriptors (k<sub>iso</sub>, k<sub>vol</sub>, k<sub>geo</sub>) from Eq. (1) are then found via solving the BRDF inversion problem, which can be stated as follows: What is the optimal BRDF model that represents the observed surface reflectances acquired under different illumination and viewing conditions? The solution to the BRDF inverse problem requires determining the wavelength-dependent BRDF descriptors using an optimal estimation approach.

Since the BRDF inversion is often an undetermined problem, noisy

<footnote n="6">6 https://lpdaac.usgs.gov/products/mcd43a1v006/</footnote>

<footnote n="7">7 https://lpdaac.usgs.gov/products/mcd43a2v006/</footnote>

<footnote n="8">8 https://e4ftl01.cr.usgs.gov/MOTA/</footnote>

<footnote n="9">9 http://www.globalbedo.org/</footnote>

<footnote n="10">10 https://aeronet.gsfc.nasa.gov/</footnote>

<footnote n="11">11 https://lpdaac.usgs.gov/products/vnp09gav001/</footnote>

<footnote n="12">12 https://e4ftl01.cr.usgs.gov/</footnote>

<footnote n="13">13 https://lpdaac.usgs.gov/products/vnp43ia4v001/</footnote>
problem (Combal et al., 2003; Kimes et al., 2000). The ReBeLS algorithm uses a climatology of BRDF descriptors (k<sub>iso</sub> - Isotropic kernel weight; k<sub>vol</sub> - Volumetric kernel weight; k<sub>geo</sub> - Geometric kernel weight) built from MODIS MCD43A1<sup>6</sup> and MCD43A2<sup>7</sup> products version V006 (Schaaf et al., 2011; Schaaf et al., 2002; Wang et al., 2018) to constrain the results of the BRDF inversion, as in GlobAlbedo (Muller et al., 2013). The MCD43P dataset is here used as a priori parameter distribution that emphasises solutions that differ from climatology but that may yield marginally better fits to the observations.

As a first step in assembling the MCD43P prior database, MCD43A1 and MCD43A2 products for each Monday between January 1, 2013, and December 31, 2019, are downloaded from the MOTA LP DAAC server<sup>8</sup>. Then, we converted MCD43A1 and MCD43A2 tiles from MODIS (sinusoidal, 500 m) to Sentinel-3 OLCI (Plate-Carre, 333 m) grid projection and pixel size using GDAL (GDAL/OGR contributors, 2022). Afterwards, good-quality BRDF descriptors are screened by decoding the BRDF_Albedo_Band_Quality (for all bands) layer in the MCD43A2 product. The BRDF descriptors' mean and variance are computed as in GlobAlbedo<sup>9</sup> (Lewis et al., 2012b; Muller et al., 2011, 2013), taking into account statistical weights associated with the qualitative quality assurance flag of the MCD43A1 products. Lastly, linear interpolation is applied to gap-fill the MCD43P climatology. It should be emphasised that this final step is applied to all tiles, although predominantly affecting those at the highest latitudes. Hereafter, the MCD43P dataset is referred throughout the article as the dynamic prior. The term "dynamic" in this context means that the prior BRDF descriptors are not constant throughout the year and it should not be interpreted as regularly updated. Thus, once the dynamic prior has been assembled, the operational status of MODIS (or any other sensor used to assemble the dynamic prior) does not pose any risk to ReBeLS processing chains.

## 2.3. PROBA-V

Time series of atmospherically corrected surface reflectance for the Aerosol Robotic NETwork (AERONET)<sup>10</sup> sample were extracted from the daily synthesis of 333 m PROBA-V TOC product (Dierckx et al., 2014) from July 2018 to June 2019. The quality of the PROBA-V atmospherically corrected surface reflectances was screened using the Status Map (SM) layer included in the products and following the recommendations given in the PROBA-V product manual (Wolters et al., 2018). We only keep pixels when SM denotes a clear measurement associated with land and good radiometry quality in the NIR channel.

## 2.4. VIIRS

The VIIRS VNP09GA V001<sup>11</sup> product provides atmospherically corrected surface reflectances derived from the VIIRS sensor onboard the Suomi National Polar-Orbiting (Suomi NPP) platform. We consider the VNP09GA v001 product as a reference dataset to cross-validate the BRDF model retrieved from Sentinel-3 OLCI surface reflectance measurements - see section 4.1.1. Since the Sentinel-3 BRDF-descriptors are retrieved for each Monday of the period ranging from July 2018 to June 2019, we downloaded VNP09GA tiles comprising the spatial extent of Europe for each Monday of the period mentioned above from the LP DAAC server.<sup>12</sup> Further, using GDAL (GDAL/OGR contributors, 2022), we converted the VNP09GA tiles from sinusoidal to PlateCaree projection and resampled its original spatial resolution of 500 m to 333 m to

match Sentinel-3 OLCI data. Following the coding of the VNP09GA v001 quality flags presented in Table 10 of the VNP09 ATBD (see Roger et al., 2016), we rejected pixels if the quality flag QF1 indicates confident cloudy or night. In addition, the radiometric quality in spectral bands I1 and I2 should be categorised in the quality flag QF6 as good for pixels to be considered in the assembled reference dataset.

The VIIRS VNP43IA4 v001 nadir BRDF-adjusted reflectance (NBAR) product<sup>13</sup> derived from VIIRS surface reflectance observations is used in our analyses to explore the agreement of our Sentinel-3 BRDF-adjusted reflectance estimates with a reference BRDF-adjusted reflectance product. We applied the same methodology as described above to transform the original VNP43IA4 data into the spatial grid and pixel size of the Sentinel-3 data. By interpreting the quality layer in the VNP43IA4 product, we consider only those pixels for which the BRDF_Albedo_Band_Mandatory_Quality_I2 indicates full BRDF inversions.

# 3. Methodology

## 3.1. BRDF model inversion

The bidirectional surface reflectance (ρ) measured with the Sentinel-3 optical sensors can be approximated as the linear combination of three angular functions (BRDF kernels - f<sub>iso</sub>, f<sub>vol</sub>, f<sub>geo</sub>) weighted by the BRDF descriptors (k<sub>iso</sub>, k<sub>vol</sub>, k<sub>geo</sub>),

$$ \rho(\theta_s, \theta_v, \phi_r, \lambda) = k_{iso}(\lambda)f_{iso}(\theta_{SZA}, \theta_{VZA}, \phi_r) + k_{vol}(\lambda)f_{vol}(\theta_{SZA}, \theta_{VZA}, \phi_r) + k_{geo}(\lambda)f_{geo}(\theta_{SZA}, \theta_{VZA}, \phi_r) \quad (1) $$

where λ represents wavelength, θ<sub>SZA</sub> and θ<sub>VZA</sub> are the Sun and viewing (sensor) zenith angles, respectively. The relative azimuth angle (φ<sub>r</sub> = φ<sub>VAA</sub> - φ<sub>SAA</sub>) is the difference between viewing azimuth angle (φ<sub>VAA</sub>) and Sun azimuth angle (φ<sub>SAA</sub>). The BRDF kernels are approximations made from physically-based models (i.e. radiative transfer theory) that aim to retain some physical interpretation while allowing efficient numerical inversion. f<sub>iso</sub> is an isotropic function that accounts for bidirectional reflectance while observing at nadir viewing and the overhead Sun (f<sub>iso</sub> = 1). f<sub>vol</sub> is a function that replicates the surface's volume scattering characteristics (Ross, 1981; Roujean et al., 1992a) and f<sub>geo</sub> is a geometric function that accounts for the impact of shadows and the geometrical structure of surface protuberances (Li and Strahler, 1985).

MODIS has been observing the Earth since 2000, providing the remote sensing community with a large variety of Earth observation products. Among those, the MCD43 suite of products arises as the longest temporal record of land surface BRDF and albedo estimates at sub-km spatial resolution. The MCD43 algorithm assumes that a linear combination of the RossThickLiSparseReciprocal (commonly known as "MODIS kernels") angular functions can approximate the land surface reflectance. Full expressions for these BRDF kernels can be found in the literature (Muller et al., 2013; Wanner et al., 1995). Despite the limitations of the MODIS family of kernels to accurately reproduce the hot-spot effect (Maignan et al., 2004), it becomes natural to exploit the long-term time series of MCD43 BRDF descriptors to calculate a climatology that could serve as prior information on the land surface BRDF. For this reason, ReBeLS adopts the same family of BRDF kernels used in the MCD43 product to approximate the land surface reflectance.

The wavelength-dependent BRDF descriptors (k<sub>iso</sub>, k<sub>vol</sub>, k<sub>geo</sub>) from Eq. (1) are then found via solving the BRDF inversion problem, which can be stated as follows: What is the optimal BRDF model that represents the observed surface reflectances acquired under different illumination and viewing conditions? The solution to the BRDF inverse problem requires determining the wavelength-dependent BRDF descriptors using an optimal estimation approach.

Since the BRDF inversion is often an undetermined problem, noisy

<footnote n="6"><sup>6</sup> https://lpdaac.usgs.gov/products/mcd43a1v006/</footnote>

<footnote n="7"><sup>7</sup> https://lpdaac.usgs.gov/products/mcd43a2v006/</footnote>

<footnote n="8"><sup>8</sup> https://e4ftl01.cr.usgs.gov/MOTA/</footnote>

<footnote n="9"><sup>9</sup> http://www.globalbedo.org/</footnote>

<footnote n="10"><sup>10</sup> https://aeronet.gsfc.nasa.gov/</footnote>

<footnote n="11"><sup>11</sup> https://lpdaac.usgs.gov/products/vnp09gav001/</footnote>

<footnote n="12"><sup>12</sup> https://e4ftl01.cr.usgs.gov/</footnote>

<footnote n="13"><sup>13</sup> https://lpdaac.usgs.gov/products/vnp43ia4v001/</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Diagram illustrating the processing flow for Sentinel-3 BRDF modelling using ReBeLS.</figcaption></figure>

The steps within the Inversion and Outputs (BAR) stages are described in Appendix B and section 3.2, respectively. It should be noticed that the processing is independently applied to each pixel and spectral band.

observations and outliers (undetected clouds or shadows, etc.) can result in physically questionable solutions. A Bayesian context therefore allows us to introduce constraints (prior knowledge of the model temporal evolution and parameter values) to stabilize the inversion procedure and to quantify the effect of these constraints in the uncertainty estimation. In the ReBeLS BRDF inverse problem formulation, the measurements (surface reflectance observations) are assumed to be corrupted by additive zero mean Gaussian noise. The a priori probability density function (PDF) for the BRDF descriptors is assumed Gaussian. The goal is to estimate the a posteriori PDF of the BRDF descriptors. Given the linearity of the BRDF model and the use of Gaussian statistics (Tarantola, 2004), the posterior PDF will also be a multivariate Gaussian (defined by a posterior mean vector and a covariance matrix):

$$P(\mathbf{x}|\mathbf{R}) \sim \exp\left(-\frac{1}{2}(\mathbf{x} - \mathbf{x}_{\text{posterior}})^T \mathbf{C}_{\text{posterior}}^{-1} (\mathbf{x} - \mathbf{x}_{\text{posterior}})\right) \quad (2)$$

where $x_{posterior}$ and $C_{posterior}^{-1}$ are the mean and covariance matrix of the maximum likelihood estimator of $P(x)$ and are obtained by minimising a cost function $J(x)$ (taken to be the log of the posterior PDF expressed as follows:

$$\begin{aligned} J(\mathbf{x}) &= J_{\text{obs}}(\mathbf{x}) + J_{\text{prior}}(\mathbf{x}) + J_{\text{model}}(\mathbf{x}) \\ &= \frac{1}{2} [\mathbf{R} - H(\mathbf{x}, \Omega_{\text{sensor}}, \Omega_{\text{Sun}})]^T \mathbf{C}_{\text{obs}}^{-1} [\mathbf{R} - H(\mathbf{x}, \Omega_{\text{sensor}}, \Omega_{\text{Sun}})] \\ &\quad + \frac{1}{2} [\mathbf{x} - \mathbf{x}_{\text{prior}}]^T \mathbf{C}_{\text{prior}}^{-1} [\mathbf{x} - \mathbf{x}_{\text{prior}}] \\ &\quad + \frac{1}{2} [\mathbf{x} - \mathbf{M}(\mathbf{x})]^T \mathbf{C}_{\text{model}}^{-1} [\mathbf{x} - \mathbf{M}(\mathbf{x})] \end{aligned} \quad (3)$$

The terms on the right-hand side of the cost functions can be understood as follows: $J_{obs}(x)$ measures the mismatch between the model's predicted values $H(x, \Omega_{sensor}, \Omega_{Sun})$ of the surface reflectance and the observed ones ($R$). The kernel-driven semi-empirical RossThickLS-SparseReciprocal BRDF model is represented by $H(x, \Omega_{sensor}, \Omega_{Sun})$ where $\Omega_{sensor}$ and $\Omega_{Sun}$ are the sensor ($\theta_{VZA}, \phi_{VAA}$) and Sun geometries ($\theta_{SZA}, \phi_{SAA}$), respectively. $J_{prior}(x)$ represents the constraint introduced by the prior knowledge of the BRDF descriptors $x_{prior}$ and their associated uncertainties $C_{prior}^{-1}$ obtained from the auxiliary MCD43P dataset. Since the wide imaging swaths of the Sentinel-3 instruments have the potential to enable a successful BRDF inversion without requiring prior information (if enough cloud-free observations with sufficient angular diversity are available over a non-changing surface), one might wonder about the need of a dynamic prior of BRDF-descriptors constraint. However, our experience with sensors such as SPOT-VGT, MERIS and PROBA-V (León-Tavares et al., 2021; Lewis et al., 2012b; Muller et al.,

2013; Roujean et al., 2018), which have similar imaging swath widths and orbits as Sentinel-3, suggests that using a dynamic prior constraint helps to (i) stabilize inversions in periods with few observations, (ii) dampen the effect of outliers due to unfiltered clouds or cloud shadows and (iii) introducing some realistic temporal dynamics during long periods of persistent cloud cover. These remarks are acknowledged by other operational BRDF processing chains, such as LandSaf (Carrer et al., 2021; Carrer et al., 2018; Geiger et al., 2008), where despite the fact that the MSG-SEVIRI provides a good sampling of the BRDF (as depicted in Fig. 3 of the LandSaf albedo ATBD<fnref n="14" />), their processing chain still requires prior estimates to stabilize the kernel parameter inversion.

$J_{model}(X)$ in Eq. (3) is a constraint on the temporal evolution of the BRDF descriptors where $M$ is the dynamic model. Using a first-order difference model represented by matrix $B$

$$B = \begin{bmatrix} 1 & -1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & -1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 & 0 \\ 0 & 0 & 0 & \cdots & 1 & -1 \end{bmatrix} \quad (4)$$

the constraint to the temporal evolution of the state vector $x$ takes the form $J_{model}(x) = \gamma x^T (B^T B)x$. Here, the first-order difference model matrix is applied at a lag of one day to exert the constraint that the solution for tomorrow will be the same as the solution of today (Gómez-Dans et al., 2016; Lewis et al., 2012a; Quaife and Lewis, 2010; Zobitz et al., 2020). The term $\gamma$ is the regularisation parameter representing the uncertainty on the assumption that the difference of the BRDF descriptors between consecutive days should be zero (Quaife and Lewis, 2010). The value of the regularisation parameter is deemed constant within ReBeLS, and section 6 (Appendix A) describes the methodology implemented to find it.

The main ReBeLS algorithm described above is applied to time series of individual Sentinel-3 pixels, with processing conducted independently for each spectral band. A schematic representation of the processing sequence is delineated in Fig. 1. For a more comprehensive description of the Inversion stage and its intermediate steps (A, B1-B3 in Fig. 1), we refer the reader to section Appendix B. Subsequently, in the following subsection we elaborate on the BAR stage. In summary, the ReBeLS BRDF retrieval methodology should be understood as a statistical blending where the goodness of fit to the observations is maximized simultaneously by penalizing high-frequency components of the temporal BRDF descriptors trajectories and solutions departing from the prior estimates. Finally, ReBeLS outputs the solutions at the cadence

<footnote n="14">¹⁴ https://nextcloud.lsasvc.ipma.pt/s/SxANfxgaw5sSB63 accessed on 24April2023</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Diagram illustrating the processing flow for Sentinel-3 BRDF modelling using ReBeLS.</figcaption></figure>

The steps within the Inversion and Outputs (BAR) stages are described in Appendix B and section 3.2, respectively. It should be noticed that the processing is independently applied to each pixel and spectral band.

observations and outliers (undetected clouds or shadows, etc.) can result in physically questionable solutions. A Bayesian context therefore allows us to introduce constraints (prior knowledge of the model temporal evolution and parameter values) to stabilize the inversion procedure and to quantify the effect of these constraints in the uncertainty estimation. In the ReBeLS BRDF inverse problem formulation, the measurements (surface reflectance observations) are assumed to be corrupted by additive zero mean Gaussian noise. The a priori probability density function (PDF) for the BRDF descriptors is assumed Gaussian. The goal is to estimate the a posteriori PDF of the BRDF descriptors. Given the linearity of the BRDF model and the use of Gaussian statistics (Tarantola, 2004), the posterior PDF will also be a multivariate Gaussian (defined by a posterior mean vector and a covariance matrix):

$$P(\mathbf{x}|\mathbf{R}) \sim \exp\left(-\frac{1}{2}(\mathbf{x} - \mathbf{x}_{\text{posterior}})^T \mathbf{C}_{\text{posterior}}^{-1} (\mathbf{x} - \mathbf{x}_{\text{posterior}})\right) \quad (2)$$

where $x_{posterior}$ and $C_{posterior}^{-1}$ are the mean and covariance matrix of the maximum likelihood estimator of $P(x)$ and are obtained by minimising a cost function $J(x)$ (taken to be to the log of the posterior PDF expressed as follows:

$$\begin{aligned} J(\mathbf{x}) &= J_{\text{obs}}(\mathbf{x}) + J_{\text{prior}}(\mathbf{x}) + J_{\text{model}}(\mathbf{x}) \\ &= \frac{1}{2} [\mathbf{R} - H(\mathbf{x}, \Omega_{\text{sensor}}, \Omega_{\text{Sun}})]^T \mathbf{C}_{\text{obs}}^{-1} [\mathbf{R} - H(\mathbf{x}, \Omega_{\text{sensor}}, \Omega_{\text{Sun}})] \\ &\quad + \frac{1}{2} [\mathbf{x} - \mathbf{x}_{\text{prior}}]^T \mathbf{C}_{\text{prior}}^{-1} [\mathbf{x} - \mathbf{x}_{\text{prior}}] \\ &\quad + \frac{1}{2} [\mathbf{x} - \mathbf{M}(\mathbf{x})]^T \mathbf{C}_{\text{model}}^{-1} [\mathbf{x} - \mathbf{M}(\mathbf{x})] \end{aligned} \quad (3)$$

The terms on the right-hand side of the cost functions can be understood as follows: $J_{obs}(x)$ measures the mismatch between the model's predicted values $H(x, \Omega_{sensor}, \Omega_{Sun})$ of the surface reflectance and the observed ones ($R$). The kernel-driven semi-empirical RossThickLS-SparseReciprocal BRDF model is represented by $H(x, \Omega_{sensor}, \Omega_{Sun})$ where $\Omega_{sensor}$ and $\Omega_{Sun}$ are the sensor ($\theta_{VZA}, \phi_{VAA}$) and Sun geometries ($\theta_{SZA}, \phi_{SAA}$), respectively. $J_{prior}(x)$ represents the constraint introduced by the prior knowledge of the BRDF descriptors $x_{prior}$ and their associated uncertainties $C_{prior}^{-1}$ obtained from the auxiliary MCD43P dataset. Since the wide imaging swaths of the Sentinel-3 instruments have the potential to enable a successful BRDF inversion without requiring prior information (if enough cloud-free observations with sufficient angular diversity are available over a non-changing surface), one might wonder about the need of a dynamic prior of BRDF-descriptors constraint. However, our experience with sensors such as SPOT-VGT, MERIS and PROBA-V (León-Tavares et al., 2021; Lewis et al., 2012b; Muller et al.,

2013; Roujean et al., 2018), which have similar imaging swath widths and orbits as Sentinel-3, suggests that using a dynamic prior constraint helps to (i) stabilize inversions in periods with few observations, (ii) dampen the effect of outliers due to unfiltered clouds or cloud shadows and (iii) introducing some realistic temporal dynamics during long periods of persistent cloud cover. These remarks are acknowledged by other operational BRDF processing chains, such as LandSaf (Carrer et al., 2021; Carrer et al., 2018; Geiger et al., 2008), where despite the fact that the MSG-SEVIRI provides a good sampling of the BRDF (as depicted in Fig. 3 of the LandSaf albedo ATBD<fnref n="14" />), their processing chain still requires prior estimates to stabilize the kernel parameter inversion.

$J_{model}(X)$ in Eq. (3) is a constraint on the temporal evolution of the BRDF descriptors where $M$ is the dynamic model. Using a first-order difference model represented by matrix $B$

$$B = \begin{bmatrix} 1 & -1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & -1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & -1 & 0 \\ 0 & 0 & 0 & \cdots & 1 & -1 \end{bmatrix} \quad (4)$$

the constraint to the temporal evolution of the state vector $x$ takes the form $J_{model}(x) = \gamma x^T (B^T B)x$. Here, the first-order difference model matrix is applied at a lag of one day to exert the constraint that the solution for tomorrow will be the same as the solution of today (Gómez-Dans et al., 2016; Lewis et al., 2012a; Quaife and Lewis, 2010; Zobitz et al., 2020). The term $\gamma$ is the regularisation parameter representing the uncertainty on the assumption that the difference of the BRDF descriptors between consecutive days should be zero (Quaife and Lewis, 2010). The value of the regularisation parameter is deemed constant within ReBeLS, and section 6 (Appendix A) describes the methodology implemented to find it.

The main ReBeLS algorithm described above is applied to time series of individual Sentinel-3 pixels, with processing conducted independently for each spectral band. A schematic representation of the processing sequence is delineated in Fig. 1. For a more comprehensive description of the Inversion stage and its intermediate steps (A, B1-B3 in Fig. 1), we refer the reader to section Appendix B. Subsequently, in the following subsection we elaborate on the BAR stage. In summary, the ReBeLS BRDF retrieval methodology should be understood as a statistical blending where the goodness of fit to the observations is maximized simultaneously by penalizing high-frequency components of the temporal BRDF descriptors trajectories and solutions departing from the prior estimates. Finally, ReBeLS outputs the solutions at the cadence

<footnote n="14">¹⁴ https://nextcloud.lsasvc.ipma.pt/s/SxANfxgaw5sSB63 accessed on 24April2023</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. RGB image of the Oa17 BRDF descriptors retrieved for 3rd June 2019. The colour channels are as follows: Red - kiso, green - kvol and blue kgeo.</figcaption></figure>

desired. For the examples in section 4, the output cadence is set to each Monday but could be daily (see Fig. 7) or dekadal (Sánchez-Zapero et al., 2023). ReBeLS does not depend on a fixed composite window period and provides uncertainties associated with the BRDF descriptors. ReBeLS allows users to choose the most appropriate accumulation period for their needs rather than relying on a fixed composite window. An example of the benefits of this feature is demonstrated by the varying composite windows used in the production of the CGLS NDVI product<fnref n="15" />, which is computed from ReBeLS BRDF-adjusted surface reflectances and the ReBeLS-derived C3S surface albedo product.<fnref n="16" /> The former employs a 30-day composite window, while the latter uses composite windows of 30 days and 365 days for near-real-time and back-processing products, respectively.

Moreover, ReBeLS does not require a backup algorithm for periods

where the input surface reflectance time series is sparsely populated. However, it is worth mentioning that the solutions obtained during periods of poor data availability will be dominated by the dynamic prior (see section 2.2), which translates into a greater uncertainty in the retrieved BRDF descriptors.

## 3.2. BRDF-adjustment

Once ReBeLS has retrieved the BRDF descriptors – Outputs (CBD) in Fig. 1, it follows the computation of the BRDF-adjusted reflectance. To start with, ReBeLS computes the solar-zenith angle at local 10 AM ($\theta_{SZA_{10am}}$) for each pixel and sets the normalisation VZA to nadir ($\theta_{VZA} = 0^{\circ}$) – step C1 in Fig. 1. The BRDF-adjusted reflectance (also indistinctly referred to as BAR) is obtained – step C2 in Fig. 1 – using the retrieved BRDF descriptors and the RossThickLiSparse kernels computed for the normalisation geometry. This work computes BRDF adjusted reflectances at nadir viewing and solar zenith angle at local 10 AM as follows:

$$ BAR = k_{iso} + k_{vol}f_{vol}(\theta_{SZA_{10am}}, 0, 0) + k_{geo}f_{geo}(\theta_{SZA_{10am}}, 0, 0) \quad (5) $$

<footnote n="15">¹⁵ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_NDVI300m-V2_11.20.pdf</footnote>

<footnote n="16">¹⁶ https://datastore.copernicus-climate.eu/documents/satellite-albedo/C3S_COP_059_D-02_ATBD_CDR-ICDR_SA_SENTINEL3_v3.0_PRODUCTs_v1.1.pdf</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. RGB image of the Oa17 BRDF descriptors retrieved for 3rd June 2019. The colour channels are as follows: Red - kiso, green - kvol and blue kgeo.</figcaption></figure>

desired. For the examples in section 4, the output cadence is set to each Monday but could be daily (see Fig. 7) or dekadal (Sánchez-Zapero et al., 2023). ReBeLS does not depend on a fixed composite window period and provides uncertainties associated with the BRDF descriptors. ReBeLS allows users to choose the most appropriate accumulation period for their needs rather than relying on a fixed composite window. An example of the benefits of this feature is demonstrated by the varying composite windows used in the production of the CGLS NDVI product<fnref n="15" />, which is computed from ReBeLS BRDF-adjusted surface reflectances and the ReBeLS-derived C3S surface albedo product.<fnref n="16" /> The former employs a 30-day composite window, while the latter uses composite windows of 30 days and 365 days for near-real-time and back-processing products, respectively.

Moreover, ReBeLS does not require a backup algorithm for periods

where the input surface reflectance time series is sparsely populated. However, it is worth mentioning that the solutions obtained during periods of poor data availability will be dominated by the dynamic prior (see section 2.2), which translates into a greater uncertainty in the retrieved BRDF descriptors.

## 3.2. BRDF-adjustment

Once ReBeLS has retrieved the BRDF descriptors – Outputs (CBD) in Fig. 1, it follows the computation of the BRDF-adjusted reflectance. To start with, ReBeLS computes the solar-zenith angle at local 10 AM ($\theta_{SZA_{10am}}$) for each pixel and sets the normalisation VZA to nadir ($\theta_{VZA} = 0°$) – step C1 in Fig. 1. The BRDF-adjusted reflectance (also indistinctly referred to as BAR) is obtained – step C2 in Fig. 1 – using the retrieved BRDF descriptors and the RossThickLiSparse kernels computed for the normalisation geometry. This work computes BRDF adjusted reflectances at nadir viewing and solar zenith angle at local 10 AM as follows:

$$ BAR = k_{iso} + k_{vol}f_{vol}(\theta_{SZA_{10am}}, 0, 0) + k_{geo}f_{geo}(\theta_{SZA_{10am}}, 0, 0) \quad (5) $$

<footnote n="15">¹⁵ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_NDVI300m-V2_11.20.pdf</footnote>

<footnote n="16">¹⁶ https://datastore.copernicus-climate.eu/documents/satellite-albedo/C3S_COP_059_D-02_ATBD_CDR-ICDR_SA_SENTINEL3_v3.0_PRODUCTS_v1.1.pdf</footnote>
<table><caption>**Table 1**

List of stations used to assess the performance of the BRDF modelling.</caption><thead><tr><th>Site_ID</th><th>Site_name</th><th>IGBP classification</th><th>Site details</th></tr></thead><tbody><tr><td>PT-Mi1</td><td>Mitra II (Evora)</td><td>Ever Broadleaf Forests (EBF)</td><td>http://www.europe-fluxdata.eu/home/site-details?id=PT-Mi1</td></tr><tr><td>IT-Isp</td><td>Ispra ABC-IS</td><td>Deciduous Broadleaf forests (DBF)</td><td>http://www.europe-fluxdata.eu/home/site-details?id=IT-Isp</td></tr><tr><td>DE-Hai</td><td>Hainich</td><td>DBF</td><td>http://www.europe-fluxdata.eu/home/site-details?id=DE-Hai</td></tr></tbody></table>

For a full description of the IGBP classes, see https://fluxnet.org/data/badm-dat-a-templates/igbp-classification/.

By assuming there is no covariance between BRDF descriptors, we propagate the BRDF retrievals uncertainties through the BAR calculations to estimate the uncertainty associated with the BAR computations. To this end, we make use of the following expression.

$$ \sigma BAR = \sqrt{\sigma_{kiso}^2 + \sigma_{kvol}^2 f_{vol}^2(\theta_{SZA_{10am}}, 0, 0) + \sigma_{kgeo}^2 f_{geo}^2(\theta_{SZA_{10am}}, 0, 0)} \quad (6) $$

## 4. Results and discussion

Here we evaluate the BRDF descriptors retrievals obtained with ReBeLS by inspecting their spatial consistency and ability to capture the intrinsic land surface reflectance anisotropy. We assembled mosaics from the ReBeLS BAR product and compared them against mosaics assembled with the MVC technique and those derived from the VNP43 A4 products. To this end, we ran ReBeLS at a European scale using surface reflectance observations acquired with Sentinel-3 A and Sentinel-3 B OLCI sensors from June 2018 to June 2019. A weekly dataset of BRDF descriptors for OLCI spectral channels Oa08 ($\lambda_{centre} = 665nm$) and Oa17 ($\lambda_{centre} = 865nm$) was obtained, and for brevity, we focus our analyses on the results obtained for Oa17 spectral band. But as it can be seen in section 4.1, the discussed results hold for band Oa08 as well. It should be emphasised that the ReBeLS algorithm can process data from any optical spectral band from the SLSTR instrument, as requested by the Copernicus Climate Change Service albedo product (Sánchez-Zapero et al., 2023).

ReBeLS is designed to retrieve BRDF descriptors on a daily basis throughout a user-defined accumulation period, which, as previously mentioned, in this article has been set to one year. Subsequently, the discretion lies with the user to determine which of the daily BRDF descriptors are to be recorded as output data. It is noteworthy that in our analyses, we have opted to generate ReBeLS BRDF descriptors on a weekly schedule, aligning with the temporal frequency of the MCD43P dataset. Nevertheless, the temporal frequency for BRDF descriptor retrieval can be tailored to various intervals as per the user's requirements. For instance, it can be configured to operate on a daily basis, as exemplified in Fig. 7, or at dekadal intervals, in line with the Copernicus Global Land BRDF-corrected NDVI²³ and Copernicus Climate Change Surface Albedo²⁴ products.

### 4.1. BRDF descriptors retrievals

Fig. 2 displays an RGB image of the Oa17 BRDF descriptors corresponding to the 3rd June 2019. The colour channels are as follows: Red - $k_{iso}$, green - $k_{vol}$ and blue $k_{geo}$. The spatial consistency of the ReBeLS BRDF retrievals at the Oa17 band can be appraised from the spatial patterns coinciding with most spatial features associated with land cover which is consistent with the assumption that there is a clear connection between BRDF descriptors and land cover types belonging to the same categorical class. Any difference within the same land cover class will likely be about the magnitude of the BRDF rather than its shape (Strugnell et al., 2001). Nonetheless, based on the authors' experience, for some land cover types, variations in stem density can induce

<table><caption>**Table 2**

The number of observations used in the ReBeLS BRDF inversion and goodness of the model fits represented by the average ($<z_{sc}$>) and standard deviation ($\sigma_{z_{sc}}$) of the zeta-score for each of the stations in Table 1.</caption><thead><tr><th></th><th>Oa08</th><th>Oa17</th><th>Oa08</th><th>Oa17</th></tr><tr><th></th><th>$N_{obs}$</th><th>$N_{obs}$</th><th>$<z_{sc}$>, $\sigma_{z_{sc}}$</th><th>$<z_{sc}$>, $\sigma_{z_{sc}}$</th></tr></thead><tbody><tr><th>PT-Mi1</th><td>191</td><td>194</td><td>-0.03, 1.24</td><td>-0.07, 0.78</td></tr><tr><th>IT-Isp</th><td>170</td><td>170</td><td>-0.14, 1.49</td><td>0.11, 0.90</td></tr><tr><th>DE-Hai</th><td>102</td><td>120</td><td>0.26, 0.86</td><td>-0.07, 0.84</td></tr></tbody></table>

different BRDF shapes within the same land unit.

The ability of ReBeLS to reproduce the temporal evolution of the Sentinel-3 OLCI surface reflectance observations while maintaining a physically acceptable temporal evolution of BRDF-descriptors is examined for three European FLUXNET<fnref n="17" /> stations (see Table 1). The selected sites have latitudes ranging from 38.5° to 51.0° and exhibit temporal variability to test the temporal trajectories of the BRDF retrievals.

To estimate the goodness of the BRDF modelling, we compute the zeta-score ($z_{sc}$) as follows:

$$ z_{sc} = \frac{(\text{reflectance}_{\text{obs}} - \text{reflectance}_{\text{fwd}})}{\left( \sigma_{\text{reflectance}_{\text{obs}}}^2 + \sigma_{\text{reflectance}_{\text{fwd}}}^2 \right)^{1/2}} \quad (7) $$

where $\text{reflectance}_{\text{obs}}$ and $\text{reflectance}_{\text{fwd}}$ are the observed and model fitted (using retrieved BRDF model in forward mode) surface reflectances, and their associated uncertainties are denoted by $\sigma$. The zeta-score's absolute value indicates whether observed and model-fitted surface reflectances agree within uncertainties. A zeta-score with an absolute value above three should be regarded as a poor fit, while an absolute value of a zeta-score lower than two (three) shall be considered a good (an adequate) fit to the data. For each of the stations listed in Table 1, we report the number of observations used throughout the ReBeLS inversion, the average zeta-score value and its standard deviation for Oa08 and Oa17 spectral bands in Table 2.

The top-left panel of Fig. 3 shows the observed (open circles) and model-fitted (filled circles) Sentinel-3 surface reflectance time series in spectral band Oa17 for the pixel closest to the FLUXNET site PT-Mi1 (Portugal). The distribution of zeta-scores obtained for the Portuguese station is displayed in the bottom left panel of Fig. 3, where the colour bands are coded as follows: green denotes a good fit ($|z_{sc}| < 2$), yellow represents an adequate fit ($2 < |z_{sc}| < 3$), while points lying in the red band should be considered as a poor fit ($|z_{sc}| > 3$). The temporal trajectories of the retrieved Sentinel-3 Oa17 BRDF descriptors trajectories obtained by the ReBeLS processor are shown on the right panels (top: $k_{iso}$, middle: $k_{vol}$, bottom: $k_{geo}$) of Fig. 3. The MODIS Band2 MCD43A1 product for each Monday of the year considered is overplotted (red triangles in the right panels of Fig. 3) for comparison.

From the bottom left panel of Fig. 3, most of the model fits lie within the green zeta-score band, suggesting that our BRDF model successfully reproduces the Sentinel-3 observations of the PT-Mi1 station at spectral band Oa17. Fig. 4 displays the ReBeLS results for PT-Mi1 in the spectral band Oa08, and from the zeta-score panel, it can be seen that the large majority of the model fits are within the $\pm 2\sigma$ (green) band, although it is noticeable that few model fits are scattered across the yellow and red bands of the z-score panel which is reflected as an increase in the standard deviation of the zeta-score distribution as shown in Table 2.

The BRDF model fits and descriptor trajectories for IT-Isp and DE-Hai are shown in Figs. 5 and 6, respectively. From the top left panel of Fig. 6 and Table 2, it can be noticed that the availability of clear surface reflectance observations near the DE-Hai station is affected by a high

<footnote n="17">¹⁷ https://fluxnet.org/</footnote>
<table><caption>**Table 1**

List of stations used to assess the performance of the BRDF modelling.</caption><thead><tr><th>Site_ID</th><th>Site_name</th><th>IGBP classification</th><th>Site details</th></tr></thead><tbody><tr><td>PT-Mi1</td><td>Mitra II (Evora)</td><td>Ever Broadleaf Forests (EBF)</td><td>http://www.europe-fluxdata.eu/home/site-details?id=PT-Mi1</td></tr><tr><td>IT-Isp</td><td>Ispra ABC-IS</td><td>Deciduous Broadleaf forests (DBF)</td><td>http://www.europe-fluxdata.eu/home/site-details?id=IT-Isp</td></tr><tr><td>DE-Hai</td><td>Hainich</td><td>DBF</td><td>http://www.europe-fluxdata.eu/home/site-details?id=DE-Hai</td></tr></tbody></table>

For a full description of the IGBP classes, see https://fluxnet.org/data/badm-dat-a-templates/igbp-classification/.

By assuming there is no covariance between BRDF descriptors, we propagate the BRDF retrievals uncertainties through the BAR calculations to estimate the uncertainty associated with the BAR computations. To this end, we make use of the following expression.

$$ \sigma BAR = \sqrt{\sigma_{kiso}^2 + \sigma_{kvol}^2 f_{vol}^2(\theta_{SZA_{10am}}, 0, 0) + \sigma_{kgeo}^2 f_{geo}^2(\theta_{SZA_{10am}}, 0, 0)} \quad (6) $$

## 4. Results and discussion

Here we evaluate the BRDF descriptors retrievals obtained with ReBeLS by inspecting their spatial consistency and ability to capture the intrinsic land surface reflectance anisotropy. We assembled mosaics from the ReBeLS BAR product and compared them against mosaics assembled with the MVC technique and those derived from the VNP43 A4 products. To this end, we ran ReBeLS at a European scale using surface reflectance observations acquired with Sentinel-3 A and Sentinel-3 B OLCI sensors from June 2018 to June 2019. A weekly dataset of BRDF descriptors for OLCI spectral channels Oa08 ($\lambda_{centre} = 665nm$) and Oa17 ($\lambda_{centre} = 865nm$) was obtained, and for brevity, we focus our analyses on the results obtained for Oa17 spectral band. But as it can be seen in section 4.1, the discussed results hold for band Oa08 as well. It should be emphasised that the ReBeLS algorithm can process data from any optical spectral band from the SLSTR instrument, as requested by the Copernicus Climate Change Service albedo product (Sánchez-Zapero et al., 2023).

ReBeLS is designed to retrieve BRDF descriptors on a daily basis throughout a user-defined accumulation period, which, as previously mentioned, in this article has been set to one year. Subsequently, the discretion lies with the user to determine which of the daily BRDF descriptors are to be recorded as output data. It is noteworthy that in our analyses, we have opted to generate ReBeLS BRDF descriptors on a weekly schedule, aligning with the temporal frequency of the MCD43P dataset. Nevertheless, the temporal frequency for BRDF descriptor retrieval can be tailored to various intervals as per the user's requirements. For instance, it can be configured to operate on a daily basis, as exemplified in Fig. 7, or at dekadal intervals, in line with the Copernicus Global Land BRDF-corrected NDVI²³ and Copernicus Climate Change Surface Albedo²⁴ products.

### 4.1. BRDF descriptors retrievals

Fig. 2 displays an RGB image of the Oa17 BRDF descriptors corresponding to the 3rd June 2019. The colour channels are as follows: Red - $k_{iso}$, green - $k_{vol}$ and blue $k_{geo}$. The spatial consistency of the ReBeLS BRDF retrievals at the Oa17 band can be appraised from the spatial patterns coinciding with most spatial features associated with land cover which is consistent with the assumption that there is a clear connection between BRDF descriptors and land cover types belonging to the same categorical class. Any difference within the same land cover class will likely be about the magnitude of the BRDF rather than its shape (Strugnell et al., 2001). Nonetheless, based on the authors' experience, for some land cover types, variations in stem density can induce

<table><caption>**Table 2**

The number of observations used in the ReBeLS BRDF inversion and goodness of the model fits represented by the average ($< z_{sc} >$) and standard deviation ($\sigma_{z_{sc}}$) of the zeta-score for each of the stations in Table 1.</caption><thead><tr><th></th><th>Oa08<br>$N_{obs}$</th><th>Oa17<br>$N_{obs}$</th><th>Oa08<br>$< z_{sc} >, \sigma_{z_{sc}}$</th><th>Oa17<br>$< z_{sc} >, \sigma_{z_{sc}}$</th></tr></thead><tbody><tr><td>PT-Mi1</td><td>191</td><td>194</td><td>-0.03, 1.24</td><td>-0.07, 0.78</td></tr><tr><td>IT-Isp</td><td>170</td><td>170</td><td>-0.14, 1.49</td><td>0.11, 0.90</td></tr><tr><td>DE-Hai</td><td>102</td><td>120</td><td>0.26, 0.86</td><td>-0.07, 0.84</td></tr></tbody></table>

different BRDF shapes within the same land unit.

The ability of ReBeLS to reproduce the temporal evolution of the Sentinel-3 OLCI surface reflectance observations while maintaining a physically acceptable temporal evolution of BRDF-descriptors is examined for three European FLUXNET<fnref n="17" /> stations (see Table 1). The selected sites have latitudes ranging from 38.5° to 51.0° and exhibit temporal variability to test the temporal trajectories of the BRDF retrievals.

To estimate the goodness of the BRDF modelling, we compute the zeta-score ($z_{sc}$) as follows:

$$ z_{sc} = \frac{(\text{reflectance}_{\text{obs}} - \text{reflectance}_{\text{fwd}})}{\left( \sigma_{\text{reflectance}_{\text{obs}}}^2 + \sigma_{\text{reflectance}_{\text{fwd}}}^2 \right)^{1/2}} \quad (7) $$

where $\text{reflectance}_{\text{obs}}$ and $\text{reflectance}_{\text{fwd}}$ are the observed and model fitted (using retrieved BRDF model in forward mode) surface reflectances, and their associated uncertainties are denoted by $\sigma$. The zeta-score's absolute value indicates whether observed and model-fitted surface reflectances agree within uncertainties. A zeta-score with an absolute value above three should be regarded as a poor fit, while an absolute value of a zeta-score lower than two (three) shall be considered a good (an adequate) fit to the data. For each of the stations listed in Table 1, we report the number of observations used throughout the ReBeLS inversion, the average zeta-score value and its standard deviation for Oa08 and Oa17 spectral bands in Table 2.

The top-left panel of Fig. 3 shows the observed (open circles) and model-fitted (filled circles) Sentinel-3 surface reflectance time series in spectral band Oa17 for the pixel closest to the FLUXNET site PT-Mi1 (Portugal). The distribution of zeta-scores obtained for the Portuguese station is displayed in the bottom left panel of Fig. 3, where the colour bands are coded as follows: green denotes a good fit ($|z_{sc}| < 2$), yellow represents an adequate fit ($2 < |z_{sc}| < 3$), while points lying in the red band should be considered as a poor fit ($|z_{sc}| > 3$). The temporal trajectories of the retrieved Sentinel-3 Oa17 BRDF descriptors trajectories obtained by the ReBeLS processor are shown on the right panels (top: $k_{iso}$, middle: $k_{vol}$, bottom: $k_{geo}$) of Fig. 3. The MODIS Band2 MCD43A1 product for each Monday of the year considered is overplotted (red triangles in the right panels of Fig. 3) for comparison.

From the bottom left panel of Fig. 3, most of the model fits lie within the green zeta-score band, suggesting that our BRDF model successfully reproduces the Sentinel-3 observations of the PT-Mi1 station at spectral band Oa17. Fig. 4 displays the ReBeLS results for PT-Mi1 in the spectral band Oa08, and from the zeta-score panel, it can be seen that the large majority of the model fits are within the $\pm 2\sigma$ (green) band, although it is noticeable that few model fits are scattered across the yellow and red bands of the z-score panel which is reflected as an increase in the standard deviation of the zeta-score distribution as shown in Table 2.

The BRDF model fits and descriptor trajectories for IT-Isp and DE-Hai are shown in Figs. 5 and 6, respectively. From the top left panel of Fig. 6 and Table 2, it can be noticed that the availability of clear surface reflectance observations near the DE-Hai station is affected by a high

<footnote n="17">¹⁷ https://fluxnet.org/</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3. BRDF model retrieval results for station PT-Mi1 in Portugal.</figcaption></figure>

Top left panel: Observed and model fitted Sentinel-3 Oa17 surface reflectance time series. Bottom left panel: Zeta-score distribution over time; see text for description of colour bands. The temporal evolution of the BRDF descriptors for this station are shown in the right panels: isotropic (top), volumetric (middle) and geometric (bottom). BRDF descriptors retrieved from Sentinel-3 Oa17 data and MCD43A product are colour and symbol coded as shown in the legend.

<figure><img src="image_7.png" /><figcaption>Fig. 4. Same layout as in Fig. 3 but for spectral band Oa08.</figcaption></figure>

frequency of cloud occurrence compared to the Italian (IT-Isp, Fig. 5) and Portuguese (PT-Mi1, Fig. 3) stations. Nevertheless, the retrieved BRDF model for DE-Hai provides a good fit for the observations that is

consistent with the results for It-Isp and PT-Mi1 stations.

The goodness of fit parameters for the three selected stations, see Table 2, show that throughout the year selected, the time-evolving
<figure><img src="image_6.png" /><figcaption>Fig. 3. BRDF model retrieval results for station PT-Mi1 in Portugal.</figcaption></figure>

Top left panel: Observed and model fitted Sentinel-3 Oa17 surface reflectance time series. Bottom left panel: Zeta-score distribution over time; see text for description of colour bands. The temporal evolution of the BRDF descriptors for this station are shown in the right panels: isotropic (top), volumetric (middle) and geometric (bottom). BRDF descriptors retrieved from Sentinel-3 Oa17 data and MCD43A product are colour and symbol coded as shown in the legend.

<figure><img src="image_7.png" /><figcaption>Fig. 4. Same layout as in Fig. 3 but for spectral band Oa08.</figcaption></figure>

frequency of cloud occurrence compared to the Italian (IT-Isp, Fig. 5) and Portuguese (PT-Mi1, Fig. 3) stations. Nevertheless, the retrieved BRDF model for DE-Hai provides a good fit for the observations that is

consistent with the results for It-Isp and PT-Mi1 stations.

The goodness of fit parameters for the three selected stations, see Table 2, show that throughout the year selected, the time-evolving
<figure><img src="image_8.png" /><figcaption>Fig. 5. Same layout as in Fig. 3 but for station It-ISP, Italy.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Same layout as in Fig. 3 but for station De-Hai, Germany.</figcaption></figure>

ReBeLS BRDF models can reproduce their observations on average within ±0.14σ and ±0.08σ, for Oa08 and Oa17, respectively. The goodness of fit relation with wavelength is expected as the effects of the atmosphere are more significant in the Oa08 spectral band. Thus, in the presence of surface reflectance observations likely affected by the atmospheric correction procedure (or undetected outliers), the statistical

blend at the core of the ReBeLS algorithm will sacrifice goodness of fit and try to come up with a solution that does not depart significantly from the priors. The latter can be gleaned from Table 2, where it is noticeable that the standard deviations of the zeta-score in the Oa08 spectral band are higher than those in Oa17. Nevertheless, the average zeta-scores for both spectral bands are distributed around zero with
<figure><img src="image_8.png" /><figcaption>Fig. 5. Same layout as in Fig. 3 but for station It-ISP, Italy.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Same layout as in Fig. 3 but for station De-Hai, Germany.</figcaption></figure>

ReBeLS BRDF models can reproduce their observations on average within ±0.14σ and ±0.08σ, for Oa08 and Oa17, respectively. The goodness of fit relation with wavelength is expected as the effects of the atmosphere are more significant in the Oa08 spectral band. Thus, in the presence of surface reflectance observations likely affected by the atmospheric correction procedure (or undetected outliers), the statistical

blend at the core of the ReBeLS algorithm will sacrifice goodness of fit and try to come up with a solution that does not depart significantly from the priors. The latter can be gleaned from Table 2, where it is noticeable that the standard deviations of the zeta-score in the Oa08 spectral band are higher than those in Oa17. Nevertheless, the average zeta-scores for both spectral bands are distributed around zero with
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Daily ReBeLS BRDF adjusted surface reflectance for the Hainich station. The observations assembled during the extent of the time series and used in the BRDF inversion are also shown. For comparison, the NBAR time series from MODIS is overplotted. The datasets are symbol and colour coded, as shown in the legend. In the context of the MCD43 dataset, the visual representation includes all data points (–all: without regard to any quality flag) denoted by gray squares, and those characterised by a quality flag value of 0 (–qa = 0: good quality) are overplotted as purple squares. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Top panel: Original and predicted VIIRS VNP09GA surface reflectance in spectral channel I2. The datasets are symbol and colour coded, as shown in the legend. Bottom panel: Temporal evolution of the zeta-score computed from the original and predicted time series.</figcaption></figure>

standard deviations < 1.5σ meaning that on average the BRDF models approximated by ReBeLS can reproduce well the surface reflectance observations in both spectral bands within uncertainties.

In addition to the overall good fit of ReBeLS with the observations, the temporal trajectories of the BRDF descriptors retrieved are much less noisy than those of the MODIS MCD43A product, whereas their variations tend to resemble to variations that are typical of a vegetation growth (Jin and Eklundh, 2014). This result is consistent with previous findings (Quaife and Lewis, 2010) and can be taken as evidence of the advantage of using a regularisation approach throughout the BRDF inversion. However, within the temporal regularisation approach implemented in ReBeLS, the regularisation parameter has been assumed as constant. Although our cross-validation exercise ensures that the chosen regularisation parameter value creates a reasonable trade-off between smoothing and goodness of fit, the assumption of a generic regularisation parameter for all pixels in the global Sentinel-3 image may lead to the omission of genuine abrupt changes on the land surface. This can be considered as a caveat to the algorithm as discussed in section 4.3.

Since ReBeLS solves the BRDF inversion for each day comprised in the input reflectance time series, the weekly output cadence of the BRDF descriptors displayed in the above figures does not contribute to smoothing out the BRDF trajectories. To back up the latter argument, Fig. 7 shows the Sentinel-3 Oa17 surface reflectance observations used for the BRDF inversion (black circles) and the corresponding BRDF-adjusted surface reflectance (blue circles – see section 3.2) obtained from the daily retrieved BRDF descriptors, along with the nadir BRDF-adjusted time series from the MCD43A4 product (squares) for the pixel closest to the DE-Hai station. The daily ReBeLS BRDF-adjusted surface reflectance time series exhibit a phenological signal consistent with MCD43A4.

As seen in Fig. 7 during the autumn-winter period, no full-inversion was achievable in the MCD43A4 product, forcing the MCD43 algorithm to rely on the backup algorithm (shown as gray squares). On the other hand, ReBeLS was able to retrieve daily BRDF descriptors, allowing the computation of BRDF-adjusted surface reflectance, even in periods of persistent cloud coverage. For example, from November 2018 to February 2019, when the surface reflectance observations are very
<figure><img src="image_10.png" /><figcaption>Fig. 7. Daily ReBeLS BRDF adjusted surface reflectance for the Hainich station.</figcaption></figure>

The observations assembled during the extent of the time series and used in the BRDF inversion are also shown. For comparison, the NBAR time series from MODIS is overplotted. The datasets are symbol and colour coded, as shown in the legend. In the context of the MCD43 dataset, the visual representation includes all data points (–all: without regard to any quality flag) denoted by gray squares, and those characterised by a quality flag value of 0 (–qa = 0: good quality) are overplotted as purple squares. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

<figure><img src="image_11.png" /><figcaption>Fig. 8. Top panel: Original and predicted VIIRS VNP09GA surface reflectance in spectral channel I2.</figcaption></figure>

The datasets are symbol and colour coded, as shown in the legend. Bottom panel: Temporal evolution of the zeta-score computed from the original and predicted time series.

standard deviations < 1.5σ meaning that on average the BRDF models approximated by ReBeLS can reproduce well the surface reflectance observations in both spectral bands within uncertainties.

In addition to the overall good fit of ReBeLS with the observations, the temporal trajectories of the BRDF descriptors retrieved are much less noisy than those of the MODIS MCD43A product, whereas their variations tend to resemble to variations that are typical of a vegetation growth (Jin and Eklundh, 2014). This result is consistent with previous findings (Quaife and Lewis, 2010) and can be taken as evidence of the advantage of using a regularisation approach throughout the BRDF inversion. However, within the temporal regularisation approach implemented in ReBeLS, the regularisation parameter has been assumed as constant. Although our cross-validation exercise ensures that the chosen regularisation parameter value creates a reasonable trade-off between smoothing and goodness of fit, the assumption of a generic regularisation parameter for all pixels in the global Sentinel-3 image may lead to the omission of genuine abrupt changes on the land surface. This can be considered as a caveat to the algorithm as discussed in section 4.3.

Since ReBeLS solves the BRDF inversion for each day comprised in the input reflectance time series, the weekly output cadence of the BRDF descriptors displayed in the above figures does not contribute to smoothing out the BRDF trajectories. To back up the latter argument, Fig. 7 shows the Sentinel-3 Oa17 surface reflectance observations used for the BRDF inversion (black circles) and the corresponding BRDF-adjusted surface reflectance (blue circles – see section 3.2) obtained from the daily retrieved BRDF descriptors, along with the nadir BRDF-adjusted time series from the MCD43A4 product (squares) for the pixel closest to the DE-Hai station. The daily ReBeLS BRDF-adjusted surface reflectance time series exhibit a phenological signal consistent with MCD43A4.

As seen in Fig. 7 during the autumn-winter period, no full-inversion was achievable in the MCD43A4 product, forcing the MCD43 algorithm to rely on the backup algorithm (shown as gray squares). On the other hand, ReBeLS was able to retrieve daily BRDF descriptors, allowing the computation of BRDF-adjusted surface reflectance, even in periods of persistent cloud coverage. For example, from November 2018 to February 2019, when the surface reflectance observations are very
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** The zeta-score map was obtained by comparing the original and predicted (using the Sentinel-3 OLCI Oa17 BRDF model) VNP09GA observations. Each pixel represents the median z-score of a spatial extent of 3km × 3km. Levels of median zeta score are colour coded as displayed in the legend.</figcaption></figure>

sparse, the ReBeLS-retrieved BRDF parameters are close to the dynamic prior, and their uncertainties increase (reflecting a lack of knowledge of the surface) until new observations become available. A quenching of the uncertainties indicates that the BRDF model has been updated with observations, hence, becoming less dominated by the dynamic prior. The Bayesian optimisation framework used to solve the BRDF inversion problem in ReBeLS obviates the need for a backup algorithm, and increases the uncertainty in situations where the estimate cannot rely heavily on observations. Note that although historical prior is used, the smoothness constraint in ReBeLS ensures that the retrieved parameters are also consistent with other more data rich periods.

### 4.1.1. Capturing intrinsic land surface reflectance anisotropy

Since the model adopted by ReBeLS to approximate the BRDF of the land surface is semi-empirical, the BRDF retrievals from ReBeLS cannot be compared directly with in-situ measurements nor to BRDF descriptors retrievals obtained with another sensor. Nevertheless, if the BRDF model derived from Sentinel-3 properly describes a target's intrinsic surface reflectance anisotropy, it should mimic its directional signature observed with another sensor. For such, we selected VIIRS as a reference sensor because its spatial resolution is close to OLCI.

Fig. 8 displays the temporal evolution of the VNP09GA band I2 (0.85 – 0.88 µm) surface reflectance for a pixel near the PT-Mi1 station. Since no surface reflectance measurement uncertainties are provided in the VNP09GA products, we have taken the conservative approach to assign a surface reflectance uncertainty of 5% to VNP09GA surface reflectance measurements. The criterion to perform the anisotropy evaluation for spectral bands Oa17 (Sentinel-3 OLCI) and I2 (VIIRS) is supported by the fact that surface reflectances in the NIR spectral range are less impacted by the atmospheric correction. Therefore, we estimate predictions of VIIRS observed reflectances in the spectral band I2 by running the Oa17 Sentinel-3 BRDF model in forward mode, assuming the same acquisition geometry as VIIRS. In other words, we compute the RossThickLiSparseReciprocal angular functions using the VIIRS geometry of acquisition and, together with the Sentinel-3 retrieved BRDF

descriptors, we compute the predicted surface reflectance by using Eq. (1). These predictions are shown in the top panel of Fig. 8 as blue squares. Their associated uncertainties have been estimated by propagating the uncertainties of the retrieved BRDF descriptors.

A good agreement between the original and the predicted time series can be gleaned at first glance from the top panel of Fig. 8. Moreover, we compute the zeta score (see Eq. (7)) to assess the agreement among time series quantitatively. The bottom panel of Fig. 8 displays the temporal evolution of the zeta-score being contained within the ±2σ band yielding a median zeta-score $\mu_{z_c} = -0.10$. This statistical results allow us to conclude that since the differences between the datasets fall within uncertainties, the ReBeLS Sentinel-3 BRDF model can depict properly the surface reflectance temporal variations associated with the target's phenological changes and its intrinsic surface reflectance anisotropy.

To investigate whether the level of agreement between VIIRS VNP09GA and ReBeLS BRDF predicted surface reflectance time series could be achieved at a larger spatial scale, we have assembled a dataset at the European scale containing VNP09GA surface reflectances in VIIRS band I2 and ReBeLS Sentinel-3 Oa17 BRDF model to apply the recipe outlined above. After computing the zeta-score for all pairs of original and predicted measurements for the period considered in this article, we compute the mean zeta score for a macro pixel of spatial extent 3km × 3km, thereby facilitating a comprehensive assessment of the overall agreement between datasets. This median zeta-score image degraded from 333 m to 3 km pixel size is shown in Fig. 9. Note that the full range of zeta-scores displayed in Fig. 9 are well contained within a ±1 σ, albeit for some coastal regions or parts of the tile encompassing ocean where zeta-scores are beyond 1 (dark red regions).

It is worth outlining that none of the VIIRS VNP09GA surface reflectance observations have been used in the Sentinel-3 BRDF model retrieval, neither spectral and spatial adjustments between VIIRS and Sentinel-3 were achieved. Therefore, the good agreement between observed and predicted VIIRS VNP09GA surface reflectances supports evidence that the ReBeLS BRDF can successfully reproduce the intrinsic surface reflectance anisotropy of any land surface.
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** The zeta-score map was obtained by comparing the original and predicted (using the Sentinel-3 OLCI Oa17 BRDF model) VNP09GA observations. Each pixel represents the median z-score of a spatial extent of 3km × 3km. Levels of median zeta score are colour coded as displayed in the legend.</figcaption></figure>

sparse, the ReBeLS-retrieved BRDF parameters are close to the dynamic prior, and their uncertainties increase (reflecting a lack of knowledge of the surface) until new observations become available. A quenching of the uncertainties indicates that the BRDF model has been updated with observations, hence, becoming less dominated by the dynamic prior. The Bayesian optimisation framework used to solve the BRDF inversion problem in ReBeLS obviates the need for a backup algorithm, and increases the uncertainty in situations where the estimate cannot rely heavily on observations. Note that although historical prior is used, the smoothness constraint in ReBeLS ensures that the retrieved parameters are also consistent with other more data rich periods.

### 4.1.1. Capturing intrinsic land surface reflectance anisotropy

Since the model adopted by ReBeLS to approximate the BRDF of the land surface is semi-empirical, the BRDF retrievals from ReBeLS cannot be compared directly with in-situ measurements nor to BRDF descriptors retrievals obtained with another sensor. Nevertheless, if the BRDF model derived from Sentinel-3 properly describes a target's intrinsic surface reflectance anisotropy, it should mimic its directional signature observed with another sensor. For such, we selected VIIRS as a reference sensor because its spatial resolution is close to OLCI.

Fig. 8 displays the temporal evolution of the VNP09GA band I2 (0.85 – 0.88 µm) surface reflectance for a pixel near the PT-Mi1 station. Since no surface reflectance measurement uncertainties are provided in the VNP09GA products, we have taken the conservative approach to assign a surface reflectance uncertainty of 5% to VNP09GA surface reflectance measurements. The criterion to perform the anisotropy evaluation for spectral bands Oa17 (Sentinel-3 OLCI) and I2 (VIIRS) is supported by the fact that surface reflectances in the NIR spectral range are less impacted by the atmospheric correction. Therefore, we estimate predictions of VIIRS observed reflectances in the spectral band I2 by running the Oa17 Sentinel-3 BRDF model in forward mode, assuming the same acquisition geometry as VIIRS. In other words, we compute the RossThickLiSparseReciprocal angular functions using the VIIRS geometry of acquisition and, together with the Sentinel-3 retrieved BRDF

descriptors, we compute the predicted surface reflectance by using Eq. (1). These predictions are shown in the top panel of Fig. 8 as blue squares. Their associated uncertainties have been estimated by propagating the uncertainties of the retrieved BRDF descriptors.

A good agreement between the original and the predicted time series can be gleaned at first glance from the top panel of Fig. 8. Moreover, we compute the zeta score (see Eq. (7)) to assess the agreement among time series quantitatively. The bottom panel of Fig. 8 displays the temporal evolution of the zeta-score being contained within the ±2σ band yielding a median zeta-score $\mu_{z_c} = -0.10$. This statistical results allow us to conclude that since the differences between the datasets fall within uncertainties, the ReBeLS Sentinel-3 BRDF model can depict properly the surface reflectance temporal variations associated with the target's phenological changes and its intrinsic surface reflectance anisotropy.

To investigate whether the level of agreement between VIIRS VNP09GA and ReBeLS BRDF predicted surface reflectance time series could be achieved at a larger spatial scale, we have assembled a dataset at the European scale containing VNP09GA surface reflectances in VIIRS band I2 and ReBeLS Sentinel-3 Oa17 BRDF model to apply the recipe outlined above. After computing the zeta-score for all pairs of original and predicted measurements for the period considered in this article, we compute the mean zeta score for a macro pixel of spatial extent 3km × 3km, thereby facilitating a comprehensive assessment of the overall agreement between datasets. This median zeta-score image degraded from 333 m to 3 km pixel size is shown in Fig. 9. Note that the full range of zeta-scores displayed in Fig. 9 are well contained within a ±1 σ, albeit for some coastal regions or parts of the tile encompassing ocean where zeta-scores are beyond 1 (dark red regions).

It is worth outlining that none of the VIIRS VNP09GA surface reflectance observations have been used in the Sentinel-3 BRDF model retrieval, neither spectral and spatial adjustments between VIIRS and Sentinel-3 were achieved. Therefore, the good agreement between observed and predicted VIIRS VNP09GA surface reflectances supports evidence that the ReBeLS BRDF can successfully reproduce the intrinsic surface reflectance anisotropy of any land surface.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** BRDF-adjusted surface reflectance temporal evolution for a pixel near Gandia, Spain. This site was affected by a wildfire in the summer of 2018. The ReBeLS BRDF-adjusted surface reflectance (Oa17), NBAR MODIS (Band2) and MCD43P predicted surface reflectance time series are colour and symbol encoded as shown in the legend.</figcaption></figure>

Since there appears to exist such a close similarity in terms of performance and biases between MODIS and VIIRS BRDF products (Liu et al., 2017), one can argue that the good agreement between VNP09GA and the predicted surface reflectance from ReBeLS BRDF descriptors is driven by the MCD43P prior used in ReBeLS. However, the ReBeLS way of solving the BRDF inversion should be understood as a balanced statistical blend where the availability of surface reflectance observations and their associated uncertainties determine when and how much the algorithm trusts the prior and observations, respectively. For example, in Fig. 10, we show the surface reflectance time series for a pixel near Gandia in Spain. This region was affected by a wildfire in the summer of 2018. The predicted surface reflectance from the MCD43P prior dataset is shown in yellow squares. Since the MCD43P dataset is a climatology built from six years of data, it represents the average anisotropy of the surface. It does not reproduce the sharp decline in surface reflectance around mid-August, associated with the fire seen in the MCD43A NBAR and ReBeLS BRDF-adjusted surface reflectance time series as shown in Fig. 10. The ReBeLS algorithm captured the rapid change on the surface timely, although producing a slight smoothing of the sharp drop in reflectance caused by the fire.

## 4.2. BRDF-adjusted surface reflectance

### 4.2.1. Mosaicking

The Maximum Value Composite [MVC] (Holben, 1986) is a technique that allows the mosaicking of surface reflectance images. There are several major reasons to rethink MVC (Viovy et al., 2007; Zeng et al., 2020), as it is still in use (Dierckx et al., 2014). The bottom panel of Fig. 11 shows the MVC synthesis of Sentinel-3 atmospherically corrected surface reflectance for a geographic region comprising Europe. Since the MVC procedure retains clear pixels with the maximum NDVI values, MVC daily composites combine the different Sentinel-3 A and Sentinel-3B overpasses during the day. The combination of Sentinel-3 overpasses with different Sun-sensor acquisition geometry becomes conspicuous as sharp changes in daily MVC surface reflectance, as seen from the bottom panel of Fig. 11. If a BRDF adjustment (BRDF-correction) successfully quenches the directional dependence of the land surface reflectance, then the mosaicking pattern (i.e. sharp discontinuities across the image) introduced by MVC should not be present in a BRDF-adjusted mosaic. As shown in the top panel of Fig. 11, no discontinuity in the ReBeLS BRDF-adjusted surface reflectance image is seen across the orbital segment direction.

### 4.2.2. Comparison with VIIRS VNP43IA4

Fig. 12 shows the temporal evolution of surface reflectance acquired by the Sentinel-3 constellation in spectral band Oa17 for a pixel near Grosseto, Italy, where cultivated areas dominate the land cover. The distinctive large variations shown by the surface reflectance between

consecutive acquisitions are to do with changes in the view/illumination geometries, and the underlying anisotropy of the surface (Bréon and Vermote, 2012; Bréon et al., 2015; Roujean et al., 1992b; Vermote et al., 2009). In the same figure (Fig. 12), we have plotted the ReBeLS normalised nadir BRDF for band Oa17, as well as the equivalent estimate from the VNP43A4 product for band I2. We can see a striking similarity between the two normalised products, even ignoring the difference in spectral, spatial and temporal integration used in the different products.

To investigate whether the level of agreement between Sentinel-3 and VIIRS BRDF-adjusted surface reflectance holds at a larger scale, in the top left panel of Fig. 13, we show the VIIRS VNP43IA4 mosaic in spectral band I2 for a day during the winter of 2019 (day of year 43). The right panel of Fig. 13 shows the ReBeLS Sentinel-3 BRDF-adjusted surface reflectance in the Oa17 band for the same geographical region as the left panel. Following the quality layer information of the VIIRS VNP43IA4 product, we have retained only good-quality pixels as described in section 2.2. A visual inspection of the top panels in Fig. 13 reveals that most Scandinavia and North-East Europe pixels are flagged out in the VIIRS (left panel) and Sentinel-3 (right panel) BRDF-adjusted reflectance products. The lack of clear observations due to persistent cloud coverage combined with periods of natural darkness at northern latitudes during the winter has noticeably impacted the completeness of both sensors' BRDF-adjusted surface reflectances products.

Despite the spectral mismatch between Sentinel-3 Oa17 and VIIRS I2 spectral bands, our goal is to compare both BRDF-adjusted reflectance levels directly. For this purpose, the bottom panel in Fig. 13 shows the relationship between BRDF-adjusted surface reflectance for all pixels in the mosaics and the best-fitted bisector line modelling the VIIRS NBAR–Sentinel-3 BAR relation is shown in the bottom panel of Fig. 13 as a dashed green line. The corresponding slope and offset are $\alpha = 0.957 \pm 2.307 \times 10^{-4}$ and $\beta = 0.016 \pm 9.488 \times 10^{-5}$, respectively. The non-parametric Spearman ($r$) and Kendall ($\tau$) correlation coefficients of 0.95 and 0.83, respectively, suggest a statistically significant ($p \ll 0.05$) tight relation between VIIRS and Sentinel-3 NBAR. Note that the best-fitted relation is close to the identity relation represented in the bottom panel of Fig. 13 with a solid white line.

The same exercise is repeated for mosaics of VIIRS NBAR and Sentinel-3 BRDF-adjusted products for a day during the summer. Following the same layout as in Fig. 13, the results for day 154 of 2019 are shown in Fig. 14. We notice that the number of good-quality pixels in both BRDF-adjusted surface reflectance products is significantly higher for the summer mosaic than for the winter one. Such an increase in the population of good-quality pixels is expected due to the reduction of cloud coverage, being most noticeable in the northern regions. The modelled linear relation between VIIRS NBAR and Sentinel-3 TOCR mosaics in the NIR spectral range for the selected day in the summer shows a relation ($\alpha = 1.017 \pm 2.94 \times 10^{-4}$, $\beta = -0.005 \pm 1.257 \times 10^{-4}$) that is fairly tight ($r = 0.93$, $\tau = 0.79$ with $p \ll 0.05$) close to the
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** BRDF-adjusted surface reflectance temporal evolution for a pixel near Gandia, Spain. This site was affected by a wildfire in the summer of 2018. The ReBeLS BRDF-adjusted surface reflectance (Oa17), NBAR MODIS (Band2) and MCD43P predicted surface reflectance time series are colour and symbol encoded as shown in the legend.</figcaption></figure>

Since there appears to exist such a close similarity in terms of performance and biases between MODIS and VIIRS BRDF products (Liu et al., 2017), one can argue that the good agreement between VNP09GA and the predicted surface reflectance from ReBeLS BRDF descriptors is driven by the MCD43P prior used in ReBeLS. However, the ReBeLS way of solving the BRDF inversion should be understood as a balanced statistical blend where the availability of surface reflectance observations and their associated uncertainties determine when and how much the algorithm trusts the prior and observations, respectively. For example, in Fig. 10, we show the surface reflectance time series for a pixel near Gandia in Spain. This region was affected by a wildfire in the summer of 2018. The predicted surface reflectance from the MCD43P prior dataset is shown in yellow squares. Since the MCD43P dataset is a climatology built from six years of data, it represents the average anisotropy of the surface. It does not reproduce the sharp decline in surface reflectance around mid-August, associated with the fire seen in the MCD43A NBAR and ReBeLS BRDF-adjusted surface reflectance time series as shown in Fig. 10. The ReBeLS algorithm captured the rapid change on the surface timely, although producing a slight smoothing of the sharp drop in reflectance caused by the fire.

## 4.2. BRDF-adjusted surface reflectance

### 4.2.1. Mosaicking

The Maximum Value Composite [MVC] (Holben, 1986) is a technique that allows the mosaicking of surface reflectance images. There are several major reasons to rethink MVC (Viovy et al., 2007; Zeng et al., 2020), as it is still in use (Dierckx et al., 2014). The bottom panel of Fig. 11 shows the MVC synthesis of Sentinel-3 atmospherically corrected surface reflectance for a geographic region comprising Europe. Since the MVC procedure retains clear pixels with the maximum NDVI values, MVC daily composites combine the different Sentinel-3 A and Sentinel-3B overpasses during the day. The combination of Sentinel-3 overpasses with different Sun-sensor acquisition geometry becomes conspicuous as sharp changes in daily MVC surface reflectance, as seen from the bottom panel of Fig. 11. If a BRDF adjustment (BRDF-correction) successfully quenches the directional dependence of the land surface reflectance, then the mosaicking pattern (i.e. sharp discontinuities across the image) introduced by MVC should not be present in a BRDF-adjusted mosaic. As shown in the top panel of Fig. 11, no discontinuity in the ReBeLS BRDF-adjusted surface reflectance image is seen across the orbital segment direction.

### 4.2.2. Comparison with VIIRS VNP43IA4

Fig. 12 shows the temporal evolution of surface reflectance acquired by the Sentinel-3 constellation in spectral band Oa17 for a pixel near Grosseto, Italy, where cultivated areas dominate the land cover. The distinctive large variations shown by the surface reflectance between

consecutive acquisitions are to do with changes in the view/illumination geometries, and the underlying anisotropy of the surface (Bréon and Vermote, 2012; Bréon et al., 2015; Roujean et al., 1992b; Vermote et al., 2009). In the same figure (Fig. 12), we have plotted the ReBeLS normalised nadir BRDF for band Oa17, as well as the equivalent estimate from the VNP43A4 product for band I2. We can see a striking similarity between the two normalised products, even ignoring the difference in spectral, spatial and temporal integration used in the different products.

To investigate whether the level of agreement between Sentinel-3 and VIIRS BRDF-adjusted surface reflectance holds at a larger scale, in the top left panel of Fig. 13, we show the VIIRS VNP43IA4 mosaic in spectral band I2 for a day during the winter of 2019 (day of year 43). The right panel of Fig. 13 shows the ReBeLS Sentinel-3 BRDF-adjusted surface reflectance in the Oa17 band for the same geographical region as the left panel. Following the quality layer information of the VIIRS VNP43IA4 product, we have retained only good-quality pixels as described in section 2.2. A visual inspection of the top panels in Fig. 13 reveals that most Scandinavia and North-East Europe pixels are flagged out in the VIIRS (left panel) and Sentinel-3 (right panel) BRDF-adjusted reflectance products. The lack of clear observations due to persistent cloud coverage combined with periods of natural darkness at northern latitudes during the winter has noticeably impacted the completeness of both sensors' BRDF-adjusted surface reflectances products.

Despite the spectral mismatch between Sentinel-3 Oa17 and VIIRS I2 spectral bands, our goal is to compare both BRDF-adjusted reflectance levels directly. For this purpose, the bottom panel in Fig. 13 shows the relationship between BRDF-adjusted surface reflectance for all pixels in the mosaics and the best-fitted bisector line modelling the VIIRS NBAR–Sentinel-3 BAR relation is shown in the bottom panel of Fig. 13 as a dashed green line. The corresponding slope and offset are $\alpha = 0.957 \pm 2.307 \times 10^{-4}$ and $\beta = 0.016 \pm 9.488 \times 10^{-5}$, respectively. The non-parametric Spearman ($r$) and Kendall ($\tau$) correlation coefficients of 0.95 and 0.83, respectively, suggest a statistically significant ($p \ll 0.05$) tight relation between VIIRS and Sentinel-3 NBAR. Note that the best-fitted relation is close to the identity relation represented in the bottom panel of Fig. 13 with a solid white line.

The same exercise is repeated for mosaics of VIIRS NBAR and Sentinel-3 BRDF-adjusted products for a day during the summer. Following the same layout as in Fig. 13, the results for day 154 of 2019 are shown in Fig. 14. We notice that the number of good-quality pixels in both BRDF-adjusted surface reflectance products is significantly higher for the summer mosaic than for the winter one. Such an increase in the population of good-quality pixels is expected due to the reduction of cloud coverage, being most noticeable in the northern regions. The modelled linear relation between VIIRS NBAR and Sentinel-3 TOCR mosaics in the NIR spectral range for the selected day in the summer shows a relation ($\alpha = 1.017 \pm 2.94 \times 10^{-4}$, $\beta = -0.005 \pm 1.257 \times 10^{-4}$) that is fairly tight ($r = 0.93$, $\tau = 0.79$ with $p \ll 0.05$) close to the
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Top panel: ReBeLS Sentinel-3 BRDF-adjusted surface reflectance for OLCI band Oa17. Bottom panel: MVC mosaic of Sentinel-3 OLCI Oa17 surface reflectance for the 3rd of June 2019. MVC mosaicking patterns, introduced by merging orbital segments with different acquisition geometries, are evident.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Top panel: ReBeLS Sentinel-3 BRDF-adjusted surface reflectance for OLCI band Oa17. Bottom panel: MVC mosaic of Sentinel-3 OLCI Oa17 surface reflectance for the 3rd of June 2019. MVC mosaicking patterns, introduced by merging orbital segments with different acquisition geometries, are evident.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Observed (black circles) and BRDF-adjusted (blue circles) surface reflectance in the Sentinel-3 Oa17 spectral band for a pixel near Grosseto.</figcaption></figure>

The ReBeLS BRDF-adjusted surface reflectance is overplotted with blue circles. The nadir BRDF-adjusted surface reflectance from the VIIRS product VNP43IA4 in spectral channel I2 is shown as red triangles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

identity line. Given that the computation of NBAR is obtained via the retrieved BRDF descriptors and assuming a normalisation acquisition geometry – see section 3.2, the statistically significant correlation between VIIRS and Sentinel-3 NBAR can be taken as evidence that both BRDF models, despite being retrieved via distinct sensors and algorithmic approaches, are able to reproduce the land surface intrinsic directional signature. Although a slight tendency towards overestimation of NBAR for ReBeLS can be noticed in cases of high surface reflectances, this could be associated with differences in spectral bands and atmospheric processing methodologies across datasets. Therefore, the tests described above give us additional confidence in the reliability of our BRDF-descriptors retrieval and BRDF-adjustment algorithms implemented in the ReBeLS processor to normalise Sentinel-3 surface reflectance data into a common Sun-sensor geometry.

### 4.3. Caveats

The prior of the BRDF descriptors is central to the ReBeLS methodology for retrieving Sentinel-3 BRDF descriptors. Hence it is worth noticing that the resampling (from 500 m to 333 m) and reprojection (from sinusoidal to Plate-Carré projection) of the MCD43A1 and MCD43A2 products to match the Sentinel-3 OLCI data grid may have resulted in some intrinsic modification. However, a visual examination of the MCD43P dataset did not reveal any artefact associated with the MCD43P preprocessing stages. The MCD43P dataset has a weekly cadence, albeit the suite of MCD43A products is generated for a daily time grid. Nevertheless, storing daily MCD43A1 and MCD43A2 global products over six years would rapidly result in a sizeable storage capacity that was not accessible when designing the ReBeLS processor. Thus, a weekly cadence for the MCD43P dataset emerged as a reasonable trade-off between the high temporal resolution of the MCD43A products and the dekadal (10 days) frequency requested for the operational products derived from ReBeLS (i.e. Copernicus Global Land Service's BRDF corrected NDVI and Copernicus Climate Change Service's surface albedo (Sánchez-Zapero et al., 2023)). The adoption of Mondays follows the approach used in the operational University of Natural Resources and Applied Life Sciences (BOKU) NDVI product – based on MODIS (Klisch and Atzberger, 2016) and successfully exploited for operational timely critical vegetation assessment purposes (Meroni et al., 2019).

It should also be noted that while assembling the MCD43P dataset, no attempt was made to address the spectral harmonisation between Sentinel-3 OLCI and MODIS. The spectral band mismatch between these sensors introduces uncertainty into the BRDF descriptor retrieval procedure, albeit to a lesser extent than the uncertainties associated with the MCD43P prior kernel weights. Therefore, a spectral band adjustment is unlikely to have a great effect on the current BRDF retrievals or their uncertainty. Nevertheless, a spectral band mismatch exercise might be

instructive, and future developments of the ReBeLS processor should consider such spectral harmonisation between Sentinel-3 and MODIS, VIIRS, etc.

Although the MCD43A2 product offers information on the snow pixel classification, we have not considered it while assembling the MCD43P dataset. The latter is driven by the presumption that accounting for the mean and variance of a dataset made up of pixels with and without snow must account for the significant variation in kernel weights. Therefore, a large variance associated with the prior BRDF descriptors can drastically loosen the inversion constraints, enabling the temporal development of the BRDF descriptors to change according to the true surface state. Unfortunately, due to the significant limitations on the cloud mask identification<fnref n="18" /> of the Sentinel-3 atmospherically corrected surface reflectance input data, the above assumption has not been confirmed in the present version of the ReBeLS processor.

The MCD43 suite of products version was V006 when the MCD43P dataset was assembled. However, an updated version, V006.1, has been made available since 2021, and any future updates of the ReBeLS processor must update the MCD43P dataset with the most recent version of MCD43A products or using a BRDF prior derived from Sentinel-3.

It is important to acknowledge that the smoothness in ReBeLS BRDF descriptors and BRDF-adjusted surface reflectance trajectories is achieved by imposing a first-order difference to the cost function solved in the BRDF inversion problem. Further, we assume a constant regularisation parameter (which can be understood as how much we trust the first-order difference) across time and space. As such, the results of the cross-validation exercise described in Appendix A are pivotal to achieving the appropriate balance between smoothing and accuracy, allowing us to ensure that the chosen regularisation parameter value creates a reasonable trade-off between smoothing and goodness of fit. Nonetheless, treating the regularisation parameter as a constant may lead to the omission of genuine high-frequency variations in the temporal trajectories of the BRDF descriptors, a negligible issue if these variations are covered by the uncertainty envelope. On the other hand, significant and rapid changes on the land surface, such as a transitory snowfall – snowmelt event, can potentially be missed by the algorithm. Conversely, abrupt non-transient changes like those associated with the remote sensing signal of a large wild fire will be detected by ReBeLS as shown in Fig. 10, although a small degree of temporal blurring and increased uncertainty may occur.

Thus, while the regularisation approach implemented in ReBeLS can have an impact on the shape of the BRDF-descriptor trajectories, the

<footnote n="18">¹⁸ Snow pixels are mistakenly identified as clouds even in geographic areas where snow is a constant feature in the landscape. For details see https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_QAR_S3-CloudMask_II.10.pdf.</footnote>
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Observed (black circles) and BRDF-adjusted (blue circles) surface reflectance in the Sentinel-3 Oa17 spectral band for a pixel near Grosseto. The ReBeLS BRDF-adjusted surface reflectance is overplotted with blue circles. The nadir BRDF-adjusted surface reflectance from the VIIRS product VNP43IA4 in spectral channel I2 is shown as red triangles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

identity line. Given that the computation of NBAR is obtained via the retrieved BRDF descriptors and assuming a normalisation acquisition geometry – see section 3.2, the statistically significant correlation between VIIRS and Sentinel-3 NBAR can be taken as evidence that both BRDF models, despite being retrieved via distinct sensors and algorithmic approaches, are able to reproduce the land surface intrinsic directional signature. Although a slight tendency towards overestimation of NBAR for ReBeLS can be noticed in cases of high surface reflectances, this could be associated with differences in spectral bands and atmospheric processing methodologies across datasets. Therefore, the tests described above give us additional confidence in the reliability of our BRDF-descriptors retrieval and BRDF-adjustment algorithms implemented in the ReBeLS processor to normalise Sentinel-3 surface reflectance data into a common Sun-sensor geometry.

### 4.3. Caveats

The prior of the BRDF descriptors is central to the ReBeLS methodology for retrieving Sentinel-3 BRDF descriptors. Hence it is worth noticing that the resampling (from 500 m to 333 m) and reprojection (from sinusoidal to Plate-Carré projection) of the MCD43A1 and MCD43A2 products to match the Sentinel-3 OLCI data grid may have resulted in some intrinsic modification. However, a visual examination of the MCD43P dataset did not reveal any artefact associated with the MCD43P preprocessing stages. The MCD43P dataset has a weekly cadence, albeit the suite of MCD43A products is generated for a daily time grid. Nevertheless, storing daily MCD43A1 and MCD43A2 global products over six years would rapidly result in a sizeable storage capacity that was not accessible when designing the ReBeLS processor. Thus, a weekly cadence for the MCD43P dataset emerged as a reasonable trade-off between the high temporal resolution of the MCD43A products and the dekadal (10 days) frequency requested for the operational products derived from ReBeLS (i.e. Copernicus Global Land Service's BRDF corrected NDVI and Copernicus Climate Change Service's surface albedo (Sánchez-Zapero et al., 2023)). The adoption of Mondays follows the approach used in the operational University of Natural Resources and Applied Life Sciences (BOKU) NDVI product – based on MODIS (Klisch and Atzberger, 2016) and successfully exploited for operational timely critical vegetation assessment purposes (Meroni et al., 2019).

It should also be noted that while assembling the MCD43P dataset, no attempt was made to address the spectral harmonisation between Sentinel-3 OLCI and MODIS. The spectral band mismatch between these sensors introduces uncertainty into the BRDF descriptor retrieval procedure, albeit to a lesser extent than the uncertainties associated with the MCD43P prior kernel weights. Therefore, a spectral band adjustment is unlikely to have a great effect on the current BRDF retrievals or their uncertainty. Nevertheless, a spectral band mismatch exercise might be

instructive, and future developments of the ReBeLS processor should consider such spectral harmonisation between Sentinel-3 and MODIS, VIIRS, etc.

Although the MCD43A2 product offers information on the snow pixel classification, we have not considered it while assembling the MCD43P dataset. The latter is driven by the presumption that accounting for the mean and variance of a dataset made up of pixels with and without snow must account for the significant variation in kernel weights. Therefore, a large variance associated with the prior BRDF descriptors can drastically loosen the inversion constraints, enabling the temporal development of the BRDF descriptors to change according to the true surface state. Unfortunately, due to the significant limitations on the cloud mask identification<fnref n="18" /> of the Sentinel-3 atmospherically corrected surface reflectance input data, the above assumption has not been confirmed in the present version of the ReBeLS processor.

The MCD43 suite of products version was V006 when the MCD43P dataset was assembled. However, an updated version, V006.1, has been made available since 2021, and any future updates of the ReBeLS processor must update the MCD43P dataset with the most recent version of MCD43A products or using a BRDF prior derived from Sentinel-3.

It is important to acknowledge that the smoothness in ReBeLS BRDF descriptors and BRDF-adjusted surface reflectance trajectories is achieved by imposing a first-order difference to the cost function solved in the BRDF inversion problem. Further, we assume a constant regularisation parameter (which can be understood as how much we trust the first-order difference) across time and space. As such, the results of the cross-validation exercise described in Appendix A are pivotal to achieving the appropriate balance between smoothing and accuracy, allowing us to ensure that the chosen regularisation parameter value creates a reasonable trade-off between smoothing and goodness of fit. Nonetheless, treating the regularisation parameter as a constant may lead to the omission of genuine high-frequency variations in the temporal trajectories of the BRDF descriptors, a negligible issue if these variations are covered by the uncertainty envelope. On the other hand, significant and rapid changes on the land surface, such as a transitory snowfall – snowmelt event, can potentially be missed by the algorithm. Conversely, abrupt non-transient changes like those associated with the remote sensing signal of a large wild fire will be detected by ReBeLS as shown in Fig. 10, although a small degree of temporal blurring and increased uncertainty may occur.

Thus, while the regularisation approach implemented in ReBeLS can have an impact on the shape of the BRDF-descriptor trajectories, the

<footnote n="18">¹⁸ Snow pixels are mistakenly identified as clouds even in geographic areas where snow is a constant feature in the landscape. For details see https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_QAR_S3-CloudMask_II.10.pdf.</footnote>
year: 2019; doy: 042

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Left top panel: VIIRS VNP43IA4 mosaic in spectral band I2 for day 43 of 2019 (winter). Right top panel: Sentinel-3 ReBeLS BRDF-adjusted mosaic in Oa17 spectral band for the same day. Bottom panel: Scatter plot between VNP43IA4 and ReBeLS BRDF-adjusted surface reflectance. The identity line is shown as a solid white line, while the best linear model fit is displayed with a green dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

cross-validation exercises, together with the dynamic prior, ensure that the resulting BRDF-descriptors trajectories are physically plausible and can predict the Sentinel-3 reflectance observations within their uncertainty, as shown in the bottom-left panels of Figs. 4 to 7. Moreover, dedicated investigations evaluating the ReBeLS-derived products have demonstrated favourable temporal consistency with datasets obtained from other satellites and in-situ data. For instance, in the study conducted by Sánchez-Zapero et al. (2023), the Sentinel-3 broad-band surface albedo product, computed from the ReBeLS BRDF descriptors, exhibited spatial and temporal consistency comparable to surface albedo

products derived from MODIS and PROBA-V. Direct comparison with broad-band surface albedo in-situ data showed an accuracy of 0.005, a precision of 0.016, and an uncertainty of 0.032 (Sánchez-Zapero et al., 2023), in line with in-situ validation metrics for the MODIS surface albedo product MCD43 A3 C6.

The ReBeLS algorithm offers an option to split the region of interest into small chunks to be processed in parallel within a Hadoop cluster.<fnref n="19" />

<footnote n="19">¹⁹ https://hadoop.apache.org/</footnote>
year: 2019; doy: 042

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Left top panel: VIIRS VNP43IA4 mosaic in spectral band I2 for day 43 of 2019 (winter). Right top panel: Sentinel-3 ReBeLS BRDF-adjusted mosaic in Oa17 spectral band for the same day. Bottom panel: Scatter plot between VNP43IA4 and ReBeLS BRDF-adjusted surface reflectance. The identity line is shown as a solid white line, while the best linear model fit is displayed with a green dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

cross-validation exercises, together with the dynamic prior, ensure that the resulting BRDF-descriptors trajectories are physically plausible and can predict the Sentinel-3 reflectance observations within their uncertainty, as shown in the bottom-left panels of Figs. 4 to 7. Moreover, dedicated investigations evaluating the ReBeLS-derived products have demonstrated favourable temporal consistency with datasets obtained from other satellites and in-situ data. For instance, in the study conducted by Sánchez-Zapero et al. (2023), the Sentinel-3 broad-band surface albedo product, computed from the ReBeLS BRDF descriptors, exhibited spatial and temporal consistency comparable to surface albedo

products derived from MODIS and PROBA-V. Direct comparison with broad-band surface albedo in-situ data showed an accuracy of 0.005, a precision of 0.016, and an uncertainty of 0.032 (Sánchez-Zapero et al., 2023), in line with in-situ validation metrics for the MODIS surface albedo product MCD43 A3 C6.

The ReBeLS algorithm offers an option to split the region of interest into small chunks to be processed in parallel within a Hadoop cluster.<fnref n="19" />

<footnote n="19">¹⁹ https://hadoop.apache.org/</footnote>
<figure><img src="image_17.png" /><figcaption>Fig. 14. The layout is the same as Fig. 13 but for day 154 of 2019 (summer).</figcaption></figure>

The parallelised mode of ReBeLS has been successfully integrated into the operational processing chain of the Copernicus Global Land and Climate Change Services, which generates operational Sentinel-3-based BRDF-adjusted surface reflectances – further used to compute an NDVI product – and surface albedo products, respectively (León-Tavares, 2022; Sánchez-Zapero et al., 2023).

At the time of writing this article, ReBeLS remains confined to proprietary access, thereby granting the scientific community access to ReBeLS-derived products encompassing NDVI and surface albedo via the Copernicus Global Land and Copernicus Climate change services, respectively. It is important to acknowledge that the quality of ReBeLS outcomes is intimately linked to the quality and temporal extent of the surface reflectance time series Nonetheless, we acknowledge that the computational processing of long-term time series, such as those spanning multiple years, could potentially introduce technical challenges with respect to memory allocation and computation. Consequently, the

discretion to designate the accumulation period is delegated to the user, enabling the accommodation of data availability and constraints imposed by computational resources. In any case, the robustness of all ReBeLS-derived products is amenable to quantitative evaluation through an analysis of the associated uncertainties. Thus, providing an objective assessment of the quality and reliability of ReBeLS-derived products.

## 5. Conclusions

This paper describes the ReBeLS algorithm designed and implemented in response to the Copernicus programme's (through its Global Land Service, Climate Change Service and the COPA project) need for a state-of-the-art algorithm to retrieve the land surface BRDF from Sentinel-3 A and Sentinel-3B surface reflectance observations.

We show that ReBeLS overcomes various shortcomings present in the
<figure><img src="image_17.png" /><figcaption>Fig. 14. The layout is the same as Fig. 13 but for day 154 of 2019 (summer).</figcaption></figure>

The parallelised mode of ReBeLS has been successfully integrated into the operational processing chain of the Copernicus Global Land and Climate Change Services, which generates operational Sentinel-3-based BRDF-adjusted surface reflectances – further used to compute an NDVI product – and surface albedo products, respectively (León-Tavares, 2022; Sánchez-Zapero et al., 2023).

At the time of writing this article, ReBeLS remains confined to proprietary access, thereby granting the scientific community access to ReBeLS-derived products encompassing NDVI and surface albedo via the Copernicus Global Land and Copernicus Climate change services, respectively. It is important to acknowledge that the quality of ReBeLS outcomes is intimately linked to the quality and temporal extent of the surface reflectance time series Nonetheless, we acknowledge that the computational processing of long-term time series, such as those spanning multiple years, could potentially introduce technical challenges with respect to memory allocation and computation. Consequently, the

discretion to designate the accumulation period is delegated to the user, enabling the accommodation of data availability and constraints imposed by computational resources. In any case, the robustness of all ReBeLS-derived products is amenable to quantitative evaluation through an analysis of the associated uncertainties. Thus, providing an objective assessment of the quality and reliability of ReBeLS-derived products.

## 5. Conclusions

This paper describes the ReBeLS algorithm designed and implemented in response to the Copernicus programme's (through its Global Land Service, Climate Change Service and the COPA project) need for a state-of-the-art algorithm to retrieve the land surface BRDF from Sentinel-3 A and Sentinel-3B surface reflectance observations.

We show that ReBeLS overcomes various shortcomings present in the
current global BRDF products, namely, the need for a backup algorithm that dominates the frequency of BRDF inversion solutions at northern latitudes, the lack of uncertainty quantification and the accumulation of surface reflectance observations over a fixed period. These advantages are at the cost of ReBeLS requiring a dynamic prior BRDF database – currently built from 7 years of MODIS MCD43 product, but any other sensor and period could be used – and assuming a temporal regularisation parameter.

We found that for a region comprising Europe, the BRDF model retrieved with ReBeLS can reproduce within uncertainties (|zsc| < 1σ), the directional effects in the VIIRS bidirectional reflectance dataset VNP09GA. ReBeLS generates BRDF-adjusted reflectance images where the directional effects visible in daily MVC mosaics are efficiently removed. Additionally, the levels of BRDF-adjusted reflectance are consistent – fitted slopes approaching one with statistically significant (p<0.05) and high correlations coefficients (r > 0.9) – with the nadir BRDF-adjusted reflectance product from VIIRS (VNP43IA4).

Since Sentinel-3's radiometers are foreseen to operate for more than a decade, the ReBeLS algorithm within Copernicus Global Land and Climate Change services operational chains appears as a reliable methodology to generate BRDF descriptors at a global scale. It allows the generation of BRDF-corrected vegetation (León-Tavares, 2022; León-Tavares et al., 2021) and phenological (Jin and Eklundh, 2014) indices along with surface albedo products (Sánchez-Zapero et al., 2023). The asset of ReBeLS is that it provides BRDF-descriptors with an associated uncertainty. This seems particularly relevant to state-of-the-art methodologies for atmospheric correction (Yin et al., 2022) and retrievals (Lorente et al., 2018; Qin et al., 2019; Tilstra et al., 2021). Moreover, the uncertainty quantification exercised within ReBeLS is in line with the current ESA's calibration and validation strategy to enhance Earth observation products traceability (Gorroño et al., 2017; Niro et al., 2021) where the upcoming TRUTHS mission (Fox and Green, 2020) is expected to play a pivotal role.

Given that the methodology implemented in ReBeLS is sensor
agnostic, a scenario wherein both Sentinel-3 satellites are compromised
would not pose overwhelming challenges. This can be attributed to the
structural design of the ReBeLS algorithm, which allows for a seamless
transition to alternative wide imaging swath satellite data streams with
minimal implementation complexities. The adaptable nature of ReBeLS
ensures that integrating data from alternate sources would not necessi-
tate a substantial investment of resources and effort, thereby effectively
mitigating the potential for prolonged disruptions in the provision of
near-real-time services. Future work may extend ReBeLS to a multi-
sensor approach, where harmonised datasets can be used concomi-
tantly to retrieve surface anisotropy. Additionally, a separate treatment
for snow pixels will be integrated. Finally, the dynamic BRDF descriptors

prior may be extended to comprise the full MODIS MCD43 V006.1
archive.

CRediT authorship contribution statement

**Jonathan León-Tavares:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. **Jose Gómez-Dans:** Conceptualization, Writing – review & editing. **Jean-Louis Roujean:** Writing – review & editing. **Véronique Bruniquel:** Funding acquisition, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgments

We acknowledge the anonymous referees for a very careful reading of the manuscript and comments that significantly improved this manuscript. JLT would like to thank Prof. Crystal Schaaf and Dr. Zhuosen Wang for their insightful comments and suggestions while reviewing the Copernicus Global Land Service ReBeLS ATBD. JGD would like to acknowledge financial support from the European Union Horizon 2020 research and innovation programme under grant agreement No 687320 MULTIPLY (MULTIscale SENTINEL land surface information retrieval Platform). JGD was supported by the Natural Environment Research Council's (NERC) National Centre for Earth Observation (NCEO) (project number 525861). The Sentinel-3 atmospherically corrected surface reflectance datasets were generated by the Global component of the Land Service of Copernicus, the Earth Observation program of the European Commission. The algorithm design and implementation of ReBeLS were performed in the framework of the Copernicus Global Land Service and Copernicus Climate Change Service. The validation of the ReBeLS algorithm was carried out within the Copernicus 4 Core Products Algorithm Studies (COPA) project funded by the EU and ESA contract No. 4000133521/20/I-BG. JLT appreciates the technical assistance from Dirk Daems and Yves Verheijen at the VITO computing centre throughout the development of the ReBeLS algorithm.

Appendix A. Determination of the temporal regularisation parameter

The smoothness of the BRDF descriptors' time evolution constraint is defined by the regularisation parameter, γ, representing our confidence in this smoothness constraint. If the regularisation parameter is too large, the variability of the BRDF descriptors' temporal trajectories will be over-smoothed. If, on the other hand, too little regularisation is imposed on the solution by adopting a small regularisation parameter value, then the fit will be good, and the improvement on conditioning reduced. The quality of the solution relies largely on the quality of the observations, with possible low model fit residuals but unstable trajectories (i.e. rapid and step temporal variations due to processing artefacts, depending on the BRDF descriptors retrievals).

The approach adopted in ReBeLS to solve the BRDF inverse problem, as stated in Eq. (3), requires the knowledge of the temporal regularisation parameter. However, there is no consensus on determining the optimal value for the regularisation parameter (Quaife and Lewis, 2010; Samain et al., 2008; Wang et al., 2007; Yin et al., 2022). Operational constraints hamper a suitable analytical approach (Zobitz et al., 2020). Nevertheless, estimating the regularisation parameter in an inverse problem can be guided using heuristics (Lewis et al., 2012a; Valentine and Sambridge, 2018). In ReBeLS, we determine the regularisation parameter using a cross-validation approach, leveraging on the BRDF descriptors' ability to predict the observations of an independent sensor. In other words, we sweep over a grid of regularisation parameter values and predict surface reflectance observations from another sensor acquired over the same area and time with a similar spectral configuration but differing in its view/illumination geometry. Finally, we select the value from the regularisation parameter grid that delivers the best predictive performance.

We have selected PROBA-V as a reference sensor because one of its three different spatial resolutions (Dierckx et al., 2014; Wolters et al., 2018) matches the Sentinel-3 OLCI pixel size of 333 m. Hereafter, we simply refer to the PROBA-V 333 m observations as the PROBA-V dataset. We can
current global BRDF products, namely, the need for a backup algorithm that dominates the frequency of BRDF inversion solutions at northern latitudes, the lack of uncertainty quantification and the accumulation of surface reflectance observations over a fixed period. These advantages are at the cost of ReBeLS requiring a dynamic prior BRDF database – currently built from 7 years of MODIS MCD43 product, but any other sensor and period could be used – and assuming a temporal regularisation parameter.

We found that for a region comprising Europe, the BRDF model retrieved with ReBeLS can reproduce within uncertainties (|zsc| < 1σ), the directional effects in the VIIRS bidirectional reflectance dataset VNP09GA. ReBeLS generates BRDF-adjusted reflectance images where the directional effects visible in daily MVC mosaics are efficiently removed. Additionally, the levels of BRDF-adjusted reflectance are consistent – fitted slopes approaching one with statistically significant (p<0.05) and high correlations coefficients (r > 0.9) – with the nadir BRDF-adjusted reflectance product from VIIRS (VNP43IA4).

Since Sentinel-3's radiometers are foreseen to operate for more than a decade, the ReBeLS algorithm within Copernicus Global Land and Climate Change services operational chains appears as a reliable methodology to generate BRDF descriptors at a global scale. It allows the generation of BRDF-corrected vegetation (León-Tavares, 2022; León-Tavares et al., 2021) and phenological (Jin and Eklundh, 2014) indices along with surface albedo products (Sánchez-Zapero et al., 2023). The asset of ReBeLS is that it provides BRDF-descriptors with an associated uncertainty. This seems particularly relevant to state-of-the-art methodologies for atmospheric correction (Yin et al., 2022) and retrievals (Lorente et al., 2018; Qin et al., 2019; Tilstra et al., 2021). Moreover, the uncertainty quantification exercised within ReBeLS is in line with the current ESA's calibration and validation strategy to enhance Earth observation products traceability (Gorroño et al., 2017; Niro et al., 2021) where the upcoming TRUTHS mission (Fox and Green, 2020) is expected to play a pivotal role.

Given that the methodology implemented in ReBeLS is sensor
agnostic, a scenario wherein both Sentinel-3 satellites are compromised
would not pose overwhelming challenges. This can be attributed to the
structural design of the ReBeLS algorithm, which allows for a seamless
transition to alternative wide imaging swath satellite data streams with
minimal implementation complexities. The adaptable nature of ReBeLS
ensures that integrating data from alternate sources would not necessi-
tate a substantial investment of resources and effort, thereby effectively
mitigating the potential for prolonged disruptions in the provision of
near-real-time services. Future work may extend ReBeLS to a multi-
sensor approach, where harmonised datasets can be used concomi-
tantly to retrieve surface anisotropy. Additionally, a separate treatment
for snow pixels will be integrated. Finally, the dynamic BRDF descriptors

prior may be extended to comprise the full MODIS MCD43 V006.1
archive.

CRediT authorship contribution statement

**Jonathan León-Tavares:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. **Jose Gómez-Dans:** Conceptualization, Writing – review & editing. **Jean-Louis Roujean:** Writing – review & editing. **Véronique Bruniquel:** Funding acquisition, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgments

We acknowledge the anonymous referees for a very careful reading of the manuscript and comments that significantly improved this manuscript. JLT would like to thank Prof. Crystal Schaaf and Dr. Zhuosen Wang for their insightful comments and suggestions while reviewing the Copernicus Global Land Service ReBeLS ATBD. JGD would like to acknowledge financial support from the European Union Horizon 2020 research and innovation programme under grant agreement No 687320 MULTIPLY (MULTIscale SENTINEL land surface information retrieval Platform). JGD was supported by the Natural Environment Research Council's (NERC) National Centre for Earth Observation (NCEO) (project number 525861). The Sentinel-3 atmospherically corrected surface reflectance datasets were generated by the Global component of the Land Service of Copernicus, the Earth Observation program of the European Commission. The algorithm design and implementation of ReBeLS were performed in the framework of the Copernicus Global Land Service and Copernicus Climate Change Service. The validation of the ReBeLS algorithm was carried out within the Copernicus 4 Core Products Algorithm Studies (COPA) project funded by the EU and ESA contract No. 4000133521/20/I-BG. JLT appreciates the technical assistance from Dirk Daems and Yves Verheijen at the VITO computing centre throughout the development of the ReBeLS algorithm.

Appendix A. Determination of the temporal regularisation parameter

The smoothness of the BRDF descriptors' time evolution constraint is defined by the regularisation parameter, γ, representing our confidence in this smoothness constraint. If the regularisation parameter is too large, the variability of the BRDF descriptors' temporal trajectories will be over-smoothed. If, on the other hand, too little regularisation is imposed on the solution by adopting a small regularisation parameter value, then the fit will be good, and the improvement on conditioning reduced. The quality of the solution relies largely on the quality of the observations, with possible low model fit residuals but unstable trajectories (i.e. rapid and step temporal variations due to processing artefacts, depending on the BRDF descriptors retrievals).

The approach adopted in ReBeLS to solve the BRDF inverse problem, as stated in Eq. (3), requires the knowledge of the temporal regularisation parameter. However, there is no consensus on determining the optimal value for the regularisation parameter (Quaife and Lewis, 2010; Samain et al., 2008; Wang et al., 2007; Yin et al., 2022). Operational constraints hamper a suitable analytical approach (Zobitz et al., 2020). Nevertheless, estimating the regularisation parameter in an inverse problem can be guided using heuristics (Lewis et al., 2012a; Valentine and Sambridge, 2018). In ReBeLS, we determine the regularisation parameter using a cross-validation approach, leveraging on the BRDF descriptors' ability to predict the observations of an independent sensor. In other words, we sweep over a grid of regularisation parameter values and predict surface reflectance observations from another sensor acquired over the same area and time with a similar spectral configuration but differing in its view/illumination geometry. Finally, we select the value from the regularisation parameter grid that delivers the best predictive performance.

We have selected PROBA-V as a reference sensor because one of its three different spatial resolutions (Dierckx et al., 2014; Wolters et al., 2018) matches the Sentinel-3 OLCI pixel size of 333 m. Hereafter, we simply refer to the PROBA-V 333 m observations as the PROBA-V dataset. We can
quantify the agreement between datasets by assuming a linear form of the relation between the original and predicted (from the retrieved Sentinel-3 ReBeLS BRDF) PROBA-V surface reflectance time series (July 2018 to June 2019). Atmospherically corrected surface reflectances in the NIR and Oa17 spectral channels – less sensitive to aerosol contamination and inadequate atmospheric correction than visible channels - from PROBA-V and Sentinel-3, respectively, are considered to perform the cross-validation exercise.

<figure><img src="image_18.png" /><figcaption>**Fig. A1.** Top panel: Relation between observed and predicted PROBA-V surface reflectances in the NIR spectral channel. The predicted reflectances are obtained by using the Sentinel-3 BRDF-descriptors retrieved using different values of γ. Bottom panel: The slope-γ relation obtained for the AERONET Harvard forest station. The different values of γ within the grid considered are colour coded.</figcaption></figure>

The top panel of Fig. A1 shows the relations between the original and predicted PROBA-V surface reflectance time series for different values of γ for a pixel near the AERONET Harvard forest station (lat = 42.532, lon = -72.188). We assume a linear relationship of the form $ρ_{predicted} = slope × ρ_{observed} + offset$ between the original and predicted surface reflectance datasets. The best model that fitted an orthogonal slope considering uncertainties in both datasets (Akritas and Bershady, 1996; Nemmen et al., 2012) is shown as a function of the regularisation parameter, γ, value in the bottom panel of Fig. A1. The closer a slope value approaches one, the better the correspondence between observed and predicted datasets.

The slopes and predicted datasets associated with each regularisation parameter value in the grid are colour coded in Fig. A1, where the lowest and highest regularisation parameter values in the explored grid are displayed in dark blue and dark red colours, respectively. We find, for this situation, that the best-fit slope increases with the regularisation parameter until it reaches values close to one around $γ \sim 5 × 10^{-3}$ and after that, the relation flattens, becoming proportionally inverse for $γ \ge 5e-2$. For low values of the regularisation parameter, the estimated BRDF descriptors result in a poor model of the observed PROBA-V reflectance. The latter is due to the inversion resulting in unstable estimates of the BRDF parameters due to poor conditioning. Consequently, these overfitted BRDF descriptors cannot reproduce the directional effects in the PROBA-V dataset. For large
quantify the agreement between datasets by assuming a linear form of the relation between the original and predicted (from the retrieved Sentinel-3 ReBeLS BRDF) PROBA-V surface reflectance time series (July 2018 to June 2019). Atmospherically corrected surface reflectances in the NIR and Oa17 spectral channels – less sensitive to aerosol contamination and inadequate atmospheric correction than visible channels - from PROBA-V and Sentinel-3, respectively, are considered to perform the cross-validation exercise.

<figure><img src="image_18.png" /><figcaption>**Fig. A1.** Top panel: Relation between observed and predicted PROBA-V surface reflectances in the NIR spectral channel. The predicted reflectances are obtained by using the Sentinel-3 BRDF-descriptors retrieved using different values of γ. Bottom panel: The slope-γ relation obtained for the AERONET Harvard forest station. The different values of γ within the grid considered are colour coded.</figcaption></figure>

The top panel of Fig. A1 shows the relations between the original and predicted PROBA-V surface reflectance time series for different values of γ for a pixel near the AERONET Harvard forest station (lat = 42.532, lon = -72.188). We assume a linear relationship of the form $ρ_{predicted} = slope × ρ_{observed} + offset$ between the original and predicted surface reflectance datasets. The best model that fitted an orthogonal slope considering uncertainties in both datasets (Akritas and Bershady, 1996; Nemmen et al., 2012) is shown as a function of the regularisation parameter, γ, value in the bottom panel of Fig. A1. The closer a slope value approaches one, the better the correspondence between observed and predicted datasets.

The slopes and predicted datasets associated with each regularisation parameter value in the grid are colour coded in Fig. A1, where the lowest and highest regularisation parameter values in the explored grid are displayed in dark blue and dark red colours, respectively. We find, for this situation, that the best-fit slope increases with the regularisation parameter until it reaches values close to one around $γ \sim 5 × 10^{-3}$ and after that, the relation flattens, becoming proportionally inverse for $γ \ge 5e-2$. For low values of the regularisation parameter, the estimated BRDF descriptors result in a poor model of the observed PROBA-V reflectance. The latter is due to the inversion resulting in unstable estimates of the BRDF parameters due to poor conditioning. Consequently, these overfitted BRDF descriptors cannot reproduce the directional effects in the PROBA-V dataset. For large
regularisation parameter values, the temporal dynamics are over-smoothed and eventually collapse into a horizontal line.

As shown in Fig. A1, the low (flat) slope values obtained while using either too low or too high $\gamma$ values can be taken as evidence that the derived BRDF retrievals trajectories do not adequately represent the intrinsic anisotropy of the land surface reflectance; hence a disagreement with PROBA-V observations arises. On the other hand, from the bottom panel of Fig. A1, it is noticeable that a range of $\gamma$ values spanning approximately from $5 \times 10^{-3}$ to $5 \times 10^{-2}$ deliver predicted datasets for which the slopes are steeper (approaching one) and do not change significantly. This range of regularisation parameters corresponds to a trade-off region, and selecting the lowest $\gamma$ value within this range is a conservative choice that minimises the amount of regularisation while still having an adequate predictive capability.

<figure><img src="image_19.png" /><figcaption>Fig. A2. The slope-$\gamma$ distribution for a subset of 879 stations from the AERONET sample. Individual profiles are shown as solid lines. The mode of the slope distribution and the 68% ($\pm 1\sigma$) range are displayed as filled circles error bars, respectively.</figcaption></figure>

Fig. A2 shows the slope-$\gamma$ relation derived by repeating the cross-validation exercise performed on the Harvard forest station to an ensemble of 879 stations drawn from the AERONET sample. The slope-$\gamma$ profiles for each station in the sample are shown as solid lines, while the slope distribution mode is displayed as a filled circle where the error bars represent the 16% and 84% percentiles – region defining 68% ($\pm 1\sigma$) of all slopes. Fig. A2 shows that most of the stations in our sample follow the slope-$\gamma$ relation obtained for the Harvard forest station shown in Fig. A1, namely that the lowest values for slopes are obtained when either $\gamma$ is too low or too high and that a well-defined range of slopes close to one can be found approximately in the range $\gamma = [5 \times 10^{-3}, 1 \times 10^{-1}]$. Based on these experiments, we find statistical evidence that adding a regularisation parameter higher than $5 \times 10^{-3}$ does not significantly modify favourably the slope between the observed and predicted PROBA-V datasets. We select a regularisation parameter $\gamma = 5 \times 10^{-3}$, based on the observations made above.

## Appendix B. ReBeLS inversion implementation details

The steps that ReBeLS carries out within the Inversion stage (see Fig. 1) can be summarised as listed below (all symbols are defined in Eq. (3)):

Step A. Input is a surface reflectance (atmospherically corrected) time series along with associated uncertainties, acquisition geometries ($\Omega_{\text{sensor}}$, $\Omega_{\text{Sun}}$) and quality layers. There is no constraint on the extent of the input time series; for the examples presented in section 4, only one year (July 2018 to June 2019) is used. Similar to any EO land surface processor, the algorithm begins by excluding pixels identified as cloud, shadow, or of poor quality based on the spectral band quality layers. Despite the aforementioned filtering, certain outliers, often linked to undetected shadows or clouds, may persist. These outliers are discarded from subsequent processing via an outlier detection step. A comprehensive description of the filtering an outlier identification can be found in the CGLS ReBeLS ATBD.<fnref n="20" /> ReBeLS thus rearranges the selected surface reflectance dataset data to build vector $\mathbf{R}$ and matrix $C_{\text{obs}}^{-1}$. Further, the RossThickLiSparseReciprocal angular functions are computed to populate the operator matrix $\mathbf{H}$ allowing to assemble the $J_{\text{obs}}(\mathbf{x})$ term.

Step B1. ReBeLS loads the MCD43P BRDF descriptors ($x_{\text{prior}}, C_{\text{prior}}^{-1}$) corresponding to the week(s) of the year comprised in the surface reflectance time series to build the $J_{\text{prior}}(\mathbf{x})$ component.

1. Step B2. According to the temporal extent of the surface reflectance time series, ReBeLS assembles the $\mathbf{B}$ matrix, where the number of rows is defined by the number of days in the surface reflectance time series. Together with the static regularisation parameter $\gamma$, the $J_{\text{model}}(\mathbf{x})$ term is assembled.

2. Step B3. Once the three components of the cost functions are assembled, the maximum a posteriori solution which is the value of $\mathbf{x}$ that yields the maximum of $P(\mathbf{x}|\mathbf{R})$, or equivalently the minimum of the cost function $J(\mathbf{x})$ is obtained by solving

$$ (\mathbf{H}^T \mathbf{C}_{\text{obs}}^{-1} \mathbf{H} + \mathbf{C}_{\text{prior}}^{-1} + \gamma \mathbf{B}^T \mathbf{B}) \mathbf{x} = (\mathbf{H}^T \mathbf{C}_{\text{obs}}^{-1} \mathbf{R} + \mathbf{C}_{\text{prior}}^{-1} \mathbf{x}_{\text{prior}}) \quad (8) $$

where the covariance matrix $\mathbf{C}_{\text{x}}$ is given by

<footnote n="20">²⁰ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_BRDFCorrection300m-V1_11.10.pdf</footnote>
regularisation parameter values, the temporal dynamics are over-smoothed and eventually collapse into a horizontal line.

As shown in Fig. A1, the low (flat) slope values obtained while using either too low or too high $\gamma$ values can be taken as evidence that the derived BRDF retrievals trajectories do not adequately represent the intrinsic anisotropy of the land surface reflectance; hence a disagreement with PROBA-V observations arises. On the other hand, from the bottom panel of Fig. A1, it is noticeable that a range of $\gamma$ values spanning approximately from $5 \times 10^{-3}$ to $5 \times 10^{-2}$ deliver predicted datasets for which the slopes are steeper (approaching one) and do not change significantly. This range of regularisation parameters corresponds to a trade-off region, and selecting the lowest $\gamma$ value within this range is a conservative choice that minimises the amount of regularisation while still having an adequate predictive capability.

<figure><img src="image_19.png" /><figcaption>Fig. A2. The slope-$\gamma$ distribution for a subset of 879 stations from the AERONET sample. Individual profiles are shown as solid lines. The mode of the slope distribution and the 68% ($\pm 1\sigma$) range are displayed as filled circles error bars, respectively.</figcaption></figure>

Fig. A2 shows the slope-$\gamma$ relation derived by repeating the cross-validation exercise performed on the Harvard forest station to an ensemble of 879 stations drawn from the AERONET sample. The slope-$\gamma$ profiles for each station in the sample are shown as solid lines, while the slope distribution mode is displayed as a filled circle where the error bars represent the 16% and 84% percentiles – region defining 68% ($\pm 1\sigma$) of all slopes. Fig. A2 shows that most of the stations in our sample follow the slope-$\gamma$ relation obtained for the Harvard forest station shown in Fig. A1, namely that the lowest values for slopes are obtained when either $\gamma$ is too low or too high and that a well-defined range of slopes close to one can be found approximately in the range $\gamma = [5 \times 10^{-3}, 1 \times 10^{-1}]$. Based on these experiments, we find statistical evidence that adding a regularisation parameter higher than $5 \times 10^{-3}$ does not significantly modify favourably the slope between the observed and predicted PROBA-V datasets. We select a regularisation parameter $\gamma = 5 \times 10^{-3}$, based on the observations made above.

## Appendix B. ReBeLS inversion implementation details

The steps that ReBeLS carries out within the Inversion stage (see Fig. 1) can be summarised as listed below (all symbols are defined in Eq. (3)):

Step A. Input is a surface reflectance (atmospherically corrected) time series along with associated uncertainties, acquisition geometries ($\Omega_{\text{sensor}}$, $\Omega_{\text{Sun}}$) and quality layers. There is no constraint on the extent of the input time series; for the examples presented in section 4, only one year (July 2018 to June 2019) is used. Similar to any EO land surface processor, the algorithm begins by excluding pixels identified as cloud, shadow, or of poor quality based on the spectral band quality layers. Despite the aforementioned filtering, certain outliers, often linked to undetected shadows or clouds, may persist. These outliers are discarded from subsequent processing via an outlier detection step. A comprehensive description of the filtering an outlier identification can be found in the CGLS ReBeLS ATBD.<fnref n="20" /> ReBeLS thus rearranges the selected surface reflectance dataset data to build vector $\mathbf{R}$ and matrix $C_{\text{obs}}^{-1}$. Further, the RossThickLiSparseReciprocal angular functions are computed to populate the operator matrix $\mathbf{H}$ allowing to assemble the $J_{\text{obs}}(\mathbf{x})$ term.

Step B1. ReBeLS loads the MCD43P BRDF descriptors ($x_{\text{prior}}, C_{\text{prior}}^{-1}$) corresponding to the week(s) of the year comprised in the surface reflectance time series to build the $J_{\text{prior}}(\mathbf{x})$ component.

1. Step B2. According to the temporal extent of the surface reflectance time series, ReBeLS assembles the $\mathbf{B}$ matrix, where the number of rows is defined by the number of days in the surface reflectance time series. Together with the static regularisation parameter $\gamma$, the $J_{\text{model}}(\mathbf{x})$ term is assembled.

2. Step B3. Once the three components of the cost functions are assembled, the maximum a posteriori solution which is the value of $\mathbf{x}$ that yields the maximum of $P(\mathbf{x}|\mathbf{R})$, or equivalently the minimum of the cost function $J(\mathbf{x})$ is obtained by solving

$$ (\mathbf{H}^T \mathbf{C}_{\text{obs}}^{-1} \mathbf{H} + \mathbf{C}_{\text{prior}}^{-1} + \gamma \mathbf{B}^T \mathbf{B}) \mathbf{x} = (\mathbf{H}^T \mathbf{C}_{\text{obs}}^{-1} \mathbf{R} + \mathbf{C}_{\text{prior}}^{-1} \mathbf{x}_{\text{prior}}) \quad (8) $$

where the covariance matrix $\mathbf{C}_{\text{x}}$ is given by

<footnote n="20">²⁰ https://land.copernicus.eu/global/sites/cgls.vito.be/files/products/CGLOPS1_ATBD_BRDFCorrection300m-V1_11.10.pdf</footnote>
$$C_x = (H^T C_{obs}^{-1} H + C_{prior}^{-1} + \gamma B^T B)^{-1} \qquad (9)$$

Solving Eq. 8 for *x* yields a vector comprising the BRDF descriptors for each day of the surface reflectance time series extent. The variances associated to the BRDF descriptors retrievals are obtained from the diagonal elements of the covariance matrix $C_x$.

## References

Akritas, M.G., Bershady, M.A., 1996. Linear regression for astronomical data with measurement errors and intrinsic scatter. Astrophys. J. https://doi.org/10.1086/177901.

Breon, F.M., Maignan, F., 2017. A BRDF-BPDF database for the analysis of earth target reflectances. Earth Syst. Sci. Data 9, 31–45. https://doi.org/10.5194/essd-9-31-2017.

Bréon, F., Vermote, E., 2012. Remote sensing of environment correction of MODIS surface reflectance time series for BRDF effects. Remote Sens. Environ. 125, 1–9. https://doi.org/10.1016/j.rse.2012.06.025.

Bréon, F.M., Vermote, E., Murphy, E.F., Franch, B., 2015. Measuring the directional variations of land surface reflectance from MODIS. IEEE Trans. Geosci. Remote Sens. 53, 4638–4649. https://doi.org/10.1109/TGRS.2015.2405344.

Carrer, D., Moparthy, S., Lellouch, G., Ceamanos, X., Pinault, F., Freitas, S.C., Trigo, I.F., 2018. Land surface albedo derived on a ten daily basis from meteosat second generation observations: the NRT and climate data record collections from the EUMETSAT LSA SAF. Remote Sens. 10 https://doi.org/10.3390/RS10081262 page 1262 10, 1262.

Carrer, D., Pinault, F., Lellouch, G., Trigo, I.F., Benhadj, I., Camacho, F., Ceamanos, X., Moparthy, S., Munoz-sabater, J., Schüller, L., Sánchez-zapero, J., 2021. Surface albedo retrieval from 40-years of earth observations through the eumetsat/lsa saf and eu/c3s programmes: the versatile algorithm of pyalus. Remote Sens. 13, 1–32. https://doi.org/10.3390/rs13030372.

Combal, B., Baret, F., Weiss, M., Trubuil, A., Macé, D., Pragnère, A., Myneni, R., Knyazikhin, Y., Wang, L., 2003. Retrieval of canopy biophysical variables from bidirectional reflectance using prior information to solve the ill-posed inverse problem. Remote Sens. Environ. 84, 1–15. https://doi.org/10.1016/S0034-4257(02)00035-4.

Dierckx, W., Sterckx, S., Benhadj, I., Livens, S., Duhoux, G., van Achteren, T., Francois, M., Mellab, K., Saint, G., 2014. PROBA-V Mission for Global Vegetation Monitoring: Standard Products and Image Quality, 35, pp. 2589–2614. https://doi.org/10.1080/014311614.2014.883097.

Fox, N., Green, P., 2020. Traceable radiometry underpinning terrestrial- and Helios- studies (TRUTHS): an element of a space-based climate and calibration observatory. Remote Sens. 12 https://doi.org/10.3390/RS12152400 page 2400 12, 2400.

Gao, F., Jin, Y., Li, X., Schaaf, C.B., 2002. Bidirectional NDVI and Atmospherically Resistant BRDF Inversion for Vegetation Canopy, 40, pp. 1269–1278.

Gastellu-Etchegorry, J.P., Lauret, N., Tavares, L., Lamquin, N., Bruniquel, V., Roujean, J. L., Hagolle, O., et al., 2022. Correction of 3D Directional Effects in Sentinel-2 and 3 Images with Sentinel-3 Time Series and Dart 3D Radiative Transfer Model. IEEE, pp. 4563–4566.

GDAL/OGR contributors, 2022. GDAL/OGR Geospatial Data Abstraction software Library. https://doi.org/10.5281/zenodo.5884351.

Geiger, B., Carrer, D., Franchistéguy, L., Roujean, J.L., Meurey, C., 2008. Land surface albedo derived on a daily basis from meteosat second generation observations. IEEE Trans. Geosci. Remote Sens. 46, 3841–3856. https://doi.org/10.1109/TGRS.2008.2001798.

Gómez-Dans, J.L., Lewis, P.E., Disney, M., 2016. Efficient emulation of radiative transfer codes using gaussian processes and application to land surface parameter inferences. Remote Sens. 8, 1–32. https://doi.org/10.3390/rs8020119.

Gorroño, J., Fomferra, N., Peters, M., Gascon, F., Underwood, C.I., Fox, N.P., Kirches, G., Brockmann, C., 2017. A radiometric uncertainty tool for the Sentinel 2 mission. Remote Sens. 9 https://doi.org/10.3390/RS9020178 page 178 9, 178.

Holben, B.N., 1986. Characteristics of maximum-value composite images from temporal AVHRR data. Int. J. Remote Sens. 7, 1417–1434. https://doi.org/10.1080/01431168608948945.

Jin, H., Eklundh, L., 2014. A physically based vegetation index for improved monitoring of plant phenology. Remote Sens. Environ. 152, 512–525. https://doi.org/10.1016/J.RSE.2014.07.010.

Kimes, D.S., Knyazikhin, Y., Privette, J.L., Abuelgasim, A.A., Gao, F., 2000. Inversion methods for physically-based models. Remote Sens. Rev. 18, 381–439. https://doi.org/10.1080/02757250009532396.

Klisch, A., Atzberger, C., 2016. Operational drought monitoring in Kenya using MODIS NDVI time series. Remote Sens. 8 https://doi.org/10.3390/RS8040267 page 267 8, 267.

León-Tavares, J., 2022. “CGLOPS-1” Framework Service Contract N° 199494 (JRC) Algorithm Theoretical Basis Document BRDF Model Retrieval From Sentinel-3 Collection 300M Version 1 Issue II.10. Copernicus Global Land Service ATBD.

León-Tavares, J., Roujean, J.L., Smets, B., Wolters, E., Toté, C., Swinnen, E., 2021. Correction of directional effects in vegetation ndvi time-series. Remote Sens. 13, 1–23. https://doi.org/10.3390/rs13061130.

Lewis, P., Gómez-Dans, J., Kaminski, T., Settle, J., Quaife, T., Gobron, N., Styles, J., Berger, M., 2012a. An earth observation land data assimilation system (EO-LDAS). Remote Sens. Environ. 120, 219–235. https://doi.org/10.1016/j.rse.2011.12.027.

Lewis, P., Guanter, L., Lopez Saldana, G., Muller, J.P., Watson, G., Shane, N., Kennedy, T., Fisher, J., Domenech, C., Preusker, R., North, P., Heckel, A., Danne, O.,

Kramer, U., Zuhlke, M., Fomferra, N., Brockmann, C., Schaaf, C., 2012b. The ESA globAlbedo project: algorithm. Int. Geosci. Remote Sens. Symp. 5745–5748. https://doi.org/10.1109/IGARSS.2012.6352306.

Li, X., Strahler, A.H., 1985. Geometric-optical bidirectional reflectance modeling of a conifer forest canopy. IEEE Trans. Geosci. Remote Sens. GE-24, 705–721. https://doi.org/10.1109/tgrs.1986.289706.

Liu, Y., Wang, Z., Sun, Q., Erb, A.M., Li, Z., Schaaf, C.B., Zhang, X., Román, M.O., Scott, R.L., Zhang, Q., Novick, K.A., Syndonia Bret-Harte, M., Petroy, S., SanClements, M., 2017. Evaluation of the VIIRS BRDF, albedo and NBAR products suite and an assessment of continuity with the long term MODIS record. Remote Sens. Environ. 201, 256–274. https://doi.org/10.1016/j.rse.2017.09.020.

Lorente, A., Folkert Boersma, K., Stammes, P., Gijsbert Tilstra, L., Richter, A., Yu, H., Kharbouche, S., Muller, J.P., 2018. The importance of surface reflectance anisotropy for cloud and NO2 retrievals from GOME-2 and OMI. Atmos. Meas. Tech. 11, 4509–4529. https://doi.org/10.5194/AMT-11-4509-2018.

Lucht, W., Lewis, P., 2000. Theoretical noise sensitivity of BRDF and albedo retrieval from the EOS-MODIS and MISR sensors with respect to angular sampling. Int. J. Remote Sens. 21, 81–98. https://doi.org/10.1080/014311600211000.

Maignan, F., Breon, F.-M., Lacaze, R., 2004. Bidirectional Reflectance of Earth Targets: Evaluation of Analytical Models Using a Large Set of Spaceborne Measurements with Emphasis on the Hot Spot, 90, pp. 210–220. https://doi.org/10.1016/j.rse.2003.12.006.

Meroni, M., Fasbender, D., Rembold, F., Atzberger, C., Klisch, A., 2019. Near real-time vegetation anomaly detection with MODIS NDVI: timeliness vs. accuracy and effect of anomaly computation options. Remote Sens. Environ. 221, 508–521. https://doi.org/10.1016/j.RSE.2018.11.041.

Muller, J., Lewis, P., Fischer, J., North, P., Framer, U., 2011. The ESA GlobAlbedo project for mapping the earth's land surface albedo for 15 years from European sensors. Geophys. Res. Abstr. 13, 4–5.

Muller, J.-P., Kennedy, T., Lewis, P., Kharbouche, S., North, P., Fisher, J., Preusker, R., Brockmann, C., Kramer, U., Danne, O., Fomferra, N., Leonard, O., Pinnock, S., 2013. GlobAlbedo Algorithm Theoretical Basis Document, 313.

Nemmen, R.S., Georgopoulos, M., Guiriec, S., Meyer, E.T., Gehrels, N., Sambruna, R. M., 2012. A universal scaling for the energetics of relativistic jets from black hole systems. Science 1979 (338), 1445–1448. https://doi.org/10.1126/science.1227416/SUPPL_FILE/1227416.NEMMEN.SM.REVISION.1.PDF.

Nicodemus, F.E., Richmond, J.C., Hsia, J.J., Ginsberg, I.W., Limeris, T., 1977. Geometrical considerations and nomenclature for reflectance. Natl. Bur. Stand. (US) Monogr. 1–52.

Niro, F., Goryl, P., Dransfeld, S., Boccia, V., Gascon, F., Adams, J., Themann, B., Scifoni, S., Doxani, G., 2021. European Space Agency (ESA) calibration/validation strategy for optical land-imaging satellites and pathway towards interoperability. Remote Sens. 13 https://doi.org/10.3390/RS13153003. Page 3003 13, 3003.

Privette, J.L., Eck, T.F., Deering, D.W., 1997. Estimating spectral albedo and nadir reflectance through inversion of simple BRDF models with AVHRR/MODIS-like data. J. Geophys. Res. Atmos. 102 (D24), 29529–29542.

Qin, W., Fasnacht, Z., Haffner, D., Vasilikov, A., Joiner, J., Krotkov, N., Fisher, B., Spurr, R., 2019. A geometry-dependent surface Lambertian-equivalent reflectivity product for UV-vis retrievals - part 1: evaluation over land surfaces using measurements from OMI at 466 nm. Atmos. Meas. Tech. 12, 3997–4017. https://doi.org/10.5194/AMT-12-3997-2019.

Quaife, T., Lewis, P., 2010. Temporal constraints on linear BRDF model parameters. IEEE Trans. Geosci. Remote Sens. 48, 2445–2450. https://doi.org/10.1109/TGRS.2009.2038901.

Rahman, H., Dedieu, G., 1994. SMAC: a simplified method for the atmospheric correction of satellite measurements in the solar spectrum. Int. J. Remote Sens. 15, 123–143. https://doi.org/10.1080/01431169408954055.

Roger, J.C., Vermote, E.F., Devadiga, S., Ray, J.P., 2016. Suomi-NPP VIIRS Surface Reflectance User's Guide. V1 Re-processing (NASA Land SIPS).

Ross, J.K., 1981. The radiation regime and architecture of plant stands. J. Ecol. https://doi.org/10.2307/2259995.

Roujean, J.L., 2017. Inversion of lumped parameters using BRDF kernels. In: Comprehensive Remote Sensing. Elsevier, pp. 23–34. https://doi.org/10.1016/B978-0-12-409548-9.10346-X.

Roujean, J.L., Leroy, M., Deschamps, P.Y., 1992a. A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data. J. Geophys. Res. Atmos. 97, 20455–20468. https://doi.org/10.1029/92JD01411.

Roujean, J.L., Leroy, M., Podaire, A., Deschamps, P.Y., 1992b. Evidence of surface reflectance bidirectional effects from a noaa/avhrr multi-thermal data set. Int. J. Remote Sens. 13, 685–698. https://doi.org/10.1080/01431169208904146.

Roujean, J., Leon-Tavares, J., Smets, B., Claes, P., Camacho, F., Coca, D., Sanchez-zapero, J., 2018. Remote sensing of environment surface albedo and toc-r 300 m products from PROBA-V instrument in the framework of copernicus global land service. Remote Sens. Environ. 215, 57–73. https://doi.org/10.1016/j.rse.2018.05.015.

Samain, O., Roujean, J.L., Geiger, B., 2008. Use of a Kalman filter for the retrieval of surface BRDF coefficients with a time-evolving model based on the ECOCLIMAP land
$$C_x = (H^T C_{obs}^{-1} H + C_{prior}^{-1} + \gamma B^T B)^{-1} \qquad (9)$$

Solving Eq. 8 for *x* yields a vector comprising the BRDF descriptors for each day of the surface reflectance time series extent. The variances associated to the BRDF descriptors retrievals are obtained from the diagonal elements of the covariance matrix *C*<sub>*x*</sub>.

## References

Akritas, M.G., Bershady, M.A., 1996. Linear regression for astronomical data with measurement errors and intrinsic scatter. Astrophys. J. https://doi.org/10.1086/177901.

Breon, F.M., Maignan, F., 2017. A BRDF-BPDF database for the analysis of earth target reflectances. Earth Syst. Sci. Data 9, 31–45. https://doi.org/10.5194/essd-9-31-2017.

Bréon, F., Vermote, E., 2012. Remote sensing of environment correction of MODIS surface reflectance time series for BRDF effects. Remote Sens. Environ. 125, 1–9. https://doi.org/10.1016/j.rse.2012.06.025.

Bréon, F.M., Vermote, E., Murphy, E.F., Franch, B., 2015. Measuring the directional variations of land surface reflectance from MODIS. IEEE Trans. Geosci. Remote Sens. 53, 4638–4649. https://doi.org/10.1109/TGRS.2015.2405344.

Carrer, D., Moparthy, S., Lellouch, G., Ceamanos, X., Pinault, F., Freitas, S.C., Trigo, I.F., 2018. Land surface albedo derived on a ten daily basis from meteosat second generation observations: the NRT and climate data record collections from the EUMETSAT LSA SAF. Remote Sens. 10 https://doi.org/10.3390/RS10081262 page 1262 10, 1262.

Carrer, D., Pinault, F., Lellouch, G., Trigo, I.F., Benhadj, I., Camacho, F., Ceamanos, X., Moparthy, S., Munoz-sabater, J., Schüller, L., Sánchez-zapero, J., 2021. Surface albedo retrieval from 40-years of earth observations through the eumetsat/lsa saf and eu/c3s programmes: the versatile algorithm of pyalus. Remote Sens. 13, 1–32. https://doi.org/10.3390/rs13030372.

Combal, B., Baret, F., Weiss, M., Trubuil, A., Macé, D., Pragnère, A., Myneni, R., Knyazikhin, Y., Wang, L., 2003. Retrieval of canopy biophysical variables from bidirectional reflectance using prior information to solve the ill-posed inverse problem. Remote Sens. Environ. 84, 1–15. https://doi.org/10.1016/S0034-4257(02)00035-4.

Dierckx, W., Sterckx, S., Benhadj, I., Livens, S., Duhoux, G., van Achteren, T., Francois, M., Mellab, K., Saint, G., 2014. PROBA-V Mission for Global Vegetation Monitoring: Standard Products and Image Quality, 35, pp. 2589–2614. https://doi.org/10.1080/014311614.2014.883097.

Fox, N., Green, P., 2020. Traceable radiometry underpinning terrestrial- and Helios- studies (TRUTHS): an element of a space-based climate and calibration observatory. Remote Sens. 12 https://doi.org/10.3390/RS12152400 page 2400 12, 2400.

Gao, F., Jin, Y., Li, X., SchAAF, C.B., 2002. Bidirectional NDVI and Atmospherically Resistant BRDF Inversion for Vegetation Canopy, 40, pp. 1269–1278.

Gastellu-Etchegorry, J.P., Lauret, N., Tavares, L., Lamquin, N., Bruniquel, V., Roujean, J. L., Hagolle, O., et al., 2022. Correction of Directional Effects in Sentinel-2 and 3 Images with Sentinel-3 Time Series and Dart 3D Radiative Transfer Model. IEEE, pp. 4563–4566.

GDAL/OGR contributors, 2022. GDAL/OGR Geospatial Data Abstraction software Library. https://doi.org/10.5281/zenodo.5884351.

Geiger, B., Carrer, D., Franchistéguy, L., Roujean, J.L., Meurey, C., 2008. Land surface albedo derived on a daily basis from meteosat second generation observations. IEEE Trans. Geosci. Remote Sens. 46, 3841–3856. https://doi.org/10.1109/TGRS.2008.2001798.

Gómez-Dans, J.L., Lewis, P.E., Disney, M., 2016. Efficient emulation of radiative transfer codes using gaussian processes and application to land surface parameter inferences. Remote Sens. 8, 1–32. https://doi.org/10.3390/rs8020119.

Gorroño, J., Fomferra, N., Peters, M., Gascon, F., Underwood, C.I., Fox, N.P., Kirches, G., Brockmann, C., 2017. A radiometric uncertainty tool for the Sentinel 2 mission. Remote Sens. 9 https://doi.org/10.3390/RS9020178 page 178 9, 178.

Holben, B.N., 1986. Characteristics of maximum-value composite images from temporal AVHRR data. Int. J. Remote Sens. 7, 1417–1434. https://doi.org/10.1080/01431168608948945.

Jin, H., Eklundh, L., 2014. A physically based vegetation index for improved monitoring of plant phenology. Remote Sens. Environ. 152, 512–525. https://doi.org/10.1016/J.RSE.2014.07.010.

Kimes, D.S., Knyazikhin, Y., Privette, J.L., Abuelgasim, A.A., Gao, F., 2000. Inversion methods for physically-based models. Remote Sens. Rev. 18, 381–439. https://doi.org/10.1080/02757250009532396.

Klisch, A., Atzberger, C., 2016. Operational drought monitoring in Kenya using MODIS NDVI time series. Remote Sens. 8 https://doi.org/10.3390/RS8040267 page 267 8, 267.

León-Tavares, J., 2022. “CGLOPS-1” Framework Service Contract N° 199494 (JRC) Algorithm Theoretical Basis Document BRDF Model Retrieval From Sentinel-3 Collection 300M Version 1 Issue II.10. Copernicus Global Land Service ATBD.

León-Tavares, J., Roujean, J.L., Smets, B., Wolters, E., Toté, C., Swinnen, E., 2021. Correction of directional effects in vegetation ndvi time-series. Remote Sens. 13, 1–23. https://doi.org/10.3390/rs13061130.

Lewis, P., Gómez-Dans, J., Kaminski, T., Settle, J., Quaife, T., Gobron, N., Styles, J., Berger, M., 2012a. An earth observation land data assimilation system (EO-LDAS). Remote Sens. Environ. 120, 219–235. https://doi.org/10.1016/j.rse.2011.12.027.

Lewis, P., Guanter, L., Lopez Saldana, G., Muller, J.P., Watson, G., Shane, N., Kennedy, T., Fisher, J., Domenech, C., Preusker, R., North, P., Heckel, A., Danne, O.,

Kramer, U., Zuhlke, M., Fomferra, N., Brockmann, C., SchAAF, C., 2012b. The ESA globAlbedo project: algorithm. Int. Geosci. Remote Sens. Symp. 5745–5748. https://doi.org/10.1109/IGARSS.2012.6352306.

Li, X., Strahler, A.H., 1985. Geometric-optical bidirectional reflectance modeling of a conifer forest canopy. IEEE Trans. Geosci. Remote Sens. GE-24, 705–721. https://doi.org/10.1109/tgrs.1986.289706.

Liu, Y., Wang, Z., Sun, Q., Erb, A.M., Li, Z., SchAAF, C.B., Zhang, X., Román, M.O., Scott, R.L., Zhang, Q., Novick, K.A., Syndonia Bret-Harte, M., Petroy, S., SanClements, M., 2017. Evaluation of the VIIRS BRDF, albedo and NBAR products suite and an assessment of continuity with the long term MODIS record. Remote Sens. Environ. 201, 256–274. https://doi.org/10.1016/j.rse.2017.09.020.

Lorente, A., Folkert Boersma, K., Stammes, P., Gijsbert Tilstra, L., Richter, A., Yu, H., Kharbouche, S., Muller, J.P., 2018. The importance of surface reflectance anisotropy for cloud and NO2 retrievals from GOME-2 and OMI. Atmos. Meas. Tech. 11, 4509–4529. https://doi.org/10.5194/AMT-11-4509-2018.

Lucht, W., Lewis, P., 2000. Theoretical noise sensitivity of BRDF and albedo retrieval from the EOS-MODIS and MISR sensors with respect to angular sampling. Int. J. Remote Sens. 21, 81–98. https://doi.org/10.1080/014311600211000.

Maignan, F., Breon, F.-M., Lacaze, R., 2004. Bidirectional Reflectance of Earth Targets: Evaluation of Analytical Models Using a Large Set of Spaceborne Measurements with Emphasis on the Hot Spot, 90, pp. 210–220. https://doi.org/10.1016/j.rse.2003.12.006.

Meroni, M., Fasbender, D., Rembold, F., Atzberger, C., Klisch, A., 2019. Near real-time vegetation anomaly detection with MODIS NDVI: timeliness vs. accuracy and effect of anomaly computation options. Remote Sens. Environ. 221, 508–521. https://doi.org/10.1016/j.RSE.2018.11.041.

Muller, J.L., Lewis, P., Fischer, J., North, P., Framer, U., 2011. The ESA GlobAlbedo project for mapping the earth's land surface albedo for 15 years from European sensors. Geophys. Res. Abstr. 13, 4–5.

Muller, J.-P., Kennedy, T., Lewis, P., Kharbouche, S., North, P., Fisher, J., Preusker, R., Brockmann, C., Kramer, U., Danne, O., Fomferra, N., Leonard, O., Pinnock, S., 2013. GlobAlbedo Algorithm Theoretical Basis Document, 313.

Nemmen, R.S., Georgopoulos, M., Guiriec, S., Meyer, E.T., Gehrels, N., Sambruna, R. M., 2012. A universal scaling for the energetics of relativistic jets from black hole systems. Science 1979 (338), 1445–1448. https://doi.org/10.1126/science.1227416/SUPPL_FILE/1227416.NEMMEN.SM.REVISION.1.PDF.

Nicodemus, F.E., Richmond, J.C., Hsia, J.J., Ginsberg, I.W., Limeris, T., 1977. Geometrical considerations and nomenclature for reflectance. Natl. Bur. Stand. (US) Monogr. 1–52.

Niro, F., Goryl, P., Dransfeld, S., Boccia, V., Gascon, F., Adams, J., Themann, B., Scifoni, S., Doxani, G., 2021. European Space Agency (ESA) calibration/validation strategy for optical land-imaging satellites and pathway towards interoperability. Remote Sens. 13 https://doi.org/10.3390/RS13153003. Page 3003 13, 3003.

Privette, J.L., Eck, T.F., Deering, D.W., 1997. Estimating spectral albedo and nadir reflectance through inversion of simple BRDF models with AVHRR/MODIS-like data. J. Geophys. Res. Atmos. 102 (D24), 29529–29542.

Qin, W., Fasnacht, Z., Haffner, D., Vasilikov, A., Joiner, J., Krotkov, N., Fisher, B., Spurr, R., 2019. A geometry-dependent surface Lambertian-equivalent reflectivity product for UV-vis retrievals - part 1: evaluation over land surfaces using measurements from OMI at 466 nm. Atmos. Meas. Tech. 12, 3997–4017. https://doi.org/10.5194/AMT-12-3997-2019.

Quaife, T., Lewis, P., 2010. Temporal constraints on linear BRDF model parameters. IEEE Trans. Geosci. Remote Sens. 48, 2445–2450. https://doi.org/10.1109/TGRS.2009.2038901.

Rahman, H., Dedieu, G., 1994. SMAC: a simplified method for the atmospheric correction of satellite measurements in the solar spectrum. Int. J. Remote Sens. 15, 123–143. https://doi.org/10.1080/01431169408954055.

Roger, J.C., Vermote, E.F., Devadiga, S., Ray, J.P., 2016. Suomi-NPP VIIRS Surface Reflectance User's Guide. V1 Re-processing (NASA Land SIPS).

Ross, J.K., 1981. The radiation regime and architecture of plant stands. J. Ecol. https://doi.org/10.2307/2259995.

Roujean, J.L., 2017. Inversion of lumped parameters using BRDF kernels. In: Comprehensive Remote Sensing. Elsevier, pp. 23–34. https://doi.org/10.1016/B978-0-12-409548-9.10346-X.

Roujean, J.L., Leroy, M., Deschamps, P.Y., 1992a. A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data. J. Geophys. Res. Atmos. 97, 20455–20468. https://doi.org/10.1029/92JD01411.

Roujean, J.L., Leroy, M., Podaire, A., Deschamps, P.Y., 1992b. Evidence of surface reflectance bidirectional effects from a noaa/avhrr multi-thermal data set. Int. J. Remote Sens. 13, 685–698. https://doi.org/10.1080/01431169208904146.

Roujean, J., Leon-Tavares, J., Smets, B., Claes, P., Camacho, F., Coca, D., Sanchez-zapero, J., 2018. Remote sensing of environment surface albedo and toc-r 300 m products from PROBA-V instrument in the framework of copernicus global land service. Remote Sens. Environ. 215, 57–73. https://doi.org/10.1016/j.rse.2018.05.015.

Samain, O., Roujean, J.L., Geiger, B., 2008. Use of a Kalman filter for the retrieval of surface BRDF coefficients with a time-evolving model based on the ECOCLIMAP land
cover classification. Remote Sens. Environ. 112, 1337–1346. https://doi.org/
10.1016/j.rse.2007.07.007.

Sánchez-Zapero, J., Camacho, F., Martínez-Sánchez, E., Gorroño, J., León-Tavares, J., Benhadj, I., Toté, C., Swinnen, E., Muñoz-Sabater, J., 2023. Global estimates of surface albedo from Sentinel-3 OLCI and SLSTR data for copernicus climate change service: algorithm and preliminary validation. Remote Sens. Environ. 287, 113460 https://doi.org/10.1016/j.rse.2023.113460.

Schaaf, C.B., Strahler, A.H., Lucht, W., 2002. First Operational BRDF, Albedo Nadir Reflectance Products from MODIS.

Schaaf, C.B., Liu, J., Gao, F., Strahler, A.H., 2011. Aqua and Terra MODIS Albedo and Reflectance Anisotropy Products. Land Remote Sensing and Global Environmental Change: NASA's Earth Observing System and the Science of ASTER and MODIS, Remote Sensing and Digital Image Processing. https://doi.org/10.1007/978-1-4419-6749-7_24.

Schaepman-Strub, G., Schaepman, M.E., Painter, T.H., Dangel, S., Martonchik, J.V., 2006. Reflectance quantities in optical remote sensing-definitions and case studies. Remote Sens. Environ. 103, 27–42. https://doi.org/10.1016/j.rse.2006.03.002.

Strugnell, N.C., Lucht, W., Schaaf, C., 2001. A global albedo data set derived from AVHRR data for use in climate simulations. Geophys. Res. Lett. 28, 191–194. https://doi.org/10.1029/2000GL011580.

Tarantola, A., 2004. Inverse Problem Theory, and Methods for Model Parameter Estimation.

Tilstra, L.G., Tuinder, O.N.E., Wang, P., Stammes, P., 2021. Directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface measured by the GOME-2 satellite instruments. Atmos. Meas. Tech. 14, 4219–4238. https://doi.org/
10.5194/AMT-14-4219-2021.

Valentine, A.P., Sambridge, M., 2018. Optimal regularization for a class of linear inverse problem. Geophys. J. Int. 215, 1003–1021. https://doi.org/10.1093/GJI/GGY303.

Vermote, E., Justice, C.O., Bréon, F.M., 2009. Towards a generalized approach for correction of the BRDF effect in MODIS directional reflectances. IEEE Trans. Geosci. Remote Sens. 47, 898–908. https://doi.org/10.1109/TGRS.2008.2005977.

Viovy, N., Arino, O., Belward, A.S., 2007. The Best Index Slope Extraction (BISE): A Method for Reducing Noise in NDVI Time-Series, 13, pp. 1585–1590. https://doi.org/10.1080/01431169208904212.

Wang, Y., Li, X., Nashed, Z., Zhao, F., Yang, H., Guan, Y., Zhang, H., 2007. Regularized kernel-based BRDF model inversion method for ill-posed land surface parameter retrieval. Remote Sens. Environ. 111, 36–50. https://doi.org/10.1016/j.rse.2007.03.007.

Wang, Z., Schaaf, C.B., Sun, Q., Shuai, Y., Román, M.O., 2018. Capturing rapid land surface dynamics with collection V006 MODIS BRDF/NBAR/albedo (MCD43) products. Remote Sens. Environ. 207, 50–64. https://doi.org/10.1016/j.rse.2018.02.001.

Wanner, W., Li, X., Strahler, A.H., 1995. On the derivation of kernels for kernel-driven models of bidirectional reflectance. J. Geophys. Res. Atmos. 100, 21077–21089. https://doi.org/10.1029/95JD02371.

Wolters, E., Dierckx, W., Iordache, M.-D., Swinnen, E., 2018. PROBA-V Products User Manual. VITO.

Yin, F., Lewis, P.E., Gómez-Dans, J.L., 2022. Bayesian atmospheric correction over land: Sentinel-2/MSI and Landsat 8/OLI. Geosci. Model Dev. 15, 7933–7976. https://doi.org/10.5194/GMD-15-7933-2022.

Zeng, L., Wardlow, B.D., Xiang, D., Hu, S., Li, D., 2020. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens. Environ. 237, 111511 https://doi.org/10.1016/J.RSE.2019.111511.

Zobitz, J.M., Quaife, T., Nichols, N.K., 2020. Efficient hyper-parameter determination for regularised linear BRDF parameter retrieval. Int. J. Remote Sens. 41, 1437–1457. https://doi.org/10.1080/01431161.2019.1667552.
cover classification. Remote Sens. Environ. 112, 1337–1346. https://doi.org/
10.1016/j.rse.2007.07.007.

Sánchez-Zapero, J., Camacho, F., Martínez-Sánchez, E., Gorroño, J., León-Tavares, J., Benhadj, I., Toté, C., Swinnen, E., Muñoz-Sabater, J., 2023. Global estimates of surface albedo from Sentinel-3 OLCI and SLSTR data for copernicus climate change service: algorithm and preliminary validation. Remote Sens. Environ. 287, 113460
https://doi.org/10.1016/j.rse.2023.113460.

Schaaf, C.B., Strahler, A.H., Lucht, W., 2002. First Operational BRDF, Albedo Nadir Reflectance Products from MODIS.

Schaaf, C.B., Liu, J., Gao, F., Strahler, A.H., 2011. Aqua and Terra MODIS Albedo and Reflectance Anisotropy Products. Land Remote Sensing and Global Environmental Change: NASA's Earth Observing System and the Science of ASTER and MODIS, Remote Sensing and Digital Image Processing. https://doi.org/10.1007/978-1-4419-6749-7_24.

Schaepman-Strub, G., Schaepman, M.E., Painter, T.H., Dangel, S., Martonchik, J.V., 2006. Reflectance quantities in optical remote sensing-definitions and case studies. Remote Sens. Environ. 103, 27–42. https://doi.org/10.1016/j.rse.2006.03.002.

Strugnell, N.C., Lucht, W., Schaaf, C., 2001. A global albedo data set derived from AVHRR data for use in climate simulations. Geophys. Res. Lett. 28, 191–194.
https://doi.org/10.1029/2000GL011580.

Tarantola, A., 2004. Inverse Problem Theory, and Methods for Model Parameter Estimation.

Tilstra, L.G., Tuinder, O.N.E., Wang, P., Stammes, P., 2021. Directionally dependent Lambertian-equivalent reflectivity (DLER) of the Earth's surface measured by the GOME-2 satellite instruments. Atmos. Meas. Tech. 14, 4219–4238. https://doi.org/
10.5194/AMT-14-4219-2021.

Valentine, A.P., Sambridge, M., 2018. Optimal regularization for a class of linear inverse problem. Geophys. J. Int. 215, 1003–1021. https://doi.org/10.1093/GJI/GGY303.

Vermote, E., Justice, C.O., Bréon, F.M., 2009. Towards a generalized approach for correction of the BRDF effect in MODIS directional reflectances. IEEE Trans. Geosci. Remote Sens. 47, 898–908. https://doi.org/10.1109/TGRS.2008.2005977.

Viovy, N., Arino, O., Belward, A.S., 2007. The Best Index Slope Extraction (BISE): A Method for Reducing Noise in NDVI Time-Series, 13, pp. 1585–1590. https://doi.org/10.1080/01431169208904212.

Wang, Y., Li, X., Nashed, Z., Zhao, F., Yang, H., Guan, Y., Zhang, H., 2007. Regularized kernel-based BRDF model inversion method for ill-posed land surface parameter retrieval. Remote Sens. Environ. 111, 36–50. https://doi.org/10.1016/j.rse.2007.03.007.

Wang, Z., Schaaf, C.B., Sun, Q., Shuai, Y., Román, M.O., 2018. Capturing rapid land surface dynamics with collection V006 MODIS BRDF/NBAR/albedo (MCD43) products. Remote Sens. Environ. 207, 50–64. https://doi.org/10.1016/j.rse.2018.02.001.

Wanner, W., Li, X., Strahler, A.H., 1995. On the derivation of kernels for kernel-driven models of bidirectional reflectance. J. Geophys. Res. Atmos. 100, 21077–21089.
https://doi.org/10.1029/95JD02371.

Wolters, E., Dierckx, W., Iordache, M.-D., Swinnen, E., 2018. PROBA-V Products User Manual. VITO.

Yin, F., Lewis, P.E., Gómez-Dans, J.L., 2022. Bayesian atmospheric correction over land: Sentinel-2/MSI and Landsat 8/OLI. Geosci. Model Dev. 15, 7933–7976. https://doi.org/10.5194/GMD-15-7933-2022.

Zeng, L., Wardlow, B.D., Xiang, D., Hu, S., Li, D., 2020. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens. Environ. 237, 111511 https://doi.org/10.1016/J.RSE.2019.111511.

Zobitz, J.M., Quaife, T., Nichols, N.K., 2020. Efficient hyper-parameter determination for regularised linear BRDF parameter retrieval. Int. J. Remote Sens. 41, 1437–1457. https://doi.org/10.1080/01431161.2019.1667552.