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

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.

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,10,11,12,13,14,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,25,26,27,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 }$ $\left(j=1,\cdots ,n\right)$ 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 }_{\mathrm{vj}}$ and solar zenith angles ${\theta }_{\mathrm{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).

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 }^{\left(sensor,harm\right)}$ that would have been observed by a VGT-2 sensor on band β, by combining the reflectances $R_{\beta {}^{\prime }}^{\left(sensor\right)}$ measured by each sensor at their own spectral bands $\beta {}^{\prime }$. Equation (1) exhibits the model applied here for the harmonisation step of each sensor.

$$R_{\beta }^{\left(sensor,harm\right)}={\alpha }_{\to \beta }^{\left(sensor\right)}+{\displaystyle \sum }_{\beta {}^{\prime }}{\alpha }_{\beta {}^{\prime }\to \beta }^{\left(sensor\right)}·R_{\beta {}^{\prime }}^{\left(sensor\right)}+\epsilon$$

(1)

where ${\alpha }_{\to \beta }^{\left(sensor\right)}$ is the constant coefficient (depending on $\beta$) of the linear model for the sensor to harmonise into the VGT-2 band $\beta$, ${\alpha }_{\beta {}^{\prime }\to \beta }^{\left(sensor\right)}$ is the weight (depending on $\beta$and $\beta {}^{\prime }$) of the band $\beta {}^{\prime }$ in the linear model for the sensor to harmonise into the VGT-2 band $\beta$. In addition, ε represents the Gaussian noise of the model but it is not a coefficient as such. For the cases where a TOC reflectance $R_{\beta {}^{\prime }}^{\left(sensor\right)}$ is missing for one or several bands $\beta {}^{\prime }$, the corresponding harmonised reflectances $R_{\beta }^{\left(sensor,harm\right)}$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 }_{\to \beta }^{\left(sensor\right)}$and ${\alpha }_{ \beta {}^{\prime }\to \beta }^{\left(sensor\right)}$ 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, φ)

(2)

with $k_{\beta }={\left(k_{0\beta }, k_{1\beta }, k_{2\beta }, \dots \right)}^T$ and $f={\left(f_0, f_1, f_2, \dots \right)}^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,\beta }$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. 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.

Product	Sensor (Coverage)	Product ID (Type)	Production Frequency	Composite Window	Temporal Characteristic Time Scale	Spatial Scale	Composition Method	BRDF Model	Atmospheric Correction	Documentation and Data Link Accesses	Temporal Coverage	Product Continuity
Albedo (CDR)	NOAA-X/AVHRR (global)	C3S-V1	10-day	20-day	10-day	1–4 km	recursive	Li-Sparse Reciprocal	Hygeos	C3S-V1 access	1981–2015	C3S-V2
Albedo (CDR)	SPOT/VGT (global)	C3S-V1	10-day	20-day	10-day	1 km	recursive	Li-Sparse Reciprocal	Hygeos	C3S-V1 access	1998–2014	C3S-V2
Albedo (ICDR)	PROBA-V (global)	C3S-V0	10-day	30-day	20-day	1 km	recursive	Roujean	VITO	C3S-V0 access	2014+	C3S-V2
Albedo (NRT)	SPOT/VGT (global)	VGP-P	10-day	30-day	20-day	1 km	recursive	Roujean	VITO	CGLS access	1998–2014
Albedo (NRT)	PROBA-V (global)	PROBA-V L2A	10-day	30-day	20-day	1 km	recursive	Roujean	VITO	CGLS access	2014–2020
Albedo (NRT)	MSG/SEVIRI (Africa, Europe, South America)	LSA-101	1-day	1-day	5-day	SEVIRI grid	recursive	Roujean	Météo France	LSA-101 access	2005+	LSA-107 (MTG/FCI)
Albedo (NRT)	MSG/SEVIRI (Africa, Europe, South America)	LSA-102	10-day	30-day	30-day	SEVIRI grid	classic	Roujean	Météo France	LSA-102 access	2009+	LSA-108 (MTG/FCI)
Albedo (CDR)	MSG/SEVIRI (Africa, Europe, South America)	LSA-150	10-day	30-day	30-day	SEVIRI grid	classic	Roujean	Météo France	LSA-150 access	2005–2015
Albedo (NRT)	Metop/AVHRR (global)	LSA-103	10-day	20-day	10-day	1 km	recursive	Li-Sparse Reciprocal	Météo France	LSA-103 access	2015+	LSA-110, 111 (Metop-SG/METimage, 3MI)

. 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.

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 σ[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_jw_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^TA) k = A^T b

(3)

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

C_k = (A^TA)⁻¹

(4)

The diagonal elements C_ii of this matrix represent the variance σ²[k_i] of the respective parameters k_i. If the matrix A^TA 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 ± σ_ap [k_i]

(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₀, k₁ and k₂. 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^{\prime }=\left|\begin{array}{c}A\\ \begin{array}{ccc}{\sigma }_{ap}^{-1}\left[k_0\right]& 0& 0\\ 0& {\sigma }_{ap}^{-1}\left[k_1\right]& 0\\ 0& 0& {\sigma }_{ap}^{-1}\left[k_2\right]\end{array}\end{array}\right|$$

(6)

and extending the vector b to b′ = (b₁, …, b_n, k_0apσ_ap⁻¹[k₀], k_1apσ_ap⁻¹[k₁], k_2apσ_ap⁻¹[k₂]^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⁻¹)k = A^T b + C_ap⁻¹ k_ap

(7)

and

C_k = (A^T A + C_ap⁻¹)⁻¹

(8)

with k_ap = (k_0ap, …, k_m−1ap)^T and C_ap the covariance matrix. For uncorrelated a priori information on the parameters, C_ap = diag(σ_ap²[k₀], …, σ_ap²[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 Ck_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{array}{c}{k}_{ap}=k_{in}\\ C_{ap}=C_k^{in}·\delta ,\end{array}$$

(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.

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:

k₀ = k_0reg ± σ_reg[k₀]
k₁ = k_1reg ± σ_reg[k₁]
k₂ = k_2reg ± σ_reg[k₂]

(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₁ = 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 β, the directional-hemispherical ‘dh’ and bi-hemispherical ‘bh’ albedos are estimated from the integrals of the kernel functions $f={\left(f_0, f_1, f_2\right)}^T$ which are noted respectively $I^{dh}\left({\theta }_s\right)={\left(I_0^{dh}\left({\theta }_s\right), I_1^{dh}\left({\theta }_s\right),I_2^{dh}\left({\theta }_s\right)\right)}^T$ and $I^{bh}={\left(I_0^{bh}, I_1^{bh},I_2^{bh}\right)}^T$. This angular integration is described in [29,30] by inserting $R_{\beta }$ (from Equation (2)) in the following integral equation:

$$\begin{gathered} a_{\beta }^{dh}\left({\theta }_s\right)=\frac{1}{\pi }{{\displaystyle \int }}_{2\pi }^{}{R}_{\beta }\left({\theta }_s, {\theta }_v, \phi \right)\cos {\theta }_{v}d{\mathrm{\Omega }}_v \\ {a}_{\beta }^{dh}\left({\theta }_s\right)=k_{\beta }^T{I}^{dh}{\left({\theta }_s\right)}^{} \mathrm{and} a_{\beta }^{bh}=k_{\beta }^T{I}^{bh} \\ \sigma \lceil a_{\beta }\rceil ={\sqrt{I^T{C}_{k}I}}_{} \end{gathered}$$

(11)

The resulting narrowband albedo values serve as approximations for the spectral albedo at the central band with wavelength λ_β.

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 =\left[{\lambda }_1,{\lambda }_2\right]$ are derived from the spectral quantities by applying a linear transformation of the form

$$a_{\gamma }=c_{0\gamma }+{\displaystyle \sum }_{\beta }{c}_{\beta \gamma }{a}_{\beta }$$

(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^{\mathrm{dh}}\left({\theta }_{\mathrm{ref}}\right)$, 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 (σ²[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).

Based on these input data and the algorithm detailed in Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 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.

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 Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 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.

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 EUMETSAT. 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 Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 2.6, Metop/AVHRR-3 ETAL albedos were generated. Figure 8 shows an example of retrieval of 1 July 2015.

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 Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5 and Section 2.6, Meteosat Daily Albedos (called MDAL) were generated. Figure 9 shows an example of retrieval of the 1 July 2015.

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 Section 3.2.1 and Section 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.

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 $\delta$ 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 $\delta$ 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. VGT-2 spectral bands, full width at half maximum and centre wavelength in parentheses.

  Band	VGT-2 (µm)
  Blue (B0)	0.439–0.476 (0.458)
  Red (B2)	0.616–0.690 (0.653)
  NIR (B3)	0.783–0.892 (0.838)
  SWIR (MIR)	1.584–1.685 (1.635)

  );
• 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.

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. MSG/SEVIRI and Metop/AVHRR-3 spectral bands, full width at half maximum and centre wavelength in parentheses.

Band	SEVIRI (µm)	AVHRR-3 (µm)
Blue (B0)	-	0.439–0.476 (0.458)
Red (B2)	0.56–0.71 (0.635)	0.616–0.690 (0.653)
NIR (B3)	0.74–0.88 (0.81)	0.783–0.892 (0.838)
SWIR (MIR)	1.50–1.78 (1.64)	1.584–1.685 (1.635)

). 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.

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_clim and Ck_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.

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.