# Surface Albedo Retrieval from 40-Years of Earth Observations through the EUMETSAT/LSA SAF and EU/C3S Programmes: The Versatile Algorithm of PYALUS

Dominique Carrer, Florian Pinault, Gabriel Lellouch, Isabel F. Trigo, Iskander Benhadj, Fernando Camacho, Xavier Ceamanos, Suman Moparthy, Joaquin Munoz-Sabater, Lothar Schüller, et al.

► To cite this version:

Dominique Carrer, Florian Pinault, Gabriel Lellouch, Isabel F. Trigo, Iskander Benhadj, et al. "Surface Albedo Retrieval from 40-Years of Earth Observations through the EUMETSAT/LSA SAF and EU/C3S Programmes: The Versatile Algorithm of PYALUS." Remote Sensing, 2021, 13, 10.3390/rs13030372. insu-04833621

HAL Id: insu-04833621

https://insu.hal.science/insu-04833621v1

Submitted on 12 Dec 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>

Article

# Surface Albedo Retrieval from 40-Years of Earth Observations through the EUMETSAT/LSA SAF and EU/C3S Programmes: The Versatile Algorithm of PYALUS

Dominique Carrer<fnref n="1" />,*, Florian Pinault<fnref n="1" />, Gabriel Lellouch<fnref n="1" />, Isabel F. Trigo<fnref n="2" /> ID, Iskander Benhadj<fnref n="3" />, Fernando Camacho<fnref n="4" />, Xavier Ceamanos<fnref n="1" />, Suman Moparthy<fnref n="1" />,<fnref n="5" /> ID, Joaquin Munoz-Sabater<fnref n="6" /> ID, Lothar Schüller<fnref n="7" /> and Jorge Sánchez-Zapero<fnref n="4" /> ID

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

check for
updates

**Citation:** Carrer, D.; Pinault, F.; Lellouch, G.; Trigo, I.F.; Benhadj, I.; Camacho, F.; Ceamanos, X.; Moparthy, S.; Munoz-Sabater, J.; Schüller, L.; et al. 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.* **2021**, *13*, 372. https://doi.org/10.3390/rs13030372

Received: 8 December 2020

Accepted: 12 January 2021

Published: 21 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Land surface albedo quantifies the fraction of the sunlight reflected by the surface of the Earth. This article presents the algorithm concepts for the remote sensing of this variable based on the heritage of several developments which were performed at Méteo France over the last decade and described in several papers by Carrer et al. The scientific algorithm comprises four steps: an atmospheric correction, a sensor harmonisation (optional), a BRDF (Bidirectional Reflectance Distribution Function) inversion, and the albedo calculation. At the time being, the method has been applied to 11 sensors in the framework of two European initiatives (Satellite Application Facility on Land Surface Analysis—LSA SAF, and Copernicus Climate Change Service—C3S): NOAA-7-9-11-14-16-17/AVHRR2-3, SPOT/VGT1-2, Metop/AVHRR-3, PROBA-V, and MSG/SEVIRI. This work leads to a consistent archive of almost 40 years of satellite-derived albedo data (available in 2020). From a single sensor, up to three different albedo products with different characteristics have been developed to address the requirements of both, near real-time (NRT) (weather prediction with a demand of timeliness of 1 h) and climate communities. The evaluation of the algorithm applied to different platforms was recently made by Lellouch et al. and Sánchez Zapero et al. in 2020 which can be considered as companion papers. After a summary of the method for the retrieval of these surface albedos, this article describes the specificities of each retrieval, lists the differences, and discusses the limitations. The plan of continuity with the next European satellite missions and perspectives of improvements are introduced. For example, Metop/AVHRR-3 albedo will soon become the medium resolution sensor product with the longest NRT data record, since MODIS is approaching the end of its life-cycle. Additionally, Metop-SG/METimage will ensure its continuity thanks to consistent production of data sets guaranteed till 2050 by the member states of the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT). In the end, the common strategy which we proposed through the different programmes may offer an unprecedented opportunity to study the temporal trends affecting surface properties and to analyse human-induced climate change. Finally, the access to the source code (called PYALUS) is provided through an open access platform in order to share with the community the expertise on the satellite retrieval of this variable.

<footnote n="1">¹ Centre National de Recherches Météorologiques (CNRM), Université de Toulouse, Météo France, CNRS 42, Avenue Gaspard Coriolis, 31057 Toulouse, France; florian.pinault@meteo.fr (F.P.); gabriel.lellouch@meteo.fr (G.L.); xavier.ceamanos@meteo.fr (X.C.); sumankasyap@gmail.com (S.M.)</footnote>

<footnote n="2">² Instituto Portugues do Mar e da Atmosfera (IPMA), Rua C do Aeroporto, 1749-077 Lisbon, Portugal; isabel.trigo@ipma.pt</footnote>

<footnote n="3">³ Flemish Institute for Technological Research (VITO), Remote Sensing Unit, Boeretang 200, B-2400 Mol, Belgium; iskander.benhadj@vito.be</footnote>

<footnote n="4">⁴ Earth Observation Laboratory (EOLAB), Parc Científic University of Valencia, C/Catedràtic Agustín Escardino, 9, 46980 Paterna, Valencia, Spain; fernando.camacho@eolab.es (F.C.); jorge.sanchez@eolab.es (J.S.-Z.)</footnote>

<footnote n="5">⁵ ACRI-ST, BP 234, 06904 Sophia Antipolis CEDEX, France</footnote>

<footnote n="6">⁶ ECMWF, Shinfield Park, Reading RG2 9AX, UK; joaquin.Munoz@ecmwf.int</footnote>

<footnote n="7">⁷ EUMETSAT, Eumetsat Allee 1, 64295 Darmstadt, Germany; Lothar.Schueller@eumetsat.int</footnote>

<footnote>* Correspondence: dominique.carrer@meteo.fr</footnote>
**Keywords:** albedo; vegetation properties; MSG; VGT; EPS; AVHRR; PROBA-V; BRDF; radiative transfer; LSA SAF; EUMETSAT; Copernicus; C3S

# 1. Introduction

Land surface albedo is defined as the ratio of the solar radiation flux reflected by the surface to the total incoming solar radiation flux. In other words, land surface albedo quantifies the fraction of the sunlight reflected by the surface of the Earth. Typical values are between 0.05 for dense vegetation areas and 0.8 for areas covered by snow [1] i.e., fresh snow reflects around 80% of the incoming solar radiation in the shortwave (0.3–4.0 µm), whereas forests only reflect 5% of the same incoming quantity. In that sense, deforestations over boreal forests have been put forward as a strategy leading to a beneficial cooling effect within a climate change mitigation framework [2]. Surface albedo can be considered as fundamental quantity of the surface energy budget. For example, *Pielke and Avissar* [3] found that a 4% increase in land surface albedo would result in a decrease of around 0.7 °C in the Earth’s equilibrium temperature.

Different concepts and definitions of the land surface albedo coexist. In remote sensing, surface albedo is basically the mean value of the changes of surface brightness within the visible/shortwave domain, as observed by a remote sensor averaged over a relatively short composite period (usually between one day and one month). The main causes of surface albedo changes are the variation of viewing and illumination angular geometries, the change of the diffuse fraction of the atmosphere, the evolution of the leaves and bare soil optical properties, and the change of roughness, soil moisture or the architecture of the canopy. In the 1980s, the mathematical BRDF (Bidirectional Reflectance Distribution Function) concept was applied by Ross et al. [4] to natural surfaces in order to reconstruct their brightness for all the possible geometrical configurations (viewing and illumination), even if these configurations are very far from the observed configuration of the given sensor. The common principle of these algorithms which use a kernel-based approach is to use a set of semi-physical (or semi-empirical) kernel functions that simulate the surface brightness as series of geometrical contributions representing the mean brightness (isotropic kernel function), the scattering of the radiation within the leaves of a forest (volumetric kernel function), the shadowing effects (geometric kernel function), and other effects. The remote sensing observations are fitted on these series of kernel functions (weighted by isotropic, volumetric, geometric, and other coefficients) to estimate the type of surface (crop, forest, etc.) that is observed, and consecutively, to assess its anisotropy for all the non-observed or non-illuminated configurations. The angular integration of the series of kernel functions then provides the mean value of the surface brightness. By doing so, we incorporate all possible viewing and illumination angles, i.e., we are not limited to the geometries of the acquired observations. As a remote sensing observation is usually done at a thin discrete range of wavelength corresponding to the spectral response of the sensor, spectral integration (also called narrow to broadband conversion) is required to derive land surface albedo for the entire shortwave domain, usually defined between 0.3 and 4 µm. The knowledge of the overall brightness integrated in the shortwave domain is important as it is the optical domain in which almost all the solar energy is received. The quantity that is not reflected by the surface in the shortwave domain drives how much energy is absorbed by the surface which is one of the main input variables of the energy, carbon, and water cycles at the surface. Additionally, the most relevant surface albedo quantity for the energy budget refers usually to the total shortwave broadband (BB) interval, comprising the VISible (VIS) plus Near InfraRed (NIR) wavelength ranges. It is important to highlight that in addition to serving as an intermediate product for deriving the broadband albedo quantities, the spectral estimates contain also a wealth of information about the physical state of the surface. This information can be used for a variety of purposes such as
vegetation monitoring and land cover classification. In turn, this also constitutes important elements for setting up adequate surface modelling schemes [5,6].

As the surface is usually non-Lambertian, the brightness of the surface changes along the day according to its directional reflectivity properties. Additionally, two concepts were introduced to define the land surface albedo due to the motion of the Sun with respect to the surface: the bi-hemispherical albedo and the directional-hemispherical albedo [7]. The directional-hemispherical albedo, also called "black-sky albedo", is the integration of the BRDF over the viewing hemisphere, considering a fixed illumination direction (chosen arbitrarily by convention as the local solar noon). The quantity is assumed to be non-dependent on the sensor geometry. The bi-hemispherical albedo, also called "white-sky albedo", completes the previous integration with an additional angular integration on the zenith angles on the full illumination hemisphere. It is usually comparable to the observed value of albedo for cloudy sky conditions (i.e., under significant diffuse radiation). The Global Climate Observing System (GCOS) specified the black-sky albedo as an Essential Climate Variable [8].

The kernel-based approach was adopted for the development of albedo products from a number of multi-angular space-borne sensors. These include POLDER, SeaWiFS, VEGETATION, and MODIS (e.g., [9–15]). Other approaches were used for MISR [16] and Meteosat [17,18]. Besides, a number of initiatives including GlobAlbedo, [19,20] and Global LAnd Surface Satellites (GLASS) [21], have provided long-term albedo products from these satellite data. To go one step forward, exploiting the combination of kernel models to represent both the surface and aerosol optical directional reflectivities enable the simultaneous retrieving of aerosols in the atmosphere and surface albedos [22,23].

Nowadays, 40 years after the work proposed by Ross et al. [4], many developments have been performed to improve the atmospheric correction, the spectral harmonisation, the representation of the BRDF kernel functions, e.g., [24–28], and the inversion process [1,23,29,30]. However, the state of the art is still to use the BRDF approach for the computation of the surface albedo. It offers a pragmatic and cost-effective solution to the problem of inversion of surface brightness properties within the framework of exploitation programmes constrained by the requirement to process an extensive amount of satellite data in near real-time (NRT). The increasing use of machine learning methods [31,32] appears to open interesting areas of improvements. However, the counterpart is that these methods often need BRDF data for quite a long training period, while their robustness and capability to evolve toward a state never observed in the training period may always be questionable. The latter may be a limitation, for example, to detect long-term trends due to climate change.

The main objective of the two following important European initiatives from European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) and the European Union (EU), is to provide added-value products for the meteorological and environmental science communities, especially in the fields of climate modelling, environmental management, natural hazards management, and climate change detection. Among many other products, EUMETSAT and EU initiated the development of land surface albedos originating from different Earth Observation instruments.

The EUMETSAT maintains a number of decentralized processing centres—called Satellite Application Facilities (SAFs)—dedicated to different application areas. The Portuguese Institute for the Sea and the Atmosphere, IPMA (or ‘Instituto Português do Mar e da Atmosfera’), coordinates and hosts the Satellite Application Facility on Land Surface Analysis (LSA SAF). The primary objective of the LSA SAF is to provide added-value products describing land surface aspects for the meteorological and environmental operational and science communities [33]. Data from different instruments are routinely processed in near real-time by the LSA SAF operational system, i.e., ensuring users receive satellite products up to three hours after the last satellite observation required for their generation (https://www.eumetsat.int/lsa-saf). Long term reanalyses were also performed [30]. Presently, the delivered operational products ranging from land surface
albedo and temperature, to downwelling radiation fluxes, vegetation parameters and evapotranspiration, amongst many others. The initiative started in 1999 with research and development activities and is currently in its third *Continuous Development* and Operations Phase (CDOP 3). The LSA SAF programme is expected to continue to use instruments on board EUMETSAT satellites that were, or will be, launched between 2002 and well beyond 2030.

The objective of the Copernicus Climate Change Service (C3S) of the EU is to deliver, among others, three Essential Climate Variables (ECVs) that are fAPAR (fraction of Absorbed radiation in the PAR domain which is the photosynthetically active domain for the vegetation comprised between 0.4 and 0.7 µm), LAI (Leaf Area Index), and surface albedo. In comparison to the LSA SAF that covers the period 2004–2022, this initiative has no NRT constrain. The first phase (312a-Lot9) of the programme began in November 2016, and ended in December 2019. An overlap period occurred with the second phase that started in April 2018 and will end in June 2021. The VITO institute in Belgium is hosting the operational system. The major challenge is the harmonisation of the datasets originating from various sensors to produce a consistent timeseries of the ECV from the 1980s until now. This consistent dataset is also called a Climate Data Record (CDR). Data from the Advanced Very High Resolution Radiometer (AVHRR) 7 to 19 on board NOAA series, SPOT/VGT1 & 2, and PROBA-V sensors are exploited for this purpose. In the future, data from Sentinel-3/SLSTR & OLCI will also be included. Directly connected to C3S, the Copernicus Global Land Service (CGLS) component operates "a multi-purpose service component" that also provides series of bio-geophysical products on the status and evolution of land surface at the global scale. This initiative is focused on NRT and has users from different communities (as agriculture) than C3S (which is dedicated to climate community). Surface albedo products were developed in C3S and CGLS. In June 2020, discussions between C3S and CGLS have agreed to rationalize developments of the scientific algorithms and to use brokered data from one programme to the other. Additionally, because surface albedo is more related to climate applications, this variable will be, in the future, only produced in C3S (and not in CGLS anymore).

In the context of these important European initiatives from EUMETSAT and the EU, several surface albedo algorithms are developed by Météo France to exploit data acquired by the Spinning Enhanced Visible and Infrared Imager (SEVIRI) radiometer embarked on Meteosat Second Generation (MSG), AVHRR on board the series of Metop satellites, the AVHRR (7 to 16) series on board the NOAA satellites, the SPOT-VEGETATION series, and PROBA-V sensor. These data correspond to an almost 40-year archive of albedo products characterizing the properties of the surface. From 2022 and for the next 30 years, these will be completed by data acquired by the Flexible Combined Imager (FCI) on board Meteosat Third Generation (MTG), and by METImage and 3MI sensors on board EUMETSAT Polar System programme—Second Generation satellites (Metop-SG), all operated by EUMETSAT.

This article reviews the efforts made by Météo France for the retrieval of surface albedo through different Earth observing systems. The objective is also to share with the community the expertise on the satellite retrieval of this variable by giving the access to the source code through an open access platform. The evaluation of the algorithm applied to different platforms was recently made by Lellouch et al. [34] and Sánchez Zapero et al. [35], which can be considered as companion papers. This article focuses on describing the methodologies (Section 2), the data (Section 3), the product design (Section 4), the known limitations (Section 5), and the roadmap for the product continuity (Section 6). The access to the source code called 'PYALUS' (Python code for ALbedo retrieval Using Satellite data) is provided in Section 7. Perspectives are discussed in Section 8.

## 2. Algorithm Description

The products presented in this article all rely on a common methodology that is slightly adapted for the specificities of each programme/product. This Section 2 presents the overall
methodology used for all the NRT and CDR products. The only differences between CDR and NRT products are the timeliness, the accuracy of the input data (forecast fields from climate models for NRT, and reanalysis for CDR), and the version of the scientific algorithm (current version for NRT, and the last version available which includes all the scientific improvements done before the processing date for the CDR).

## 2.1. Method Overview

The operational processing scheme of the albedo calculation is depicted in the flow chart of Figure 1. It encompasses four successive steps (one of them, step 2, concerns only the generation of long-term time series potentially combining different optical sensors):

Step 1. Atmospheric Correction: In a first step, the reflectances measured by the sensor at the Top of the Atmosphere (TOA) are used to estimate the reflectances at the Top of the Canopy (TOC) by applying the Simplified Method for the Atmospheric Correction (called SMAC) [36] of satellite measurements in the solar spectrum. This atmospheric correction process is described in Section 2.2.

Step 2. Harmonisation (optional): This second step aims to harmonise data that are spectrally heterogeneous because they have been acquired by different sensors. The spectral harmonisation step is an optional step that allows the calculation of spectral albedos at fixed wavelengths (corresponding to a chosen reference sensor) from different sensors having different spectral characteristics (only used in C3S). First, a reference sensor is chosen (i.e., the four bands of SPOT/VGT2). The TOC reflectances from each given sensor are then harmonised into reflectances that would have been observed by the reference sensor on each of its bands. This method is further detailed in Section 2.3.

Step 3. BRDF Inversion: In the third step, the measured TOC reflectances are used to fit the coefficients of a semi-empirical kernel-based reflectance model. These coefficients allow to rebuild the complete angular dependency of the bi-directional reflectance distribution function (BRDF). More information is presented in Section 2.4.

Step 4. Albedo Computation: This step is composed of two main processes. First, the spectral albedo values, which are associated to the instrument channels, are determined by angular integration of the bi-directional reflectance factors. Second, the narrow-to-broadband conversion of albedos is performed. More information on this last step is given in Section 2.5.

## 2.2. Atmospheric Correction

The atmospheric correction is the process that allows estimating the TOC reflectances with their associated uncertainties from the TOA satellite measurements. In the first step of the algorithm, TOC reflectance estimates are derived by applying SMAC [36]. SMAC provides reflectance measurements $R_{j\beta}$ ($j = 1, \dots, n$) at different spectral channels $\beta$, and for $n$ slots belonging to a composite window (see Section 4.1) which have discrete values of the view $\theta_{vj}$ and solar zenith angles $\theta_{sj}$. SMAC is based on the 6S radiative transfer formulation of the satellite signal [37]. All the pertinent radiative quantities are parameterized as a function of auxiliary data which are input as part of the SMAC code:

* Gas content (mainly ozone and water vapour);

* Aerosol content and aerosol type;

* Molecular scattering mainly driven by the sea-level surface pressure and the surface elevation.

In the albedo algorithm, pixels are considered for atmospheric correction if they are not classified as "cloud-covered" or "cloud-contaminated" in the input cloud mask (Figure 1). This atmospheric correction is performed separately for each pixel and each band.

In the C3S, Météo France is not in charge of the atmospheric correction processing, and TOC reflectances are pre-computed by HYGEOS partner from TOA radiances prepared by the German Aerospace Center (DLR). HYGEOS used the SMAC model (same radiative transfer model used by Météo France for the atmospheric correction in the LSA SAF).
<figure><img src="image_5.png" /><figcaption>Figure 1. Flow diagram of the surface albedo retrieval algorithm. The spectral harmonisation (step 2) is optional. The boxes at the right side are the inputs: i.e., aerosol optical depth (AOD), water vapour (WV), ozone (O3), digital elevation model (DEM) etc. The boxes in grey colour are the outputs of each module (box bordered by dashed lines).</figcaption></figure>

## 2.3. Spectral Harmonisation

This spectral harmonisation step is an optional step, today, performed within the C3S, only, to homogenise the NOAA-X/AVHRR, SPOT/VGT-1 & VGT-2 and PROBA-V time series from 1981 to 2015. The reference sensor here is chosen to be SPOT/VGT-2. During the full period, the final spectral albedo products are provided with SPOT/VGT2 spectral band characteristics. Spectral harmonisation consists of estimating equivalent reflectance $R_{\beta}^{(\text{sensor,harm})}$ that would have been observed by a VGT-2 sensor on band $\beta$, by combining the reflectances $R_{\beta'}^{(\text{sensor})}$ measured by each sensor at their own spectral bands $\beta'$. Equation (1) exhibits the model applied here for the harmonisation step of each sensor.

$$R_{\beta}^{(\text{sensor,harm})} = \alpha_{\rightarrow\beta}^{(\text{sensor})} + \sum_{\beta'} \alpha_{\beta'\rightarrow\beta}^{(\text{sensor})} \cdot R_{\beta'}^{(\text{sensor})} + \varepsilon \quad (1)$$

where $\alpha_{\rightarrow\beta}^{(\text{sensor})}$ is the constant coefficient (depending on $\beta$) of the linear model for the sensor to harmonise into the VGT-2 band $\beta$, $\alpha_{\beta'\rightarrow\beta}^{(\text{sensor})}$ is the weight (depending on $\beta$ and $\beta'$) of the band $\beta'$ in the linear model for the sensor to harmonise into the VGT-2 band $\beta$. In addition, $\varepsilon$ represents the Gaussian noise of the model but it is not a coefficient as such. For the cases where a TOC reflectance $R_{\beta'}^{(\text{sensor})}$ is missing for one or several bands $\beta'$, the corresponding harmonised reflectances $R_{\beta}^{(\text{sensor,harm})}$ are set to the missing value.

More details on the calculation method and the values of these coefficients in the C3S are given in Appendix A. The linear model coefficients, also called the harmonisation coefficients, $\alpha_{\rightarrow\beta}^{(\text{sensor})}$ and $\alpha_{\beta'\rightarrow\beta}^{(\text{sensor})}$ are provided in Appendix A.
## 2.4. BRDF Inversion

### 2.4.1. BRDF Models

In cloud-free conditions, satellite observations, after atmospheric correction (Section 2.2) and optionally spectral harmonisation (Section 2.3), deliver TOC reflectance measurements for the illumination and viewing geometries that are observed by the sensor. From these observations, we are able to infer an estimate of the complete angular dependence of the BRDF, $R_β(θ_s, θ_v, φ)$. This latter is estimated in the spectral channel, $β$ of the measuring sensor, and for all angular geometries, even for those far from the illumination and viewing geometries acquired by the sensor capabilities (that means for all sun $θ_s$ and view $θ_v$ angles comprised between 0 and 90°, and for all relative azimuth $φ$ angles comprised between 0 and 180°). The surface directionality estimation consists in determining, by mathematical inversion, the coefficients of a semi-empirical kernel-based reflectance model, based on TOC reflectances (or in the case of prior harmonisation, the TOC reflectances harmonised to another reference sensor). The system to invert, composed by a linear combination of $f_i$ kernel functions, is defined as follows:

$$R_β(θ_s, θ_v, φ) = k_β^T f(θ_s, θ_v, φ) \quad (2)$$

with $k_β = (k_{0β}, k_{1β}, k_{2β}, ...)^T$ and $f = (f_0, f_1, f_2, ...)^T$ representing vectors formed by the retrieved model parameters, $k_{iβ}$, and the kernel functions, $f_i$, respectively (with $i$ numbers of kernel functions).

As mentioned in Section 1, $f_i$ are semi-physical (or semi-empirical) kernel functions that simulate the brightness of the surface as series of contributions representing isotropic, volumetric, geometric, and other effects. $k_{iβ}$ parameters are the weights associated to these kernel functions. The objective in fitting the remote sensing observations that are observed by the satellite with these series of weighted kernel function is to address anisotropy behaviours for the different surface types (crop, forest, etc.). Then, this fitting also gives access to anisotropy of these surfaces for all the non-observed (or non-illuminated) configurations.

All products for which Météo France has the responsibility are listed in Table 1. All of these (except for historical reason MSG-based products) are derived using the RossThick-LiSparse-Reciprocal model [28] which is composed of three kernel functions representing the angular distribution related to isotropic, geometric and volumetric surface scattering processes. This RossThick-LiSparse-Reciprocal model, referred to as the RTLS model in the rest of this article, was defined as a result of the combination of a volumetric kernel (RossThick model) and a geometric kernel (LiSparse-R model). The two kernel functions of the RTLS model are represented in Figure 2. The determination of the kernel coefficients that are associated to these functions is performed using a linear least-square fit as described in the next section.

### 2.4.2. Least Square Solution

Figure 3 presents an overview of the variability of the geometric conditions for the input reflectances considered for the BRDF inversion of different sensors. A linear kernel-based directional reflectance model in the form of Equation (2) is applied separately to each spectral $β$ band. In the following, the index $β$ referring to the channel is omitted for simplicity.
<table><caption>**Table 1.** List of the operational albedo products developed by the Météo France research laboratory. LSA SAF (Satellite Application Facility on Land Surface Analysis) albedo (referred as LSA in the Table) are disseminated in NRT (near real-time) (3 h after the last observation at 00UTC). The timeliness will be reduced in 2021 to 1 h. C3S (Copernicus Climate Change Service) (referenced as C3S in the Table) are reprocessing. CGLS (Copernicus Global Land Service) product, which were under responsibility of Météo France and not still maintained, are in grey colour. PROBA-V product from CGLS was reused for the generation of C3S-V0 ICDR Interim Climate Data Record). In 2021, the C3S-V2 and LSA-107 and 108 will complete this portfolio.</caption><thead><tr><th>Product</th><th>Sensor (Coverage)</th><th>Product ID (Type)</th><th>Production Frequency</th><th>Composite Window</th><th>Temporal Characteristic Time Scale</th><th>Spatial Scale</th><th>Composition Method</th><th>BRDF Model</th><th>Atmospheric Correction</th><th>Documentation and Data Link Accesses</th><th>Temporal Coverage</th><th>Product Continuity</th></tr></thead><tbody><tr><td>Albedo (CDR)</td><td>NOAA-X/AVHRR (global)</td><td>C3S-V1</td><td>10-day</td><td>20-day</td><td>10-day</td><td>1-4 km</td><td>recursive</td><td>Li-Sparse Reciprocal</td><td>Hygeos</td><td>C3S-V1 access</td><td>1981-2015</td><td>C3S-V2</td></tr><tr><td>Albedo (CDR)</td><td>SPOT/VGT (global)</td><td>C3S-V1</td><td>10-day</td><td>20-day</td><td>10-day</td><td>1 km</td><td>recursive</td><td>Li-Sparse Reciprocal</td><td>Hygeos</td><td>C3S-V1 access</td><td>1998-2014</td><td>C3S-V2</td></tr><tr><td>Albedo (ICDR)</td><td>PROBA-V (global)</td><td>C3S-V0</td><td>10-day</td><td>30-day</td><td>20-day</td><td>1 km</td><td>recursive</td><td>Roujean</td><td>VITO</td><td>C3S-V0 access</td><td>2014+</td><td>C3S-V2</td></tr><tr><td>Albedo (NRT)</td><td>SPOT/VGT (global)</td><td>VGP-P</td><td>10-day</td><td>30-day</td><td>20-day</td><td>1 km</td><td>recursive</td><td>Roujean</td><td>VITO</td><td>CGLS access</td><td>1998-2014</td><td></td></tr><tr><td>Albedo (NRT)</td><td>PROBA-V (global)</td><td>PROBA-V L2A</td><td>10-day</td><td>30-day</td><td>20-day</td><td>1 km</td><td>recursive</td><td>Roujean</td><td>VITO</td><td>CGLS access</td><td>2014-2020</td><td></td></tr><tr><td>Albedo (NRT)</td><td>MSG/SEVIRI (Africa, Europe, South America)</td><td>LSA-101</td><td>1-day</td><td>1-day</td><td>5-day</td><td>SEVIRI grid</td><td>recursive</td><td>Roujean</td><td>Météo France</td><td>LSA-101 access</td><td>2005+</td><td>LSA-107 (MTG/FCI)</td></tr><tr><td>Albedo (NRT)</td><td>MSG/SEVIRI (Africa, Europe, South America)</td><td>LSA-102</td><td>10-day</td><td>30-day</td><td>30-day</td><td>SEVIRI grid</td><td>classic</td><td>Roujean</td><td>Météo France</td><td>LSA-102 access</td><td>2009+</td><td>LSA-108 (MTG/FCI)</td></tr><tr><td>Albedo (CDR)</td><td>MSG/SEVIRI (Africa, Europe, South America)</td><td>LSA-150</td><td>10-day</td><td>30-day</td><td>30-day</td><td>SEVIRI grid</td><td>classic</td><td>Roujean</td><td>Météo France</td><td>LSA-150 access</td><td>2005-2015</td><td></td></tr><tr><td>Albedo (NRT)</td><td>Metop/AVHRR (global)</td><td>LSA-103</td><td>10-day</td><td>20-day</td><td>10-day</td><td>1 km</td><td>recursive</td><td>Li-Sparse Reciprocal</td><td>Météo France</td><td>LSA-103 access</td><td>2015+</td><td>LSA-110, 111 (Metop-SG/METimage, 3MI)</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Angular dependence of the geometric (a) and volumetric (b) scattering kernels of the RTLS (RossThick-LiSparse-Reciprocal) model. Both graphs correspond to the principal plane (i.e., conditions with a relative azimuth equal to 0° or 180°). Negative zenith angle values correspond to the back-scattering direction (relative azimuth angle = 0°) and positive zenith angle values to the forward scattering direction (180°).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Illumination (a) and observation (b) geometries corresponding to a geographical location of (47°47′ N, 10°37′ E) and an observation period between the days of year 150 and 170. The relative azimuth angle is identical for the two graphs. The convention was chosen such that the top of the graphs (=0°) corresponds to the back-scattering regime. The colours of the dots denote observations taken by different sensors as follows: Red: MSG/SEVIRI, Green: Metop/AVHRR, Blue: NOAA-16/AVHRR.</figcaption></figure>

The available TOC reflectance measurements are combined to provide a system of linear equations following the matrix form of Equation (2). In general, the number of available observations is larger than the number of unknown parameters, and because neither observations and models are perfect and in perfect agreement with each other, no exact solution for the kernel coefficient *k* exists. The observed reflectances are affected by measurement errors, and BRDF models are built based on hypotheses and approximations.
Additionally, it is convenient to search for the best solution in a more statistical sense and quantify the uncertainties of the retrieved parameter estimates. In this case, it may turn out that a considerably larger number of observations than parameters is required in order to reasonably constrain the parameter values.

The uncertainties of the individual reflectance factor measurements $R_j$ are quantified by means of weight factors $w_j$, which are related to the inverse of the standard "1-sigma" uncertainty estimates $\sigma[R_j]$ (with $j$ the index for the cloud-free observation of $n$ elements). We introduce the scaled reflectance vector $b$ with the elements $b_j = R_j w_j$ and the "design matrix" $A$ with the elements $A_{ji} = f_{ji} w_j$ (see e.g., [38] with $i$ being the index for the kernel function number $i$ of usually three elements). The linear least squares solution to the inversion problem in Equation (2) can be found by solving the following "normal equations" for the parameters $k$.

$$ (A^T A) k = A^T b \qquad (3) $$

The uncertainty covariance matrix of the retrieved model parameters is given by

$$ C_k = (A^T A)^{-1} \qquad (4) $$

The diagonal elements $C_{ii}$ of this matrix represent the variance $\sigma^2[k_i]$ of the respective parameters $k_i$. If the matrix $A^T A$ is not singular, the solution can be found by multiplying Equation (3) "from the left" by the covariance matrix $C_k$.

### 2.4.3. Addition of a Priori Constraints Using a Recursive Method

A sufficient number of observations is needed for a robust inversion of the BRDF system (at least a minimum of three uncorrelated observations in the case when three kernel functions are used). Additionally, the BRDF inversion is usually computed every day for geostationary data, and every 10 days for polar satellites in order to collect enough observations over a composite period of several days (usually from 20 to 30 for polar orbiting satellites). In addition to the TOC reflectances collected over the composite period, it was decided to use the previous BRDF inversion through the use of a Kalman Filter (KF) method in order to better constrain the inversion process by considering older directional state of the surface. This has also the advantage of reducing the computational cost compared to a classic inversion scheme that would use more aged observations. Observations of several days before the beginning date of the composite window (i.e., >21 for a use of a 20-day composite window) are considered through a single a priori vector instead of collecting in a heavy matrix an additional number of observations and processing the mathematical inversion with all of these data. The KF method is also called recursive method, in contrast to a classical one with no a priori information and no weighted observations according to their ages in the inversion process. It affords better temporal coherences and spatial completeness of the estimates in time. The advantage to add constraints on the parameters themselves in the inversion of the linear model was already shown by Li et al. [39], Hagolle et al. [40], Geiger et al. [29], Carrer et al. [23,30,41] and a related approach was also adopted by Pokrovsky et al. [42]. In the following, the important steps of this recursive method are summarised. Independent and uncorrelated a priori information on the parameters is expressed in terms of the first and second moments (average and standard deviation, respectively) of their a priori (ap) probability distribution function, that is, an estimate of the form:

$$ k_i = k_{ap} \pm \sigma_{ap} [k_i] \qquad (5) $$

To simplify the notation, we consider here an example with a BRDF model with three kernel functions ($m = 3$), and an additional constraint for the three parameters $k_0$, $k_1$ and $k_2$.
In this case, adding the constraint of Equation (5) to the system of Equation (3) corresponds to extending the (n, m)-matrix A to the (n + m, m)-matrix

$$A' = \begin{vmatrix} & A & \\ \sigma_{ap}^{-1}[k_0] & 0 & 0 \\ 0 & \sigma_{ap}^{-1}[k_1] & 0 \\ 0 & 0 & \sigma_{ap}^{-1}[k_2] \end{vmatrix} \quad (6)$$

and extending the vector b to $b' = (b_1, \dots, b_n, k_{0ap}, \sigma_{ap}^{-1}[k_0], k_{1ap}, \sigma_{ap}^{-1}[k_1], k_{2ap}, \sigma_{ap}^{-1}[k_2]^T)$. The linear least squares solution with a priori information is then obtained in the same way as above by solving the normal equations where both A and b are replaced respectively by $A'$ and $b'$. The use of an a priori information triggers the issue of the initialisation of the run at the first-time step. The next section presents the approach used to generate a meaningful a priori for the first run.

Taking into account a multivariate Gaussian a priori probability distribution function for the parameters, quantified by its first and second moments, corresponds to re-writing Equations (3) and (4) in the form

$$(A^T A + C_{ap}^{-1})k = A^T b + C_{ap}^{-1} k_{ap} \quad (7)$$

and

$$C_k = (A^T A + C_{ap}^{-1})^{-1} \quad (8)$$

with $k_{ap} = (k_{0ap}, \dots, k_{m-1ap})^T$ and $C_{ap}$ the covariance matrix. For uncorrelated a priori information on the parameters, $C_{ap} = \text{diag}(\sigma_{ap}^2[k_0], \dots, \sigma_{ap}^2[k_{m-1}])$ is diagonal. The determination of $k_{ap}$ and $C_{ap}$ is explained in [30], and also briefly in Section 2.4.4.

### 2.4.4. Initialisation—Determination of the First a Priori Information

As explained in the previous section, the information coming from the satellites is combined over a composite period in order to reduce the sensitivity of the resulting estimates to reflectance outliers. The processing scheme based on a KF method allows the recursive procedure in time. At each execution of the algorithm (typically 10 days for polar orbiting sensors), the previous parameter estimate and the corresponding uncertainty measure (called $k_{in}$ and $C_{in}$) are used. These previous estimates are then used in the following way as a priori information for the linear model inversion specified in Equations (7) and (8):

$$ \begin{aligned} k_{ap} &= k_{in} \\ C_{ap} &= C_k^{in} \delta, \end{aligned} \quad (9) $$

where $\delta$ is an inflation factor bigger than 1 (more detail in [30]).

In order to initialise the KF process with a first meaningful a priori BRDF model, the following "spin up" procedure is run for each sensor (see Figure 4):

* Step 1—Spin up run—Run is performed for one year. The BRDF estimates of the last day of the spin up period is later used in Step 2 as first guess (a priori BRDF in Equation (7)). After one year, the KF has lost memory of its initial state: the lack of initial BRDF model has no more impact on the output product. In most cases, a shorter period (around 3 months) is sufficient to initialize the KF depending on the cloudiness, but one year has been chosen to take into account the vegetation cycle.

* Step 2—Actual run—Using the latest model from step 1, consider it as the initial a priori BRDF model and generate the product for the full period with Kalman filtering enabled.

In addition to the consideration of a priori information, additional constraint of regularisation was added for pure numerical reasons. This is described hereafter.
<figure><img src="image_8.png" /><figcaption>Figure 4. Spin up and actual run.</figcaption></figure>

### 2.4.5. Regularisation

It is well-known that matrix inversion can lead to numerical problems (e.g., instabilities when eigen values are close to 0 or have the same values). Numerical instabilities can also occur in the inversion process when the information content of the available observations is very scarce. This is why regularisation terms with small means and high variances are added in the BRDF inversion equations. In practice, regularisation consists in adding extra Gaussian signals:

$$ 
\begin{aligned}
k_0 &= k_{0\text{reg}} \pm \sigma_{\text{reg}}[k_0] \\
k_1 &= k_{1\text{reg}} \pm \sigma_{\text{reg}}[k_1] \\
k_2 &= k_{2\text{reg}} \pm \sigma_{\text{reg}}[k_2]
\end{aligned} 
\qquad (10) 
$$

According to the KF theory, the same methodology as detailed in Equation (6) can be applied to account for these regularisation (reg) terms. Values of $k_{reg}$ terms are equal to typical mean values for each kernel parameter. The most important thing to notice is that their associated uncertainties are fixed to very large values which ensure low impact of this information on the retrieval (for example, we have chosen the constraints of $k_1 = 0.03 +/ - 0.05$). The goal is here only to avoid numerical instabilities and the standard deviations are sufficiently large to avoid a noticeable prejudice in the inversion result. Moreover, if the a priori information from a previous retrieval is available with an associated uncertainty not too important, this regularization has clearly no (or negligible) impact.

### 2.5. Albedo Computation

Once the BRDF is determined, the albedo computation can start. The kernel approach offers a description of the angular dependence of the reflectance factor. It is applied to each instrument channel separately and provides no information on the spectral behaviour outside of the available narrow bands.

#### 2.5.1. Angular Integration

After inversion of Equations (7) and (8), and the calculation of the kernel parameters for each wavelength $\beta$, the directional-hemispherical 'dh' and bi-hemispherical 'bh' albedos are estimated from the integrals of the kernel functions $f = (f_0, f_1, f_2)^T$ which are noted respectively $I^{dh}(\theta_s) = (I_0^{dh}(\theta_s), I_1^{dh}(\theta_s), I_2^{dh}(\theta_s))^T$ and $I^{bh} = (I_0^{bh}, I_1^{bh}, I_2^{bh})^T$. This angular integration is described in [29,30] by inserting $R_\beta$ (from Equation (2)) in the following integral equation:

$$ 
\begin{aligned}
a_{\beta}^{dh}(\theta_s) &= \frac{1}{\pi} \int_{2\pi} R_{\beta}(\theta_s, \theta_v, \varphi) \cos \theta_v d\Omega_v \\
a_{\beta}^{dh}(\theta_s) &= k_{\beta}^{T} I^{dh}(\theta_s) \quad \text{and} \quad a_{\beta}^{bh} = k_{\beta}^{T} I^{bh} \\
\sigma[a_{\beta}] &= \sqrt{I^{T} C_{k} I}
\end{aligned} 
\qquad (11) 
$$

The resulting narrowband albedo values serve as approximations for the spectral albedo at the central band with wavelength $\lambda_\beta$.

#### 2.5.2. Narrow-to-Broadband Conversion

Broadband albedo is defined as the integral of spectral albedo over a certain wavelength interval weighted by the spectral irradiance. Since the integral can be approximated as a weighted sum at discrete values of the integration variable, broadband albedo may be
expressed as a linear combination of the spectral (or rather narrowband) albedo values in the available instrument channels.

The broadband albedo estimates for a given target interval $\gamma = [\lambda_1, \lambda_2]$ are derived from the spectral quantities by applying a linear transformation of the form

$$a_{\gamma} = c_{0\gamma} + \sum_{\beta} c_{\beta\gamma} a_{\beta} \qquad (12)$$

with coefficients $c_{0\gamma}$ and $c_{\beta\gamma}$ summarised in Algorithm Theoretical Basis Documents (ATBD) available on the two project websites for each sensor. Three different broadband albedo intervals are considered: the total shortwave range from 0.3 µm to 4 µm (BB), the visible wavelength range from 0.4 µm to 0.7 µm (VIS), as well as the near infrared range from 0.7 µm to 4 µm (NIR). Coefficients are different for snow-free observations and in case of snow/ice. For the calculation of these narrow to broadband conversion coefficients, a linear regression based on radiative transfer simulations is used. An extensive dataset of synthetic spectral reflectances of soil and vegetated surfaces was generated for different surface types by using the Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER) spectral library [43] and the SAIL radiative transfer model [44]. After calculating the narrowband albedo values in the instrument’s spectral bands and the broadband albedo values in the ranges of interest, the corresponding linear transformation coefficients can be determined.

In the evolutions of the C3S and LSA SAF algorithm, the linear transformation depicted by Equation (12) will be replaced by polynomial series of NDVI (Normalized Difference Vegetation Index) in order to better consider the type of surface. Performances of the narrow to broadband conversion should be particularly improved over desert area, and forests (tail end of the distribution).

The weighting with the spectral irradiance in the definition of the broadband albedo introduces a dependence on the atmospheric conditions since the spectral properties of the downward irradiance are different under direct and diffuse conditions. For the time being, this difference is not taken into account in the generation of the input datasets for the regression analysis and the same narrow to broadband conversion relations are applied for the directional-hemispherical albedo $a^{\text{dh}}(\theta_{\text{ref}})$, irrespective of the reference illumination angle, and for the bi-hemispherical albedo variant.

## 2.6. Extra-Filters—Impact of Clouds, Cloud Shadows, and Eclipses

Only clear-sky observations of the surface are interesting for the estimation of the land surface albedo in the optical visible domain (0.3 µm–4.0 µm). In order to select these clear-sky observations of the surface, external cloud mask (CMa) products are used. These CMa products provide information about the confidence of the cloud detection through associated quality flags. In the present algorithm, cloudy observations are discarded. In addition, the observations for which the respective flag of the cloud mask product CMa indicates a bad quality are penalized in the weighting scheme by multiplying the reflectance uncertainty estimates ($\sigma^2[R]$), see Section 2.4.2) by a factor 10. The same approach is adopted for the pixels that might be contaminated by cloud shadows, when positioned next to cloudy pixels and depending on the solar azimuth direction. Cloud shadow is estimated by using the sun geometry for expanding the cloud mask as described in LSA SAF albedo ATBDs (available on https://www.eumetsat.int/lsa-saf). In this way, those potentially affected observations are only significant in the inversion process if no “reliable” observations are available at all. Moreover, a solar eclipse calendar database is used to systematically remove the global file corrupted by a solar eclipse. Shadow of clouds and solar eclipse induce strong and rapid decreases of surface brightness which are not linked to intrinsic changes of land surface properties.
## 3. Data

The presented method has been applied with some differences (which will be detailed) to different sensor data in the framework of LSA SAF and C3S programmes. The algorithm requires in input: (i) TOA radiances at several spectral bands, (ii) radiometric uncertainties for each band, (iii) the sensor and sun view and azimuth angles, (iv) the image acquisition time, (v) longitude and latitude fields, (vi) a quality flag including land/water mask, (vii) a binary cloud mask and a snow/ice mask, (viii) and several auxiliary data (Section 3.1) comprising atmospheric components as ozone, water vapour, pressure and aerosols. These input data and example of output products are described in the following for each Earth Observing system.

### 3.1. Input: Auxiliary Data

#### 3.1.1. Digital Elevation Model (DEM)

The C3S and the LSA SAF albedo processing chains both rely on a DEM derived from the USGS GTOPO30 (Global 30 Arc-Second Elevation) dataset (https://doi.org/10.5066/f7DF6PQS, last consulted on 25 March 2020). This DEM is reprojected on each sensor grid using a bilinear interpolation so that the information is available for each pixel to process.

#### 3.1.2. Atmospheric Parameters

In C3S programme, we used the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2; Ronald Gelaro et al., 2017) as data source for hourly ozone content, water vapour, pressure, and aerosol products (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/, last accessed on 25 March 2020). It is a NASA atmospheric reanalysis for the satellite era using the Goddard Earth Observing System Model, Version 5 (GEOS-5) with its Atmospheric Data Assimilation System (ADAS), version 5.12.4. The datasets are available from 1st January 1980 until present. They can be downloaded from several sites, for example (https://disc.sci.gsfc.nasa.gov/daac-bin/FTPSubset2.pl). MERRA-2 was chosen as it was the only global scale re-analysis available from the 1980s.

In the EUMETSAT LSA SAF, the objective is to distribute NRT products within strict timeliness. This timeliness (target delay) is comprised between 1 and 3 h after the last satellite observations in order to deliver level 3 products to users in near real-time. The same list of atmospheric variables listed above from MERRA-2 are obtained from the ECMWF numerical weather prediction model except for aerosol optical thickness. For aerosols, a climatology is used from the Copernicus Atmosphere Monitoring Service (CAMS) over a 10-year period from 2003 to 2012. The first tercile is used instead of the mean value in order to avoid overcorrections.

Atmospheric gases apart from ozone and water vapour are also considered in the atmospheric correction but with a static concentration in the atmospheric column.

#### 3.1.3. SMAC Coefficients

These coefficients constitute the core of the SMAC algorithm. They have been calculated for all instruments and for several aerosol models. The source for these coefficients is the CNES/Cesbio laboratory. They can be found here (http://tully.ups-tlse.fr/olivier/smac-python/tree/master/COEFS, last accessed 6 August 2020).

## 3.2. Radiance Inputs and Albedo Outputs

### 3.2.1. NOAA-X/AVHRR

The data described here are used to generate the version 1 of the surface albedo product in the framework of the C3S programme (C3S-V1). Data from the AVHRR-2 sensor on board the NOAA-7, 9, 11, 14, and AVHRR-3 sensor on board the NOAA-16 and NOAA-17 platforms were exploited. The AVHRR TOA radiances are in the red (RED) band (0.580–0.68 µm for both AVHRR sensors), the near-infrared (NIR) band (0.725–1.10 µm for AVHRR-2, 0.725–1.00 µm for AVHRR-3), and, only for AVHRR-3 on board the NOAA-16 and NOAA-17, the mid-infrared (MIR) band (1.58–1.64 µm).
The NOAA-X/AVHRR revisit time allowed for the coverage of the full globe is approximately every 12 h. It is also important to underline that at the time of the generation of the v1 collection of C3S surface albedos that is presented in this article, no information on the presence of snow/ice mask was available; this induced limitations over pixels covered by snow/ice (discussed in Section 5.2.2). Snow and ice presence will, however, be available for the generation of the C3S v2 collection. Another consideration is the change in the acquisition geometry between each platform, the difference of spatial resolution and the changes of spectral responses. Figure 5 also shows the evolution of the spectral response functions from each AVHRR sensor and the spectral response of VGT-2, which is the reference used for the spectral harmonisation of all instruments (see Section 2.3).

<figure><img src="image_9.png" /><figcaption>Figure 5. *Cont.*</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>(g) Spectral response for VGT2 sensor on a (300,2500) nm range</figcaption></figure>

**Figure 5.** Spectral response on the (300 nm, 2500 nm) interval for NOAA-X/AVHRR and VGT-2 sensors. Only NOAA: (a)—7, (b)—9, (c)—11, (d)—14, (e)—16, and (f)—17 are considered, and (g) VGT2 in C3S.

Based on these input data and the algorithm detailed in Sections 2.1–2.6, NOAA-X/AVHRR C3S-V1 albedos were generated. In the next sections, examples of retrieval will be shown for 1 July 2005 or 1 July 2015 (for the most recent sensors). Figure 6 shows an example of retrieval of NOAA-X/AVHRR C3S-V1 albedo for this date.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Total shortwave directional-hemispherical albedo (BB-DH) for 1 July 2005 based on NOAA-X/AVHRR (C3S-V1).</figcaption></figure>

### 3.2.2. SPOT/VGT

The data described here are used to generate the version 1 surface albedo in the framework of the C3S (C3S-V1). The data used are the SPOT/VGT FCDR (Fundamental Climate Data Record) dataset. They are based upon the Collection 3 of SPOT/VGT TOA reflectance archive. Radiometric calibration and inter-calibration were performed between SPOT4/VGT-1 and SPOT5/VGT-2 [45]. Assessment of the absolute uncertainty in radiometric calibration over the Libya-4 calibration site using the Optical Sensor Calibration with Simulated Radiance (OSCAR) desert calibration approach [46] showed relative small differences between VGT-1 and VGT-2 at TOA level for Collection 3, with the largest bias for the red band (around 1%, depending on the period under consideration). The SPOT/VGT TOA reflectances are at the nominal resolution of approximately 1 km × 1 km (or 1/112°) and the pixel geolocation is provided for the pixel centre. The VGT-1 and VGT-2 TOA used data are bands B0 (blue, 0.43–0.47 µm), B2 (red, 0.61–0.68 µm), B3 (NIR, 0.78–0.89 µm) and SWIR (Shortwave Infrared, 1.58–1.75 µm)
Based on these input data and the algorithm detailed in Sections 2.1–2.6, SPOT/VGT C3S-V1 albedos were generated. Figure 7 shows an example of retrieval of 1 July 2005. This illustration was chosen for the same date than NOAA-X/AVHRR (C3S-V1) in Section 3.2.1.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Total shortwave directional-hemispherical albedo (BB-DH) for 1 July 2005 based on SPOT/VGT (C3S-V1).</figcaption></figure>

### 3.2.3. Metop/AVHRR-3

The Metop/AVHRR-3 data are used to generate the EPS (EUMETSAT Polar System) Ten-daily surface ALbedo (called ETAL) within the framework of the LSA SAF. The series of Metop satellites is developed by the European Space Agency (ESA) and operated by the EU-METSAT. It is composed of the satellite Metop-A which was launched on 19 October 2006), Metop-B launched on 17 September 2012 and Metop-C launched on 7 November 2018, are both flying the AVHRR-3 sensor. This sensor has similar characteristics as the one flying on board NOAA-16 and NOAA-17, (described earlier in the NOAA section). The Metop platforms are well stabilised and protected against orbital drift. The raw segments are mapped (nearest neighbour) to the geographic longitude/latitude system with the same framing/resolution (1°/112) as the one used for SPOT/VGT. The TOA input data used correspond to the AVHRR-3 bands; same as those of NOAA-16/17.

Based on these input data and the algorithm detailed in Sections 2.1–2.6, Metop/AVHRR-3 ETAL albedos were generated. Figure 8 shows an example of retrieval of 1 July 2015.

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Total shortwave directional-hemispherical albedo (BB-DH) for 1 July 2015 based on Metop/AVHRR-3 (ETAL).</figcaption></figure>
### 3.2.4. MSG/SEVIRI

The MSG/SEVIRI data are used to generate the surface albedo products in the framework of the LSA SAF. The SEVIRI radiometer embarked on the spin-stabilised MSG platform encompasses unique spectral characteristics and accuracy, with a 3 km resolution (sampling distance) at nadir (1 km for the high-resolution visible channel), an imaging-repeat cycle of 15 min and 12 spectral channels [47]. The same list of input data (water vapour, ozone content, etc.) as described earlier for the Metop/AVHRR sensor is used for calculation of MSG albedo products. The TOA radiances of SEVIRI exploited correspond to the red band (0.6 µm), NIR band (0.8 µm) and SWIR band (1.6 µm).

Based on these input data and the algorithm detailed in Sections 2.1–2.6, Meteosat Daily Albedos (called MDAL) were generated. Figure 9 shows an example of retrieval of the 1 July 2015.

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Total shortwave directional-hemispherical albedo (BB-DH) for 1 July 2015 based on MSG/SEVIRI (MDAL).</figcaption></figure>

### 3.2.5. PROBA-V

The PROBA-V satellite was launched on 6th May 2013 and was designed to bridge the gap in space-borne vegetation measurements between SPOT/VGT (March 1998–May 2014) and the Sentinel-3 satellites, with the first one launched in 2016, which will ensure the continuity of the SPOT/VGT mission. The optical design of PROBA-V/VGT-P sensor consists of three cameras measuring in 4 spectral bands comparable to the SPOT/VGT ones: blue (0.415–0.50 µm), red (0.58–0.77 µm), NIR (0.73–0.96 µm) and SWIR (1.48–1.76 µm). Each camera has two focal planes, one for the SWIR and one for the visible and near-infrared (VNIR) bands. These data will be used to generate the C3S-V2 surface albedo products in the framework of the C3S. However, PROBA-V surface albedo in C3S was, in a first step, brokered from PROBA-V surface albedo Collection 1 from Copernicus/CGLS project [48], which was generated in near real-time from sensor data since 2014 with a 1 km × 1 km spatial resolution. The CGLS PROBA-V albedos were contracted by the Joint Research Center (JRC) under the scientific responsibility of Météo France (https://land.copernicus.eu/global/themes/energy). The algorithm was originally developed in the context of the FP5/CYCLOPES project [49] for an application to SPOT/VEGETATION data. It has been then adapted [48], and is used since 2013 to produce the surface albedo products of the CGLS based on VEGETATION and PROBA-V sensor data. These brokered data offer an interim solution to the user in order to provide continuity after the end-of-life of SPOT/VGT mission. The Copernicus/CGLS brokered data plus several additional fields constitute the so-called PROBA-V ICDR (Interim Climate Data Record) surface albedo in the C3S programme. This collection is called C3S-V0. Only broadband albedos were made
available in Copernicus/CGLS. The C3S-V0 PROBA-V ICDR surface albedo brokered the CGLS PROBA-V broadband albedo, and additionally provides new estimates of spectral albedos with their associated uncertainties. The PROBA-V dataset was homogenised to SPOT-VGT2-like sensor which was chosen by the C3S as the reference sensor for all vegetation products (same protocol as in Sections 3.2.1 and 3.2.2). In that way, C3S-V1 CDR plus C3S-V0 ICDR have the ambition to offer continuous collection of surface albedo (spectral and broadband with associated uncertainties) from 1981 until today. In 2021, C3S-V2 will reprocess PROBA-V albedo as NOAA-X/AVHRR and SPOT/VGT series with the same version of algorithm. In the meantime, only C3S-V0 PROBA-V ICDR is available from 2015 until June 2020, end date of the operational PROBA-V mission. Despite this effort for the homogenisation, several inconsistencies exist including: C3S-V0 is based on a compositing window of 30 days (20 days for C3S-V1); and two collections are using two different BRDF models (Table 1). Based on these CGLS input data, C3S-V0 PROBA-V albedos were generated. Figure 10 shows an example of retrieval of 1 July 2015.

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Total shortwave directional-hemispherical albedo (BB-DH) for 1 July 2015 based on PROBA-V (C3S-V0).</figcaption></figure>

# 4. Product Design

## 4.1. EUMETSAT and C3S Albedos

Table 1 summarizes the characteristics of the different products Météo France has developed and which reached the operational status. ‘Operational’ means that the products have been quality assessed by independent reviewers, and EUMETSAT or EU respectively decided, following the recommendations of these reviewers, to disseminate the products to the community. The temporal characteristics of the different albedo products differ from one sensor to another, and also for a given sensor from one project to another.

In EUMETSAT/LSA SAF, a monthly-based surface albedo (LSA-102) was derived from MSG/SEVIRI with a classic composition method (not recursive) as described in [30]. This classic composition consists in a simple average of the daily albedo estimates—each albedo having the same weights over the 30-day composite window (if these days have exactly the same conditions of observations—otherwise, a weighting occurs according to their daily uncertainty estimation [30]). Excepting this, all the surface albedos developed by Météo France derived from geostationary (polar) satellites use observations collected over 1 day (20 days) and a characteristic time scale of 5 days (10 days) following a recursive strategy reminded in the following paragraph. Table 1 also reports characteristics of SPOT/VGT and PROBA-V albedos developed by Météo France in Copernicus/CGLS, and later used for ICDR PROBA-V C3S albedo (see Section 3.2.5).

The definition of the characteristic time scale $\delta$ in Equation (9) is recursively applied to the previous retrieval uncertainty in order to propagate in a recursive manner the last estimation of the BRDF parameters. Parameter $\delta$ is an inflation factor. The more we decide to have an estimate representative of the conditions of the production date, the less important is the value of $\delta$. However, higher is also the risk to have artefacts due to isolated
cases of poor-quality atmospheric corrections. This value of $δ$ depends on the choice of a variable $τ$ (described in [30] which represents the temporal characteristic time scale (see values in Table 1). For example, if a temporal characteristic time scale of 10 days is chosen, this gives, to observations collected over the composite window, the following weights: 50% for 10-day-old observations and only 10% for 20-day-old observations. However, an internal layer also provides the mean age of the observations used for the composite window (see Section 4.3). This age information can also be used to access more precisely the real age of the provided estimates. The values of the temporal characteristic time scale (which is moderating $δ$ inflation factor) and of the composite window (period for which data are collected for the inversion) are selected according to the periodicity of the observing system (geostationary or polar) and to the user requirements. User requirements usually request the best compromise between temporal resolution and sensitivity to remaining small-scale abnormal variations in the reflectances (which are often due to uncorrected atmospheric effects). For example, despite its non-global scale coverage, the highly frequent observations provided by geostationary satellites can be interesting over boreal regions in order to track the rapid melt of ice and snow [50]. They can also be used to detect very small shifts in the vegetation cycle (of about +2/+3 days per decade in the start of the growing season as discussed in [51,52].

## 4.2. Albedo Characteristics: Spectral and Temporal

Two quantities are retrieved at different wavelengths (spectral albedos) and at different spectral domains (broadband albedos) from each instrument:

* Bi-hemispherical ('white-sky' or 'WSA') albedo products that are representative of diffuse conditions of illumination (typically cloudy sky conditions).

* Directional-hemispherical ('black-sky' or 'BSA') albedo products that are representative of direct conditions of illumination (typically clear sky conditions with pure atmosphere). As the surface is usually non-Lambertian and has directional properties, the value of surface albedo is given at a reference angle (solar position at the local noon).

Quality flags and uncertainties are also delivered for each albedo estimate.

In the C3S programme (for V0 ICDR and V1 CDR), the spectral albedos have the native SPOT-5/VGT-2 spectral sensor characteristics whatever the sensor origin (NOAA-X/AVHRR, SPOT-4/VGT-1, or PROBA-V). The C3S-V0 and C3S-V1 surface albedo consist of two products:

* Spectral albedo for VGT-2 channels (B0, B2, B3 and MIR; Table 2);

* Broadband albedo for the visible (VIS) (0.4–0.7 µm), near-infrared (NIR) (0.7–4 µm) and the total shortwave (BB) (0.3–4 µm) spectral domains.

<table><caption>**Table 2.** VGT-2 spectral bands, full width at half maximum and centre wavelength in parentheses.</caption><thead><tr><td>Band</td><td>VGT-2 (µm)</td></tr></thead><tbody><tr><td>Blue (B0)</td><td>0.439–0.476 (0.458)</td></tr><tr><td>Red (B2)</td><td>0.616–0.690 (0.653)</td></tr><tr><td>NIR (B3)</td><td>0.783–0.892 (0.838)</td></tr><tr><td>SWIR (MIR)</td><td>1.584–1.685 (1.635)</td></tr></tbody></table>

Surface albedo products are derived from SPOT-VEGETATION data every 10 days using a compositing window of 20 days. The dates of production are the 10th, 20th and 30th of each month in C3S-V1. In C3S-V0, the date of the mean age of the observations is indicated as product name instead of the production date (CGLS heritage).

In the LSA SAF, all spectral albedos have the native spectral characteristics of the sensors (MSG/SEVIRI, and Metop/AVHRR; Table 3). No harmonisation to a reference sensor is performed. However, the same broadband albedos for VIS, NIR and BB are provided. The 10-day albedos are produced for the 5th, 15th and 25th of each month.
<table><caption>**Table 3.** MSG/SEVIRI and Metop/AVHRR-3 spectral bands, full width at half maximum and centre wavelength in parentheses.</caption><thead><tr><td>Band</td><td>SEVIRI (μm)</td><td>AVHRR-3 (μm)</td></tr></thead><tbody><tr><td>Blue (B0)</td><td>-</td><td>0.439–0.476 (0.458)</td></tr><tr><td>Red (B2)</td><td>0.56–0.71 (0.635)</td><td>0.616–0.690 (0.653)</td></tr><tr><td>NIR (B3)</td><td>0.74–0.88 (0.81)</td><td>0.783–0.892 (0.838)</td></tr><tr><td>SWIR (MIR)</td><td>1.50–1.78 (1.64)</td><td>1.584–1.685 (1.635)</td></tr></tbody></table>

### 4.3. Product Content

In C3S and LSA SAF multi-layer files are provided for broadband and spectral albedos (BH and DH). Broadband directional-hemispherical albedo files contain the following variables: AL_DH_VI (broadband directional-hemispherical albedo over the visible spectral range (0.4, 0.7 µm); AL_DH_VI_ERR (uncertainty on the broadband directional-hemispherical albedo over the visible spectral range); AL_DH_NI (broadband directional-hemispherical albedo over the near-infrared spectral range (0.7–4 µm); AL_DH_NI_ERR (uncertainty on the broadband directional-hemispherical albedo over the near-infrared spectral range); AL_DH_BB (broadband directional-hemispherical albedo over the total shortwave spectral range (0.3–4 µm); AL_DH_BB_ERR (uncertainty on the broadband directional-hemispherical albedo over the total shortwave spectral range); QFLAG (the quality flag of the product); NMOD (for C3S product only, the number of valid observations during the composite window that are used to calculate the surface albedo); and AGE (mean age, in number of days, of the observations that are used to produce the output surface albedo and belong to the composite period. For example, if all observations of the 20-day compositing window are valid, the mean age is 10).

For broadband bi-hemispherical albedo (AL_BH) product, the albedo variables and their uncertainties are named similarly, replacing AL_DH by AL_BH. Other variables (AGE, QFLAG and NMOD) are unique and valid for all the albedos (spectral and broadband, directional-hemispherical and bi-hemispherical). Spectral directional-hemispherical albedo (AL_SP_DH) files contain the spectral directional-hemispherical albedo for each band (VGT-2-like bands in C3S, or native sensor bands in LSA SAF) and their associated uncertainties. For spectral bi-hemispherical albedo (AL_SP_BH) product, the albedo variables and their uncertainties are named similarly, replacing AL_DH by AL_BH. AGE, QFLAG and NMOD are also included in the spectral albedo product files. The naming of these variables differs a little from one project to another. The user is advised to directly refer to each Product User Manual (PUM also called Product User Guide in C3S) for more information (see 11th column of Table 1 for the link access to this documentation).

## 5. Discussion and Known Issues

### 5.1. Differences between the LSA SAF and C3S Albedo Products

#### 5.1.1. Atmospheric Correction

The method used for atmospheric correction in the albedo algorithm is very similar among instruments and programmes, all making use of the SMAC code. The main difference lies in the atmospheric parameters used as input (listed in Section 3) to model the extinction processes from gases and aerosols. The most critical input parameter is the aerosol content (or AOD) which has often the largest impact on the clear-sky atmospheric transmittance.

#### 5.1.2. Spectral Harmonisation

In C3S, the TOC reflectances of each sensor are harmonised to reconstruct VGT-2-like TOC reflectances that would be observed by the four native VGT-2 spectral channels. The calculations of the BRDF parameters and spectral albedos from the NOAA-X/AVHRR and SPOT/VGT-1 sensors do not depend anymore on the native spectral characteristics. Instead, they are derived from the VGT-2-like reflectances reconstructed from each sensor.
Table 2 gives VGT-2 spectral bands and their widths. The goal of using this spectral harmonisation is to prepare the merging of TOC data from different sensors at the level of the BRDF inversion and to provide harmonised spectral albedos from the 1980s.

In EUMETSAT/LSA SAF, the spectral harmonisation step is not processed until today. The TOC reflectances at each β band are used in their native spectral sensitivity for the calculation of the BRDF parameters and spectral albedos at each β band before the estimation of the broadband albedos. The spectral albedos are depending on the spectral characteristics of the instruments (MSG/SEVIRI; Metop/AVHRR-3; see Table 3). In other words, LSA SAF does not derive a standardized set of spectral albedos, but instead adjusts a set of coefficients that directly convert the set of channel albedos of each sensor into broadband albedo values. Therefore, in contrast with C3S, an independent set of narrow-to-broadband coefficients is determined for each sensor within the LSA SAF. In C3S, a unique set of narrow-to-broadband coefficients is used but instead, several sets of spectral harmonisation coefficients to VGT-2-like are required for each sensor. The two approaches are equivalent since similar database and similar regression method are used to calculate these narrow-to-broadband, in one case, or spectral harmonisation coefficients, in the other (see Appendix A).

## 5.2. Known Issues and Limitations

The factors limiting the quality of the products are of several types.

### 5.2.1. Residual of Clouds and Subpixel Clouds

The impact of this issue is difficult to quantify. However, we can frequently observe outliers in the TOA radiances which find an explanation in a residual cloud contamination. The strategy we adopted through the use of the recursive method consists in giving less weight to observations possibly affected by clouds or cloud shadows. In case of geostationary observations, observations just before and after cloud detection are systematically considered as poor quality and potentially contaminated by clouds. This very restrictive protocol is permitted only for surface albedo retrievals based on geostationary observations which benefit from frequent scans of the surface.

### 5.2.2. Treatment of Snow Target Pixels

In the current version of the algorithm, the pixel is identified as snow-contaminated (resp. snow-free) when the majority of the observations for the composite period are tagged snow (resp. snow-free). The snow/ice information is an input data (usually embedded in the cloud mask products). If the pixel is snow-free (resp. snow-contaminated), all observations flagged with snow are discarded (and reciprocally) for the inversion of the BRDF model. Different NB-BB coefficients, depending on the snow status of the pixel, are then considered for the narrow to broadband conversions. The recursive method that used data for a composite period of several days (usually 20 days) can provide misleading information during the period of transition from snow to no snow (and reciprocally). As these transitions of the surface between snow and no snow are very often progressive in time, it is difficult to circumvent the issue. In addition, we point out an important limitation of C3S-V1 based on NOAA-X/AVHRR series. For NOAA-X/AVHRR series, snow/ice information was not available. Moreover, all the pixels are flagged as 'no snow'. This induces underestimations of albedo over snow/ice areas (see Figure 6 over Greenland). Snow characterisation will be improved in C3S-V2, and snow/ice information will be reported for the entire period of the collection.

### 5.2.3. Calibration, Radiometry, and Orbit Drift

Even if it is hard to make the distinction between normal and abnormal shifts in time of the satellite observations, we suspect series of calibration, radiometry, and orbit drift issues. For example, calibration of the oldest instrument (NOAA-X/AVHRR) used for C3S-V1 probably needs to be improved. EUMETSAT reported a bias of around 8% between
MSG/SEVIRI and MODIS red bands [53]. Hence, we cannot discard that LSA-101 and LSA-102 products are not impacted by these potential calibration issues (if we consider MODIS as a calibration reference). It is also well known that spectral band sensitivity of the sensors may have evolved in time with their age. This potential radiometric drift is, unfortunately, usually unknown (because we are not able to quantify it in time), and may largely impact the spectral harmonisation and narrow to broadband steps. Finally, orbit drifts usually occur for the end-of-life period of the satellites. During this period, orbit correction interventions are done in order to stabilise, for the longest possible period, the trajectory of the satellite. Even if we made the assumption that these orbit drifts are perfectly documented, the change of angular geometry in the viewing/illumination angles undeniably impacts the BRDF retrieval. The sampling of the geometry configurations observing the surface is changing during the orbit drifts, and angular characterisation is very poor for past and current instruments on board polar observing systems which could create slight changes of the surface albedo estimations. This is maybe less the case for PROBA-V which has three cameras, and will be no longer the case with 3MI instrument on board Metop-SG (LSA-111) which will offer multi-angular observations of the same surface based on 14 cameras (with around 9 degrees difference of view angle between each camera).

### 5.2.4. Atmospheric Correction

For the atmospheric correction, data of aerosol load and type from weather forecast atmospheric models are used. Uncertainties of these aerosol data can result in important biases of the albedo estimates (>10% [1]). In comparison, water vapour and other gaseous or molecular uncertainty have a lower impact (admittedly less than 2% of all biases cumulated). Furthermore, *Proud* et al. [54] explore the effectiveness of simplified radiative transfer models for atmospheric correction such as SMAC. They show that it does not offer a high level of accuracy under many sets of atmospheric conditions (relative errors are between a few percent and 10% for up to 20% of encountered observations).

### 5.2.5. Review Process

Finally, another limitation is also today resulting from a non-homogenised procedure for the evaluation of the product maturity. For the case of surface albedo, the Land Product Validation group of the joint Committee on Earth Observation Satellites (LPV/CEOS) issued recommendations on how to address both direct and indirect validations [55], based on a number of metrics that can be analytically calculated. Metrics recommended by LPV/CEOS for the evaluation are usually used. However, differences exist between EUMETSAT and EU review processes.

EUMETSAT follows a similar revision process as before the publication of an article in a peer review journal. EUMETSAT involves independent reviewers to assess the validation. They participate in the revision process of validation reports with a volunteer approach. Validations are done by the developers or independent people. However, EUMETSAT encourages validation done by independent people. After this revision process, EUMETSAT decides to disseminate (or not) the products.

In contrast, C3S have subcontractors being partners in the programme who are responsible for the validation and the review of the products before their dissemination. The validation report is called quality assessment report and it is systematically done by non-developers. After this revision process, EU decides to disseminate (or not) the products.

To summarise, EUMETSAT involves independent reviewers to assess the validation report, in contrast to C3S, who have subcontractors being partners in the programme. Validation is always conducted by non-developers in C3S (and not systematically in LSA SAF). This may lead to grey areas in the performances of the different products, and the causes of the discrepancies become sometimes difficult to identify.

## 6. Roadmap for Product Continuity

Table 1 indicates the list of CDR, ICDR and NRT products, the product ID, the principal
characteristics, the temporal coverage, and the access to the data. The mission which would
cover the product continuity is also reported in the last column.

In C3S, we are now preparing the version 2 called C3S-V2 which should be available
by the end of 2021. In this evolution of C3S-V2, snow/ice information over NOAA/AVHRR
period will be now used for the calculation of the albedo over snow/ice areas. Moreover,
a BRDF climatology based on 10 years of VGT-2 will be exploited to better constrain the
inversion and improve the temporal continuity from one sensor to the other. On a practical
aspect, in Equations (7) and (8), $k_{\text{clim}}$ and $Ck_{\text{clim}}$ resulting from this climatology will be
added. This will improve the continuity robustness between the various generations
of sensors. Moreover, a noticeable improvement of the atmospheric correction will be
implemented considering the anisotropy of the surface for taking into account the multi-
scattering effects. Finally, the long-term continuity will be prepared through the use of
Sentinel-3 (S3) observations (after the end-of-life of PROBA-V).

In LSA SAF, the last satellite of the Meteosat series (MSG-11) with on board SEVIRI
has an availability lifetime until 2033. Its continuity will be ensured with the launch in
2021 of MTG-I with, on board, the Flexible Combined Imager (FCI) with at least 20 years of
additional operational services. The last satellite of the series MTG (MTG I4) is planned to
be launched in 2032. The calculation of surface albedo from geostationary meteorological
satellites operated by EUMETSAT will be ensured along these series to, at least, 2041.
Concerning the European polar meteorological satellites, EPS series are composed by
three satellites: Metop-A, Metop-B, and Metop-C with, on board, AVHRR sensor. Metop-
A was launched in 2006 and its lifetime was extended from 2019 to 2022 (from 13 to
16 years). The last one, Metop-C, was launched in 2018 and we can also reasonably
expect a comparable lifetime of, at least, 13 years to 2031. Continuity will be ensured
with the launch in 2021 of Metop-SG with on board METImage (also called Visible and
Infrared Imager, VII) and Multi-viewing Multi-channel Multi-polarisation Imaging (3MI)
instruments. Three satellites will constitute the new Metop-SG A series with the first one
launched in 2023 and the last of the series (Metop-SG A3) in 2037. The calculation of
surface albedo from polar European meteorological satellites will then be ensured along
these series to, at least, 2050. Of course, the scientific algorithm will encounter evolutions
along these the next decades and regular reprocessing (yielding new CDR version) will
be required. EUMETSAT, which has funded research and development activities since
1999 has successfully managed through this LSA SAF programme to set up the scientific
baseline for the retrieval of surface albedo properties. Figure 11 shows the chronogram
of the satellite missions that we have (and will) considered (consider) for the surface
albedo retrieval.

<figure><img src="image_16.png" /><figcaption>Figure 11. Flowchart of the satellite missions used for the retrieval of the surface albedo products using the Météo France algorithm presented in this paper.</figcaption></figure>

7. Access to the Code Sources and Data Policy

The code for the generation of these albedo products is made available in open source
under the CeCILL-C license. We give it the name of 'PYALUS' (Python code for ALbedo
retrieval Using Satellite data). Users shall therefore follow the principles of this license and,
in particular, the rules for the exploitation of the code. In addition, users are kindly requested to duly acknowledge authors of the code: "The PYALUS code (Carrer et al., 2021) was made available to the community by CNRM/Météo-France thanks to the support of EUMETSAT."

The code is accessible at the following location: https://github.com/dominiquecarrer/pyalus.

## 8. Conclusions

Despite several important specificities among the list of products which have been developed by Météo France within the framework of the programmes discussed here, the proposed method for the retrieval of the surface albedo respects, overall, a common strategy with successive steps: (i) atmospheric correction, (ii) spectral harmonisation (facultative), and (iii) BRDF inversion with, in most of the cases, the use of a recursive procedure allowing uncertainty estimations and propagations. At the time being, this method has been applied to 11 sensors: NOAA-7-9-11-14-16-17/AVHRR2-3, SPOT/VGT1-2, Metop/AVHRR-3, PROBA-V, and MSG/SEVIRI. Up to three different albedo products with different characteristics exist today to meet the users' expectations (i.e., LSA-101, LSA-102, and LSA-150). Each sensor specific albedo estimate has been assessed by independent reviewers leading to its operational release. Dissemination is operated through different channels (satellite, ftp, websites) to the community with a timeliness of 3 h today for NRT products. For the climate community needs, CDR products with non-near real-time constraint were developed, and have, for an advantage, the use of atmospheric reanalysis fields, which are more accurate than forecasted ones, as input of the retrieval method. This article presents the method for retrieval of these surface albedos, describes their characteristics, and discusses the differences between the products. Furthermore, the plan of continuity with the next European satellite missions is introduced. For example, Metop/AVHRR-3 albedo will soon become the medium resolution sensor product with the longest NRT data record (supported by the EUMETSAT operational system of satellite data dissemination), since MODIS is approaching the end of its life-cycle. Metop-SG/METimage will ensure its continuity thanks to consistent production of data sets guaranteed till 2050 by the European member states.

The surface albedo products provide valuable information for many communities: climate modelling, environmental management, natural hazards management, and climate change detection. In the framework of LSA SAF and C3S projects, around 40 years of data characterizing the surface are now available. These albedos are derived from various sensors from the 1980s. These programmes also provide a great opportunity to monitor and identify human-induced climate change as a consistent production of data sets and is guaranteed until at least 2050 (arising from future European mission calendars). Moreover, the common strategy which we proposed through the different programmes may offer an unprecedented opportunity for the study of temporal trends affecting surface properties.

**Author Contributions:** Conceptualization, D.C.; software, D.C., and F.P.; formal analysis, D.C., F.P., G.L., I.F.T., I.B., F.C., X.C., S.M., L.S., and J.S.-Z.; project administration, I.F.T. and I.B.; funding acquisition, J.M.-S. and L.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors wish to acknowledge the members of the LSA SAF and C3S/312 teams, as well as Météo France, EUMETSAT and the European Commission for the funding of these activities.

**Institutional Review Board Statement:** No applicable.

**Informed Consent Statement:** No applicable.

**Data Availability Statement:** No applicable.

**Acknowledgments:** We thank C. Vincent and M. Pardé for their assistance in the implementation of the improvement in the algorithms and for their help to visualise the data. We thank E. Swinnen, R. Lacaze and D. Ramon for their collaboration in C3S project. We thank A. Verger and C. Schaaf for their guidance on spectral and BRDF harmonisation procedures.
**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A. Spectral Harmonisation

To calculate the linear coefficients used for the spectral harmonisation, we used the method developed by *Samain* et al. [56]. A three-step approach is considered:

Step 1: A set of *n* synthetic surface reflectance values are generated.

For each simulation, a radiative transfer function is generated using SAIL [44] to model the soil component and using PROSPECT [57] to model the vegetation type component (Figure A1). Then, this transfer function is applied to the solar spectrum to compute the surface response spectrum.

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Simulated surface reflectance values for various soil types. Only 100 simulated surface reflectances are shown but, in practice, 50,000 were generated.</figcaption></figure>

To define the simulation parameters, the parameter distributions must be as realistic as possible in order to generate surface reflectances as close to real-world reflectances as possible. Indeed, the quality of the spectral harmonisation is highly dependent on the quality of the simulation parameters.

The number of simulated surface reflectances has been set to *n* = 50,000 because the linear model coefficients (α) needed such a high number of simulations before stabilising. These surface reflectances have been simulated using various soil types and canopy types to account for their high variability—108 canopy types and 10 soil types have been considered. The leaf type has been randomly set to spherical (50%), planophile (25%) or erectophile (25%). The three key parameters required to simulate the surface reflectances are presented in Table A1, and are:

* The Fraction of Vegetation Cover (FCover) corresponds to the fraction of ground covered by green vegetation. It is a dimensionless parameter that takes values between 0 and 1.

* The Leaf Area Index (LAI) is another dimensionless quantity (m²/m²) that characterises plant canopies. It corresponds to the one-sided green leaf area per unit ground surface area.

* The hotspot parameter characterises the hotspot effect and was added to SAIL by [58]. It corresponds to the ratio between the radius of a single leaf and the canopy height. It is another main variable for the SAIL model [59].
<table><caption>**Table A1.** Key parameters used to simulate the surface reflectance.</caption><thead><tr><th>Name</th><th>Characteristics</th></tr></thead><tbody><tr><td>Fcover</td><td>Uniform distribution on (0,1)</td></tr><tr><td>LAI</td><td>Normal distribution with mean = 4 × Fcover and standard deviation = 2, bounded within (0,10)</td></tr><tr><td>Hotspot parameter</td><td>Randomized on a uniform distribution between 0 and (1 + LAI/8)/2</td></tr></tbody></table>

Step 2: Generating synthetic measurements.

The surface reflectances generated in the previous simulations are used to create synthetic measurements. In fact, the sensor acts as a filter on the signal reflected from the surface. The resulting signal can be expressed in the spectral (wavelength) domain as:

$$R_{\text{sensor,synt}} = \int_{\lambda} R_{\text{surf}}(\lambda) \cdot T_{\text{sensor}}(\lambda) \cdot d\lambda \quad (\text{A1})$$

where $R_{\text{sensor,synt}}$ is the synthetic measurement for the sensor, $R_{\text{surf}}(\lambda)$ is the surface reflectance for wavelength $\lambda$ and $T_{\text{sensor}}(\lambda)$ is the sensor spectral response for wavelength $\lambda$. Each simulation provides a synthetic measurement for each band of each sensor. These are the data used to perform the next step.

Step 3: Calculation of the linear model harmonisation coefficients.

Having collected $n$ values of $R_{\beta}^{(\text{ref sensor})}$ for each band $\beta$ of the reference sensor (here VGT-2) and $n$ values of $R_{\beta'}^{(\text{sensor})}$ for each band $\beta'$ of each other sensor to harmonise, linear regressions were applied. We chose to apply linear regressions using the “elastic net” [60] penalty, which is a standard linear regression with L1 norm and L2 norm penalty functions.

To understand how the elastic net regression works, let us first recall that in standard linear regressions, L2 norm is used for the error between the model and the data. Then, the objective function to optimise is the Residual Sum of Square (RSS):

$$\text{criterion} = \text{RSS} = \sum_{i=1}^{n} (\hat{y}_{i} - y_{i})^{2} \quad (\text{A2})$$

where $y_i$ is the $i$-th data point and $\hat{y}_i$ is the model estimate for the $i$-th data point.

In comparison, the linear regression with the “elastic net” penalty adds a penalty term that depends on the model parameters:

$$\text{criterion} = \frac{1}{2n} \sum_{i=1}^{n} (\hat{y}_i - y_i)^2 + w_0 \left( \frac{(1 - w_1)}{2} \cdot \|\alpha\|_2 + w_1 \cdot \|\alpha\|_1 \right) \quad (\text{A3})$$

where,

$$\|\alpha\|_1 = |\alpha_{\rightarrow \beta}^{(\text{sensor})}| + \sum_{\beta'} |\alpha_{\beta' \rightarrow \beta}^{(\text{sensor})}| \quad \text{is the L1 norm of the model parameters,}$$

$$\|\alpha\|_2 = (\alpha_{\rightarrow \beta}^{(\text{sensor})})^2 + \sum_{\beta'} (\alpha_{\beta' \rightarrow \beta}^{(\text{sensor})})^2 \quad \text{is the L2 norm of the model parameters,}$$

$w_0 = 0.0005$ and $w_1 = 0.01$ are the regularisation parameters and have been carefully chosen to alter the model without degrading its performance (in terms of R$^2$).

This technique results in the coefficients presented in Table A2. One uses the standard deviation parameter to assess the quality of the prediction.
<table><caption>**Table A2.** Coefficients ($\alpha_{\beta \to \beta}^{(AVHRRX)}$ in the first column and $\alpha_{\beta' \to \beta}^{(AVHRRX)}$) of each linear model (Equation (1)) to predict the reflectance values $R_{\beta}^{(VGT2)}$ for the VGT-2 B0, B2, B3 and MIR bands from the reflectance values $R_{\beta'}^{(AVHRRX)}$ for the NOAA-X/AVHRR bands. The standard deviation resulting from the regression is given in the last column. Only NOAA-7, 9, 11, 14, 16, 17 are considered. These coefficients are only valid for snow-free pixels. For pixels flagged with snow/ice, a different set of coefficients is used following the same approach based on snow ice spectra only (for more information see ATBD on https://cds.climate.copernicus.eu/).</caption><thead><tr><th>β'(NOAA-7/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR7)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR7)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR7)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR7)</sup><br>β'=MIR</th><th>Sigma</th></tr></thead><tbody><tr><td>B0</td><td>-0.0523</td><td>0.5813</td><td>0.0960</td><td></td><td>0.0309</td></tr><tr><td>B2</td><td>0.0097</td><td>1.0150</td><td>-0.0203</td><td></td><td>0.0146</td></tr><tr><td>B3</td><td>0.0053</td><td>-0.0746</td><td>1.0608</td><td></td><td>0.0162</td></tr><tr><td>MIR</td><td>0.0216</td><td>0.5062</td><td>0.0873</td><td></td><td>0.0639</td></tr><tr><th>β'(NOAA-9/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR9)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR9)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR9)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR9)</sup><br>β'=MIR</th><th>Sigma</th></tr><tr><td>B0</td><td>-0.0528</td><td>0.5820</td><td>0.0938</td><td></td><td>0.0311</td></tr><tr><td>B2</td><td>0.0088</td><td>1.0192</td><td>-0.0258</td><td></td><td>0.0147</td></tr><tr><td>B3</td><td>0.0066</td><td>-0.0802</td><td>1.0640</td><td></td><td>0.0162</td></tr><tr><td>MIR</td><td>0.2016</td><td>0.5101</td><td>0.0829</td><td></td><td>0.0638</td></tr><tr><th>β'(NOAA-11/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR11)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR11)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR11)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR11)</sup><br>β'=MIR</th><th>Sigma</th></tr><tr><td>B0</td><td>-0.0530</td><td>0.5807</td><td>0.0941</td><td></td><td>0.0312</td></tr><tr><td>B2</td><td>0.0084</td><td>1.0181</td><td>-0.0256</td><td></td><td>0.0146</td></tr><tr><td>B3</td><td>0.0071</td><td>-0.0802</td><td>1.0636</td><td></td><td>0.0162</td></tr><tr><td>MIR</td><td>0.2015</td><td>0.5103</td><td>0.0825</td><td></td><td>0.0638</td></tr><tr><th>β'(NOAA-14/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR14)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR14)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR14)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR14)</sup><br>β'=MIR</th><th>Sigma</th></tr><tr><td>B0</td><td>-0.0533</td><td>0.5706</td><td>0.1037</td><td></td><td>0.0312</td></tr><tr><td>B2</td><td>0.0082</td><td>0.9953</td><td>-0.0036</td><td></td><td>0.0141</td></tr><tr><td>B3</td><td>0.0068</td><td>-0.0382</td><td>1.0236</td><td></td><td>0.0155</td></tr><tr><td>MIR</td><td>0.2020</td><td>0.5036</td><td>0.0880</td><td></td><td>0.0638</td></tr><tr><th>β'(NOAA-16/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR16)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR16)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR16)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR16)</sup><br>β'=MIR</th><th>Sigma</th></tr><tr><td>B0</td><td>-0.0080</td><td>0.6869</td><td>0.1190</td><td>-0.2241</td><td>0.0274</td></tr><tr><td>B2</td><td>-0.0010</td><td>0.9766</td><td>-0.0068</td><td>0.0441</td><td>0.0135</td></tr><tr><td>B3</td><td>0.0072</td><td>-0.0287</td><td>1.0306</td><td>-0.0160</td><td>0.0155</td></tr><tr><td>MIR</td><td>0.0119</td><td>0.0413</td><td>0.0210</td><td>0.9326</td><td>0.0141</td></tr><tr><th>β'(NOAA-17/AVHRR)<br>β(VGT2)</th><th>α<sub>→β</sub><sup>(AVHRR17)</sup></th><th>α<sub>β'→β</sub><sup>(AVHRR17)</sup><br>β'=RED</th><th>α<sub>β'→β</sub><sup>(AVHRR17)</sup><br>β'=NIR</th><th>α<sub>β'→β</sub><sup>(AVHRR17)</sup><br>β'=MIR</th><th>Sigma</th></tr><tr><td>B0</td><td>-0.0346</td><td>0.6920</td><td>0.1961</td><td>-0.2742</td><td>0.0304</td></tr><tr><td>B2</td><td>-0.0032</td><td>0.9317</td><td>-0.0476</td><td>0.1430</td><td>0.0132</td></tr><tr><td>B3</td><td>-0.0015</td><td>-0.0749</td><td>1.0085</td><td>0.0793</td><td>0.0155</td></tr><tr><td>MIR</td><td>0.0710</td><td>-0.2946</td><td>-0.5074</td><td>1.8174</td><td>0.0346</td></tr></tbody></table>
## References

1. Carrer, D.; Meurey, C.; Ceamanos, X.; Roujean, J.L.; Calvet, J.C.; Liu, S. Dynamic mapping of snow-free vegetation and bare soil albedos at global 1km scale from 10-year analysis of MODIS satellite products. *Remote Sens. Environ.* **2014**, *140*, 420–432.
[CrossRef]

2. Davin, E.L.; de Noblet-Ducoudré, N. Climatic impact of global-scale deforestation: Radiative versus nonradiative processes. *J. Clim.* **2010**, *23*, 97–112. [CrossRef]

3. Pielke, R.A.; Avissar, R. Influence of landscape structure on local and regional climate. *Landsc. Ecol.* **1990**, *4*, 133–155. [CrossRef]

4. Ross, J. Role of phytometric investigations in the studies of plant stand architecture and radiation regime. In *The Radiation Regime and Architecture of Plant Stands*; Dr W. Junk Publishers: The Hague, The Netherlands, 1981; pp. 9–11.

5. Cedilnik, J.; Carrer, D.; Mahfouf, J.F.; Roujean, J.L. Impact assessment of daily satellite-derived surface albedo in a limited-area NWP model. *J. Appl. Meteorol. Climatol.* **2012**, *51*, 1835–1854. [CrossRef]

6. Carrer, D.; Lafont, S.; Roujean, J.-L.; Calvet, J.-C.; Meurey, C.; Le Moigne, P.; Trigo, I.F. Incoming Solar and Infrared Radiation Derived from METEOSAT: Impact on the Modeled Land Water and Energy Budget over France. *J. Hydrometeorol.* **2012**, *13*, 504–520. [CrossRef]

7. Schaepman-Strub, G.; Schaepman, M.E.; Painter, T.H.; Dangel, S.; Martonchik, J.V. Reflectance quantities in optical remote sensing—Definitions and case studies. *Remote Sens. Environ.* **2006**, *103*, 27–42. [CrossRef]

8. Bojinski, S.; Verstraete, M.; Peterson, T.C.; Richter, C.; Simmons, A.; Zemp, M. The Concept of Essential Climate Variables in Support of Climate Research, Applications, and Policy. *Bull. Am. Meteorol. Soc.* **2014**, *95*, 1431–1443. [CrossRef]

9. Leroy, M.; Deuzé, J.L.; Bréon, F.M.; Hautecoeur, O.; Herman, M.; Buriez, J.C.; Tanré, D.; Bouffiès, S.; Chazette, P.; Roujean, J.-L. Retrieval of atmospheric properties and surface bidirectional reflectances over the land from POLDER/ADEOS. *J. Geophys. Res.* **1997**, *102*, 17023–17037. [CrossRef]

10. Justice, C.O.; Vermote, E.; Townshend, J.R.; Defries, R.; Roy, D.P.; Hall, D.K.; Salomonson, V.V.; Privette, J.L.; Riggs, G.; Strahler, A.; et al. The moderate resolution imaging spectroradiometer (MODIS): Land remote sensing for global change research. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 1228–1249. [CrossRef]

11. Wanner, W.; Strahler, A.H.; Hu, B.; Lewis, P.; Muller, J.P.; Li, X.; Schaaf, C.B.; Barnsley, M.J. Global retrieval of BRDF and albedo over land from EOS MODIS and MISR data: Theory and algorithm. *J. Geophys. Res.* **1997**, *102*, 17143–17161. [CrossRef]

12. Strahler, A.H.; Muller, J.P. MODIS Science Team Members. MODIS BRDF/Albedo Product: Algorithm Theoretical Basis Document, Version 5.0. 1999. Available online: https://modis.gsfc.nasa.gov/data/atbd/atbd_mod09.pdf (accessed on 21 January 2021).

13. Wang, Z.; Schaaf, C.B.; Chopping, M.J.; Strahler, A.H.; Wang, J.; Roman, M.O.; Rocha, A.V.; Woodcock, C.E.; Shuai, Y. Evaluation of Moderate-resolution Imaging Spectroradiometer (MODIS) snow albedo product (MCD43A) over tundra. *Remote Sens. Environ.* **2012**, *117*, 264–280. [CrossRef]

14. Schaaf, C.B.; Gao, F.; Strahler, A.H.; Lucht, W.; Li, X.; Tsang, T.; Strugnell, N.C.; Zhang, X.; Jin, Y.; Muller, J.P.; et al. First Operational BRDF, Albedo and Nadir Reflectance Products from MODIS. *Remote Sens. Environ.* **2002**, *83*, 135–148. [CrossRef]

15. Wang, Z.; Schaaf, C.B.; Sun, Q.; Shuai, Y.; Román, M.O. Capturing rapid land surface dynamics with Collection V006 MODIS BRDF/NBAR/Albedo (MCD43) products. *Remote Sens. Environ.* **2018**, *207*, 50–64. [CrossRef]

16. Diner, D.J.; Beckert, J.C.; Reilly, T.H.; Bruegge, C.J.; Conel, J.E.; Kahn, R.A.; Martonchik, J.V.; Ackerman, T.P.; Davies, R.; Gerstl, S.A.; et al. Multi-angle imaging spectro-radiometer (MISR) in- strument description and experiment overview. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 1072–1087. [CrossRef]

17. Pinty, B.; Roveda, F.; Verstraete, M.M.; Gobron, N.; Govaerts, Y.; Martonchik, J.V.; Diner, D.J.; Kahn, R.A. Surface albedo retrieval from Meteosat. *J. Geophys. Res.* **2000**, *105*, 18099–18134. [CrossRef]

18. Pinty, B.; Lattanzio, A.; Martonchik, J.V.; Verstraete, M.M.; Gobron, N.; Taberner, M.; Widlowski, J.L.; Dickinson, R.E.; Govaerts, Y. Coupling diffuse sky radiation and surface Albedo. *J. Atmos. Sci.* **2005**, *62*, 2580–2591. [CrossRef]

19. Muller, J.P.; López, G.; Watson, G.; Shane, N.; Kennedy, T.; Yuen, P.; Lewis, P.; Fischer, J.; Guanter, L.; Domench, C.; et al. The ESA GlobAlbedo project for mapping the Earth’s land surface albedo for 15 years from European sensors. In Proceedings of the IEEE Geoscience and Remote Sensing Symposium (IGARSS), Munich, Germany, 22–27 July 2012.

20. Muller, J.P. GlobAlbedo Final Validation Report; University College London: London, UK, 2013. Available online: http://www.globalbedo.org/docs/GlobAlbedo_FVR_V1_2_web.pdf (accessed on 10 June 2020).

21. Liang, S.; Zhao, X.; Liu, S.; Yuan, W.; Cheng, X.; Xiao, Z.; Zhang, X.; Liu, Q.; Cheng, J.; Tang, H.; et al. A long-term Global LAnd Surface Satellite (GLASS) data-set for environmental studies. *Int. J. Digit. Earth* **2013**, *6*, 5–33. [CrossRef]

22. Govaerts, Y.M.; Wagner, S.; Lattanzio, A.; Watts, P. Joint retrieval of surface reflectance and aerosol optical depth from MSG/SEVIRI observations with an optimal estima- tion approach: 1. Theory. *J. Geophys. Res.* **2010**, *115*, D02203. [CrossRef]

23. Carrer, D.; Roujean, J.L.; Hautecoeur, O.; Elias, T. Daily estimates of aerosol optical thickness over land surface based on a directional and temporal analysis of SEVIRI MSG visible observations. *J. Geophys. Res.* **2010**, *115*, D10. [CrossRef]

24. Barnsley, M.J.; Strahler, A.H.; Morris, K.P.; Muller, J.P. Sampling the surface bidirectional reflectance distribution function (BRDF): Evaluation of current and future satellite sensors. *Remote Sens. Rev.* **1994**, *8*, 271–311. [CrossRef]

25. Strahler, A.H. Vegetation canopy reflectance modeling—Recent developments and remote sensing perspectives. *Remote Sens. Rev.* **1997**, *15*, 179–194. [CrossRef]

26. Hu, B.; Lucht, W.; Li, X.; Strahler, A.H. Validation of kernel-driven models for global modeling of bidirectional reflectance. *Remote Sens. Environ.* **1997**, *62*, 201–214. [CrossRef]
27. Wanner, W.; Li, X.; Strahler, A.H. On the derivation of kernels for kernel-driven models of bidirectional reflectance. *J. Geophys. Res. Atmos.* **1995**, *100*, 21077–21089. [CrossRef]

28. Lucht, W.; Schaaf, C.B.; Strahler, A.H. An algorithm for the retrieval of albedo from space using semiempirical BRDF models. *IEEE Trans. Geosci. Remote Sens.* **2000**, *38*, 977–998. [CrossRef]

29. Geiger, B.; Carrer, D.; Franchisteguy, L.; Roujean, J.L.; Meurey, C. Land Surface Albedo derived on a daily basis from Meteosat Second Generation Observations. *IEEE Trans. Geosci. Remote Sens.* **2008**, *46*, 3841–3856. [CrossRef]

30. Carrer, D.; Moparthy, S.; Lellouch, G.; Ceamanos, X.; Pinault, F.; Freitas, S.C.; Trigo, I.F. 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.* **2018**, *10*, 1262. [CrossRef]

31. Wang, D.D.; Liang, S.L.; He, T.; Yu, Y.Y. Direct Estimation of Land Surface Albedo from VIIRS Data: Algorithm Improvement and Preliminary Validation. *J. Geophys. Res. Atmos.* **2013**, *118*, 12577–12586. [CrossRef]

32. Wang, D.D.; Liang, S.L.; Zhou, Y.; He, T.; Yu, Y.Y. A New Method for Retrieving Daily Land Surface Albedo from VIIRS Data. *IEEE Trans. Geosci. Remote Sens.* **2017**, *55*, 1765–1775. [CrossRef]

33. Trigo, I.F.; Dacamara, C.C.; Viterbo, P.; Roujean, J.L.; Olesen, F.; Barroso, C.; Camacho-de-Coca, F.; Carrer, D.; Freitas, S.C.; García-Haro, J.; et al. The satellite application facility for land surface analysis. *Int. J. Remote Sens.* **2011**, *32*, 2725–2744. [CrossRef]

34. Lellouch, G.; Carrer, D.; Vincent, C.; Pardé, M.C.; Frietas, S.; Trigo, I.F. Evaluation of Two Global Land Surface Albedo Datasets Distributed by the Copernicus Climate Change Service and the EUMETSAT LSA SAF. *Remote Sens.* **2020**, *12*, 1888. [CrossRef]

35. Sanchez-Zapero, J.; Camacho, F.; Martinez-Sanchez, E.; Lacaze, R.; Carrer, D.; Pinault, F.; Benhadj, I.; Munoz-Sabater, J. Quality Assessment of PROBA-V Surface Albedo V1 for the Continuity of the Copernicus Climate Change Service. *Remote Sens.* **2020**, *12*, 2596. [CrossRef]

36. Rahman, H.; Dedieu, G. SMAC: A simplified method for the atmospheric correction of satellite measurements in the solar spectrum. *Intern. J. Remote Sens.* **1994**, *15*, 123–143. [CrossRef]

37. Vermote, E.F.; Tanré, D.; Deuze, J.L.; Herman, M.; Morcette, J.J. Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: An Overview. *IEEE Trans. Geosci. Remote Sens.* **1997**, *35*, 675–686. [CrossRef]

38. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. *Numerical Recipes in Fortran*; Cambridge University Press: Cambridge, UK, 1992.

39. Li, X.; Gao, F.; Wang, J.; Strahler, A. A priori knowledge accumulation and its application to linear BRDF model inversion. *J. Geophys. Res.* **2001**, *106*, 11925–11935. [CrossRef]

40. Hagolle, O.; Lobo, A.; Maisongrande, P.; Cabot, F.; Duchemin, B.; de Pereyra, A. Quality assessment and improvement of temporally composited products of remotely sensed imagery by combination of VEGETATION 1 & 2 images. *Remote Sens. Environ.* **2004**, *94*, 172–186.

41. Carrer, D.; Roujean, J.L.; Meurey, C. Comparing operational MSG/SEVIRI land surface albedo products from Land SAF with ground measurements and MODIS. *IEEE Trans. Geosci. Remote Sens.* **2010**, *48*, 1714–1728. [CrossRef]

42. Pokrovsky, I.O.; Pokrovsky, O.M.; Roujean, J.-L. Development of an operational procedure to estimate surface albedo from the SEVIRI/MSG observing system in using POLDER BRDF measurements. *Remote Sens. Environ.* **2003**, *87*, 198–242. [CrossRef]

43. Hook, S.J. ASTER Spectral Library. 1998. Available online: http://speclib.jpl.nasa.gov (accessed on 10 June 2020).

44. Verhoef, W. Light scattering by leaf layers with application to canopy reflectance modeling, the SAIL model. *Remote Sens. Environ.* **1984**, *16*, 125–141. [CrossRef]

45. Toté, C.; Swinnen, E.; Sterckx, S.; Clarijs, D.; Quang, C.; Maes, R. Evaluation of the SPOT/VEGETATION Collection 3 reprocessed dataset: Surface reflectances and NDVI. *Remote Sens. Environ.* **2017**, *201*, 219–233. [CrossRef]

46. Sterckx, S.; Benhadj, I.; Duhoux, G.; Livens, S.; Dierckx, W.; Goor, E.; Adriaensen, S.; Heyns, W.; Van Hoof, K.; Strackx, G.; et al. The PROBA-V mission: Image processing and calibration. *Int. J. Remote Sens.* **2014**, *35*, 2565–2588. [CrossRef]

47. Schmetz, J.; Pili, P.; Tjemkes, S.; Just, D.; Kerkmann, J.; Rota, S.; Ratier, A. An introduction to Meteosat Second Generation (MSG). *Bull. Am. Meteor. Soc.* **2002**, *83*, 977–992. [CrossRef]

48. Carrer, D.; Smets, B.; Ceamanos, X.; Roujean, W.H.J.-L.; Lacaze, R. *Copernicus Global Land SPOT/VEGETATION and PROBA-V Surface Albedo Products*—1 Km Version 1; Algorithm Theoretical Basis Document, Issue 2.11, Copernicus Global Land Operations Vegetation and Energy CGLOPS-1. *Framework Service Contract N°* 199494; Join Research Center (JRC): Ispra, Italy, 2018.

49. Geiger, B.; Samain, O. *Albedo Determination, Algorithm Theoretical Basis Document, of the CYCLOPES Project*; Version 2.0.; Météo-France/CNRM: Paris, France, 2004; p. 25.

50. Grenfell, T.C.; Perovich, D.K. Seasonal and spatial evolution of albedo in a snow ice land ocean environment. *J. Geophys. Res.* **2004**, *109*, C01001. [CrossRef]

51. Planque, C.; Carrer, D.; Roujean, J.-L. Analysis of MODIS albedo changes over steady woody covers in France during the period of 2001–2013. *Remote Sens. Environ.* **2017**, *191*, 13–29. [CrossRef]

52. Lebourgeois, F.; Pierrat, J.-C.; Perez, V.; Piedallu, C.; Cecchini, S.; Ulrich, E. Simulating phenological shifts in French temperate forests under two climatic change scenarios and four driving global circulation models. *Int. J. Biometeorol.* **2010**, *54*, 563–581. [CrossRef]

53. Meirink, J.F.; Roebeling, R.A.; Stammes, P. Inter-calibration of polar imager solar channels using SEVIRI. *Atmos. Meas. Tech.* **2013**, *6*, 2495–2508. [CrossRef]
54. Proud, S.R.; Fensholt, R.; Rasmussen, M.O.; Sandholt, I. A compar- ison of the effectiveness of 6S and SMAC in correcting for atmospheric interference of Meteosat Second Generation images. *J. Geophys. Res.* **2010**, *115*, D17209. [CrossRef]

55. Wang, Z.; Schaaf, C.; Lattanzio, A.; Carrer, D.; Grant, I.; Román, M.; Camacho, F.; Yu, Y.; Sánchez-Zapero, J.; Nickeson, J. Best Practice for Satellite Derived Land Product Validation. In *Global Surface Albedo Product Validation Best Practices Protocol*; Version, 1.0; Wang, Z., Nickeson, J., Román, M., Eds.; Land Product Validation Subgroup (WGCV/CEOS): Washington, DC, USA, 2019; p. 45. [CrossRef]

56. Samain, O.; Geiger, B.; Roujean, J.L. Spectral Normalization and Fusion of Optical Sensors for the Retrieval of BRDF and Albedo: Application to VEGETATION, MODIS, and MERIS Data Sets. *IEEE Trans. Geosci. Remote Sens.* **2006**, *44*, 3166–3179. [CrossRef]

57. Jacquemoud, S.; Baret, F. PROSPECT: A model of leaf optical properties spectra. *Remote Sens. Environ.* **1990**, *34*, 75–91. [CrossRef]

58. Kuusk, A. Determination of vegetation canopy parameters from optical measurements. *Remote Sens. Environ.* **1991**, *37*, 207–218. [CrossRef]

59. Jacquemoud, S.; Verhoef, W.; Baret, F.; Bacour, C.; Zarco-Tejada, P.J.; Asner, G.P.; François, C.; Ustin, S.L. PROSPECT+SAIL models: A review of use for vegetation characterization. *Remote Sens. Environ.* **2009**, *113*, 56–66. [CrossRef]

60. Zou, H.; Hastie, T. Regularisation and Variable Selection via the Elastic Net. *J. R. Stat. Soc. Ser. B* **2004**, *67*, 301–320. [CrossRef]