Quality Assessment of PROBA-V Surface Albedo V1 for the Continuity of the Copernicus Climate Change Service

The Copernicus Climate Change Service (C3S) includes estimates of Essential Climate Variables (ECVs) as a series of Climate Data Records (CDRs) derived from satellite data. The C3S Surface Albedo (SA) v1.0 CDR is composed of observations from National Oceanic and Atmospheric Administration (NOAA) Very High Resolution Radiometers (AVHRR) (1981–2005), and VEGETATION sensors onboard Satellites for the Observation of the Earth (SPOT/VGT) (1998–2014) and Project for Onboard Autonomy satellite (PROBA-V) (2014–2020), and will continue with Sentinel-3 (from 2020 onwards). The goal of this study is to assess the uncertainties associated with the C3S PROBA-V SA v1.0 product, with a focus on the transition from SPOT/VGT to PROBA-V. The methodology followed the good practices recommended by the Land Product Validation sub-group (LPV) of the Working Group on Calibration and Validation (WGCV) of the Committee on Earth Observing Satellites (CEOS) for the validation of satellite-derived global albedo products. Several performance criteria were evaluated, including an intercomparison with National Aeronautics and Space Agency (NASA) MCD43A3 C6 products. C3S PROBA-V SA v1.0 and MCD43A3 C6 showed similar completeness but had higher fractions of missing data than C3S SPOT/VGT SA v1.0. C3S PROBA-V SA v1.0 showed similar precision (~1%) to MCD43A3 C6, improving the results of SPOT/VGT SA v1.0 (2–3%), but C3S PROBA-V SA v1.0 provided residual noise in the near-infrared (NIR). Good spatio-temporal continuity between C3S PROBA-V and SPOT/VGT SA v1.0 products was found with a mean bias between ±2%. The comparison with MCD43A3 C6 showed positive mean biases (5%, 8%, and 12% for visible, NIR and total shortwave, respectively). The accuracy assessment with ground measurements showed a median error of 18.4% with systematic overestimation (positive bias of 11.5%). The percentage of PROBA-V retrievals complying with the C3S target requirements was 28.6%.


1.
To evaluate the product completeness and spatio-temporal consistency of C3S PROBA-V SA v1.0 products over global conditions compared to SPOT/VGT SA v1.0 (and MCD43A3 C6 for benchmarking) to verify the continuity of the C3S time series; 2.
To assess the accuracy and the associated uncertainties of the products via direct validation with ground measurements; and 4.
To assess the compliance of the product with regards to climate user requirements.
The remainder of this paper is structured as follows: Sections 2 and 3 describe the datasets used in the study, which are for satellite-based products and ground data, respectively; Section 4 presents the validation methodology; Sections 5 and 6 gather and discuss the intercomparison and validation results, respectively; and the conclusions are summarized in Section 7.

Remote Sensing Surface Albedo Products
The main features of the three SA products investigated in this work (C3S PROBA-V SA v1.0, C3S SPOT/VGT SA v1.0 and MCD43A3 C6) are summarized in Table 1. The information of the spectral characteristics of the four optical bands of the VEGETATION sensors onboard the PROBA and SPOT satellites, and its equivalent Moderate resolution Imaging Spectroradiometer (MODIS) bands, is provided in Table 2. The quality flag information of each product was used to filter low-quality pixels (Table 3). For the C3S PROBA-V SA v1.0 products, the land pixels showing an input status "out of range" or "invalid" or "saturated" in the B2 and B0 channels were discarded. The SPOT/VGT SA v1.0 pixels where the algorithm failed were not considered in the validation exercise. Additionally, two ancillary variables were also taken into account: The uncertainty (ERR) and the mean age (AGE, in number of days) of the observations that are used to produce the SA. The SPOT/VGT SA v1.0 pixels with associated uncertainty greater than 0.2 and an AGE greater than 20 were discarded. For MODIS C6, the quality flags are given in the MCD43A2 product, and those pixels with fewer than seven valid observations were discarded. Table 3. Quality flag (QFLAG) information used to filter pixels flagged as 'low-quality'.

Product
Quality PROBA-V has operated since May 2013 at an altitude of 820 km in a sun-synchronous orbit with a local overpass time at launch of 10:45 a.m. However, because the satellite has no onboard propellant, this local overpass time gradually differs from the at-launch value [3]. The payload consists of three identical cameras, which are equipped with a very compact Three Mirror Anastigmat telescope. Each camera has two focal planes, one for the short-wave infrared (SWIR) and one for the visible and near-infrared (VNIR) bands, with a slightly different off-nadir along track viewing direction. PROBA-V observes four spectral bands: Blue (called B0: centered at 0.463 µm), Red (B1: 0.655 µm), NIR (B2: 0.845 µm), and SWIR (MIR: 1.600 µm). The target on the ground is imaged at different times and with different viewing angles. This specific technical concept makes the angular configurations of the observations in the VNIR and SWIR bands different, and hence, different angular information is provided. Observations are taken at resolutions from 100 to 180 m at nadir and up to 350 m and 660 m at the swath extremes for the VNIR and SWIR channels, respectively [38]. The final PROBA-V data [39], disseminated by the PROBA-V Ground Segment at 100 m, 300 m and 1 km resolutions, have been available on [40] since the end of 2013. The C3S PROBA-V SA v1.0 products (example in Figure 1) are freely available in the C3S CDS [11]. They also provide Directional-Hemispherical albedos (AL-DH)-also called Black-Sky Albedo (BSA)-and the Bi-Hemispherical albedos (AL-BH)-also called White-Sky Albedo (WSA)-in three broadband spectral domains (visible [VI: 0.4-0.7 µm], NIR [NI: 0.7-4 µm] and total shortwave [BB: 0.3-4 µm]). These broadband albedos are taken from the CGLS [10]. In addition, the spectral AL-DH and AL-BH albedos in the four spectral channels (see Table 2) are made available in the C3S. These products have been produced every 10 days from the end of 2013 onwards, with the production dates on the 3rd, 13th, and 23rd of each month and delivered with a 12 day lag in near real-time. The spatial resolution of the grid is 1/112 • , resulting in a pixel size of approximately 1 km at the equator. Apart from the layers corresponding to albedo retrievals, the ancillary layers corresponding to their respective errors (ERRs), the associated quality flag (QFLAG) and the number of valid observations during the synthesis period (NMOD) are provided. The information of each layer is described in the Product User Guide and Specification (PUGS) document [41].
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 32 The methodology, which is described in the Algorithm Theoretical Basis Document (ATBD) [35], follows the approach of separating the atmospheric correction, directional reflectance normalization, and albedo determination. First, the top-of-atmosphere (TOA) data are processed to obtain the cloud-free top-of-canopy (TOC) reflectances. Then, the spectral TOC reflectances acquired under different solar-viewing configurations during the synthesis period are normalized by inverting a linear kernel-driven model [42]. The synthesis period is 30 day and a semi-Gaussian weighting function with the maximum weight on the last observation of the period was selected for near real-time production. Then, the spectral albedos are computed using the angular integration of the kernel functions with the retrieved parameters for each pixel. Finally, the broadband albedo is defined as a linear combination of the spectral albedo values in the available spectral channels. The narrow to broadband conversion coefficients are applied both for the black-sky and for the white-sky albedos. 2.2. C3S SA v1.0 Based on SPOT/VGT Data C3S SPOT/VGT SA v1.0, which is available from C3S CDS [11], provides a CDR of the global black-sky and white-sky albedo estimates for the period from 1999 to 2014 at a ground sampling distance of 1/112° and a temporal frequency of 10 days using a compositing window of 20 days. The dates of production are the 10 th , 20 th and 30 th of each month. As for PROBA-V, the albedo quantities are provided for the three broadband domains (visible [0.4-0.7 µm], NIR [0.7-4 µm] and total shortwave [0.3-4 µm]) and in the four spectral channels (see Table 2) of the VEGETATION sensors. The albedo products are based on the data acquired by the VEGETATION-1 and VEGETATION-2 sensors aboard the SPOT-4 and SPOT-5 satellites, respectively.
The retrieval algorithm [6,7,9] was developed based on the previous developments of the SA products [6][7][8][9] based on MSG/SEVIRI and Metop/AVHRR in the framework of the EUMETSAT LSA SAF project [4], and it was later adapted here to these other sensors in CGLS and C3S [43]. The input data were Collection 3 of SPOT/VGT, which corrects the bug previously detected regarding the incorrect calculation of the Sun-Earth distance of Collection 2, includes improved cloud and snow/ice detection, and revises the absolute radiometric calibration for the entire archive [44]. The main differences with the C3S PROBA-V SA v1.0 algorithm are related to the use of a different Bidirectional Reflectance Distribution Function (BRDF) model and temporal composite scheme. PROBA-V SA v1.0 implemented the Ross_Thick kernel for volumetric scattering and the Roujean kernel for geometric scattering [42], whereas SPOT/VGT SA v1.0 implemented the Ross_Thick kernel [42] for volumetric scattering and the Li_Sparse_Reciprocal kernel for geometrical scattering [45]. The synthesis period of PROBA-V SA v1.0 is 30 days using a semi-Gaussian weighting function. The SPOT/VGT SA v1.0 retrieval methodology uses a 20-day temporal composite, meaning that each surface albedo product is built from the valid SPOT/VGT observations corresponding to the 20-day The methodology, which is described in the Algorithm Theoretical Basis Document (ATBD) [35], follows the approach of separating the atmospheric correction, directional reflectance normalization, and albedo determination. First, the top-of-atmosphere (TOA) data are processed to obtain the cloud-free top-of-canopy (TOC) reflectances. Then, the spectral TOC reflectances acquired under different solar-viewing configurations during the synthesis period are normalized by inverting a linear kernel-driven model [42]. The synthesis period is 30 day and a semi-Gaussian weighting function with the maximum weight on the last observation of the period was selected for near real-time production. Then, the spectral albedos are computed using the angular integration of the kernel functions with the retrieved parameters for each pixel. Finally, the broadband albedo is defined as a linear combination of the spectral albedo values in the available spectral channels. The narrow to broadband conversion coefficients are applied both for the black-sky and for the white-sky albedos.
2.2. C3S SA v1.0 Based on SPOT/VGT Data C3S SPOT/VGT SA v1.0, which is available from C3S CDS [11], provides a CDR of the global black-sky and white-sky albedo estimates for the period from 1999 to 2014 at a ground sampling distance of 1/112 • and a temporal frequency of 10 days using a compositing window of 20 days. The dates of production are the 10th, 20th and 30th of each month. As for PROBA-V, the albedo quantities are provided for the three broadband domains (visible [0.4-0.7 µm], NIR [0.7-4 µm] and total shortwave [0.3-4 µm]) and in the four spectral channels (see Table 2) of the VEGETATION sensors. The albedo products are based on the data acquired by the VEGETATION-1 and VEGETATION-2 sensors aboard the SPOT-4 and SPOT-5 satellites, respectively.
The retrieval algorithm [6,7,9] was developed based on the previous developments of the SA products [6][7][8][9] based on MSG/SEVIRI and Metop/AVHRR in the framework of the EUMETSAT LSA SAF project [4], and it was later adapted here to these other sensors in CGLS and C3S [43]. The input data were Collection 3 of SPOT/VGT, which corrects the bug previously detected regarding the incorrect calculation of the Sun-Earth distance of Collection 2, includes improved cloud and snow/ice detection, and revises the absolute radiometric calibration for the entire archive [44]. The main differences with the C3S PROBA-V SA v1.0 algorithm are related to the use of a different Bidirectional Reflectance Distribution Function (BRDF) model and temporal composite scheme. PROBA-V SA v1.0 implemented the Ross_Thick kernel for volumetric scattering and the Roujean kernel for geometric scattering [42], whereas SPOT/VGT SA v1.0 implemented the Ross_Thick kernel [42] for volumetric scattering and the Li_Sparse_Reciprocal kernel for geometrical scattering [45]. The synthesis period of PROBA-V SA v1.0 is 30 days using a semi-Gaussian weighting function. The SPOT/VGT SA v1.0 retrieval methodology uses a 20-day temporal composite, meaning that each surface albedo product is built from the valid SPOT/VGT observations corresponding to the 20-day period preceding the calculation date. At the BRDF inversion step, the previous inversion result is used as a priori information. Hence, a recursive temporal composition of the information over a longer time period (approximately three times longer than the production period) can be achieved to guarantee the coherence and spatial completeness of the product. The "age" (AGE parameter in the output product) of the clear-sky observations exploited in the recursive inversion scheme is an important piece of information for potential applications, and is therefore, also made available to users. This age corresponds to the mean age with respect to the date of the calculation of the clear observations considered for the albedo calculation. More details are given in the ATBD [36].
Sánchez-Zapero et al. [46] showed that the C3S SPOT/VGT SA v1.0 products were good quality based on most of the criteria evaluated, reaching validated stage 1 in the CEOS LPV hierarchy. However, two main limitations were found: (i) Some temporal noise existed at the short-time scale; and (ii) the product was not able to capture some snow events, showing large uncertainties (>0.2) in those cases. Comparisons with ground measurements (461 samples, 2000-2005 period) from 15 FLUXNET stations showed similar overall uncertainty (RMSD = 0.05) as other satellite products (MCD43A3 C6, GlobAlbedo, and GLASS), but a positive bias (14%) was found.

NASA MCD43 C6 Based on MODIS Data
The MODIS BRDF/Albedo (MCD43A3) Collection 6 (C6) dataset, which is available from [47], provides both directional albedo at local solar noon and bi-hemispherical albedo for MODIS bands 1-7 and for three broadbands (visible [0.3-0.7 µm], NIR [0.7-5.0 µm], and total shortwave [0.3-5.0 µm]). The MCD43A3 albedo quantities are delivered at a resolution of 500 m in a sinusoidal projection. They have been produced every day since 2000 with a synthesis period of 16 days. Data from both Terra and Aqua satellites are used in the generation of this product.
The MODIS albedo algorithm uses atmospherically corrected reflectance data (MOD09 product) to establish the best fit to a linear kernel-driven BRDF model, with the exception of the observations flagged as "cloud", "cirrus high" or "aerosol high". Like for C3S SPOT/VGT SA v1.0, the parametric BRDF model uses the Ross_Thick kernel for volumetric scattering and the Li_Sparse_Reciprocal kernel for geometrical scattering [45,48]. A full retrieval of the model is attempted if there are seven or more high-quality observations well distributed over the viewing hemisphere during the 16 day synthesis period. When the number of observations is strictly less than seven and strictly greater than two, or if observations are not well sampled or do not well fit the BRDF model, a back-up algorithm with prior information is used. Then, a fill value is stored if the number of observations is strictly less than three, and the separated snow-free gap-filled products are also accessible [49]. Then, the BRDF model parameters are used for estimating spectral albedos from angular integration. The broadband albedos are computed using the spectral to broadband conversion approach [50]. The MCD43 C6 products use an improved back-up database, which is pixel-based updated from the latest full inversion as opposed to the land cover-based database used in the previous Collection 5.
MCD43A3 SA products have reached CEOS LPV validation stage 3 [51]. Sánchez-Zapero et al. [46] reported an overall uncertainty (RMSD) of 0.053 and a low negative bias (−11.9%) of MCD43A3 C6 compared with the ground data from 15 FLUXNET homogeneous sites for the 2000-2005 period (653 samples), and a slight overall uncertainty (RMSD = 0.032) was reported compared with European FLUXNET measurements over snow-free conditions [52]. Existing studies of the previous collection 5 indicate that the accuracy of the MODIS shortwave broadband albedo meets the requirements (<5%) for both snow-free and snow-covered surfaces [31][32][33][34].

Ground-Based Observations for Validation (GBOV) Database
The CGLS GBOV [53] initiative aims at facilitating the use of observations from operational ground-based monitoring networks and their comparison to Earth Observation (EO) products. In the case of the SA, GBOV provides the Reference Measurements (RMs) and upscaled Land Products (LPs) based on the RM data [54,55] from over 24 sites coming from different networks, such as BSRN, FLUXNET, or SURFRAD. Currently, the GBOV database is made available for the 2012-2018 period [53]. The RM GBOV data corresponds to blue-sky albedo, which is defined as the fraction between the downward and the upward shortwave radiative flux. The diffuse fraction is also included in the RM dataset.
The innovative approach for the generation of GBOV LPs could be very useful for the validation of satellite SA products since it can be applied to both homogeneous and heterogeneous land surfaces. However, we can expect large discrepancies over heterogeneous sites compared to homogeneous sites [54,56]. The CEOS LPV albedo validation protocol [21] recommends the use of ground values measured at the flux tower (i.e., the RM) for the direct validation of EO products. For those reasons, GBOV RMs are used in this study for the accuracy assessment. Note that RMs are provided daily in the GBOV database with a typical temporal step of 30 or 1 min depending on the station. The estimation of the uncertainty associated with a solar radiation measurement by a commercial pyranometer is around 5% (at a 95% confidence level) for daily values under ideal conditions [57]. Since satellite products are defined to provide SA retrievals at solar local noon, the RMs used in this study have been taken at the same time. For this study, the concomitant period to PROBA-V was used (i.e., 2014-2018). The main characteristics of the 20 GBOV sites providing the RMs at solar local noon for the period under study are summarized in Table 4.

Analysis of the Spatial Representativeness of Tower Measurements
The reference albedo measured from towers covers a circular footprint that varies with the tower height. It is unlikely that the footprint of the ground measurements exactly matches the satellite pixel sizes. Then, the spatial representativeness of the tower-based measurements should be evaluated to minimize the issues associated with spatial representativeness in the point-to-pixel comparison [31,32].
The semivariogram [58,59] is one of the most efficient tools for describing spatial representativeness. Semivariograms were computed for the stations under evaluation using Sentinel-2 surface reflectances at a spatial resolution of 10 m near-nadir. Two different periods (leaf-off and leaf-on) periods during the year were considered to evaluate the spatial representativeness of the region around the ground tower. The spatial attributes (e.g., the range, sill and nugget) were computed by fitting the variogram estimator to an isotropic spherical model [60], and they can reveal the spatial variability of land surfaces and the scaling effects associated with remotely sensed data [31][32][33][34][61][62][63]. The methodology adopted in this study for the estimation of the geostatistical indexes is based on the comparison of the variogram model parameters retrieved at different spatial resolutions (i.e., from 1.0 km 2 to 1.5 km 2 squared subsets). Four different geostatistical attributes were used [16,32]: The relative coefficient of variation (R CV ), the scale requirement index (R SE ), the relative strength of the spatial correlation (R ST ) and the relative proportion of the structural variation (R SV ). The four geostatistical attributes can be combined in a compact metric (ST score [31,32], see Equation (1)), which represents a spatial representativeness score using R SE as the primary marker and the others like secondary weights. When the spherical variogram model does not provide a good fit to the variogram estimator, another indicator (RAW score [31,32], see Equation (2) could be used to provide a spatial representativeness score, which is only based on the R CV .
Both scores are directly proportional to the site representativeness. We consider that sites are heterogeneous or not spatially representative when both scores are lower than 2.0 (see Table 5) because large differences, due to spatial heterogeneity, are expected for scores below this threshold [16]. Based on it, the ground measurements coming from Boulder (USA_BAO), Renon (ITA_REN), Rock Springs (USA_PSU) and Table Mountain (USA_TBL) sites are not used for the accuracy assessment during the leaf-off seasonal period. In addition, the Cabauw (NET_CAB), Goodwin Creek (USA_GCM) and Gobabeb (NAM_GOB) sites were discarded during leaf-on seasonal period. Figure 2 shows two examples of sites of stations spatially representative (top side), and two examples of not spatially representative (bottom side).

Uncertainty Requirements
The accuracy assessment results were analyzed against predefined uncertainty levels based on a review of the existing user requirements from the Global Climate Observing System (GCOS), the World Meteorological Organization (WMO) and the key performance indicators of C3S.
In the last update of the GCOS requirements [64], there is a distinction between the products targeted for "adaptation" and "modelling" applications that results in different needs for the horizontal resolution. Modelling requirements (i.e., uncertainty Max [5%; 0.0025], see Table 6) are

Uncertainty Requirements
The accuracy assessment results were analyzed against predefined uncertainty levels based on a review of the existing user requirements from the Global Climate Observing System (GCOS), the World Meteorological Organization (WMO) and the key performance indicators of C3S.
In the last update of the GCOS requirements [64], there is a distinction between the products targeted for "adaptation" and "modelling" applications that results in different needs for the horizontal resolution. Modelling requirements (i.e., uncertainty Max [5%; 0.0025], see Table 6) are the focus of this study since modelling is the main application targeted. Other requirements come from the WMO [65], which aids in the setting of the priorities to be agreed upon by WMO members and their space agencies for enhancing the space-based GCOS system. The WMO specifies the requirements for the SA for climatologic applications at three quality levels ( Table 6): Threshold (minimum requirement), breakthrough (significant improvement) and goal (optimum, no further improvement required). The stated "goal" WMO uncertainty requirement of 5% is, thus, equivalent to the GCOS requirement in relative terms. Apart from the GCOS and WMO, the Key Performance Indicator (KPI) of the maximum accuracy being between 10% and 0.01 was defined in the C3S program [66]. Based on the existing requirements, three different levels (i.e., optimal, target and threshold) are predefined in this study (Table 7) with the aim to verify whether the results fit the purpose ( Table 6). The optimal level (Max [5%, 0.0025]) was selected according to the GCOS uncertainty requirement and is equivalent to the WMO goal. The target level (Max [10%, 0.01]) is equivalent to the C3S KPI and partly equivalent to the WMO breakthrough level. Finally, the threshold level (Max [20%, 0.02]) was adopted from the WMO. Poor performances of the product correspond to values above the threshold level (minimum requirement). Figure 3 displays the selected uncertainty levels as a function of the product values. Based on the existing requirements, three different levels (i.e., optimal, target and threshold) are predefined in this study (Table 7) with the aim to verify whether the results fit the purpose ( Table 6). The optimal level (Max [5%, 0.0025]) was selected according to the GCOS uncertainty requirement and is equivalent to the WMO goal. The target level (Max [10%, 0.01]) is equivalent to the C3S KPI and partly equivalent to the WMO breakthrough level. Finally, the threshold level (Max [20%, 0.02]) was adopted from the WMO. Poor performances of the product correspond to values above the threshold level (minimum requirement). Figure 3 displays the selected uncertainty levels as a function of the product values.

Validation Methods
The validation methodology follows the CEOS LPV good practice protocol for the validation of satellite-derived albedo products [21], and the validation metrics are presented in Section 4.2.1. The different strategies for product intercomparison and direct validation are described in Sections 4.2.2 and 4.2.3, respectively. Figure 4 shows the global distribution of the sampling used in this study. The product intercomparison is evaluated over a 725-site land validation network (LANDVAL) of sites [67], which is designed to globally represent the variability of land surface types, and was used as the spatial sampling to evaluate these criteria. This network also includes 20 well-known desert calibration sites [26] for the precision evaluation, due to their high temporal stability. For the direct validation, the 20 selected GBOV homogeneous stations with available ground data (see Table 5) were used.

Validation Methods
The validation methodology follows the CEOS LPV good practice protocol for the validation of satellite-derived albedo products [21], and the validation metrics are presented in Section 4.2.1. The different strategies for product intercomparison and direct validation are described in Sections 4.2.2 and 4.2.3, respectively. Figure 4 shows the global distribution of the sampling used in this study. The product intercomparison is evaluated over a 725-site land validation network (LANDVAL) of sites [67], which is designed to globally represent the variability of land surface types, and was used as the spatial sampling to evaluate these criteria. This network also includes 20 well-known desert calibration sites [26] for the precision evaluation, due to their high temporal stability. For the direct validation, the 20 selected GBOV homogeneous stations with available ground data (see Table 5) were used.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 32 validation, the 20 selected GBOV homogeneous stations with available ground data (see Table 5) were used. The definitions of the completeness, precision, uncertainty and accuracy (Table 8), which are applicable to SA validation, are drawn from the experimental recommendations of the Joint Committee for Guides in Metrology (JCGM) regarding the expression of uncertainty in

Validation Metrics
The definitions of the completeness, precision, uncertainty and accuracy (Table 8), which are applicable to SA validation, are drawn from the experimental recommendations of the Joint Committee for Guides in Metrology (JCGM) regarding the expression of uncertainty in measurement [68] and from the GCOS [64]. Table 8. Validation metrics for product validation.

Quantity
Validation Metric

Completeness
Gap size distribution (spatial and temporal) Gap length

Accuracy
Median Error Box-plots of the bias per albedo range

Uncertainty
Root mean square deviation (RMSD) Scatter-plots of match-ups (MAR Linear models and correlation) Completeness is the proportion of valid retrievals over an observation domain at any given time. It is, therefore, mandatory to document the completeness of the product (i.e., the distribution in space and time of missing data).
Two aspects of the precision, which is also called the repeatability, are evaluated: The intra-annual and inter-annual precision. The intra-annual precision (smoothness or δ function) corresponds to the temporal noise assumed to have no serial correlation within a season. The δ function is computed for each triplet of consecutive observations [69,70] as the absolute value of the difference between the center and the corresponding linear interpolation between the two extremes, and the median of the δ values is provided as an indicator of the intra-annual precision of satellite albedo products [21]. The inter-annual precision (i.e., the dispersion of the albedo values from year to year) was assessed over the 20 desert calibration sites between two consecutive years, and the median absolute deviation is provided as an indicator of the inter-annual precision [21].
Accuracy is the degree of "closeness of the agreement between the result of a measurement and a true value of the measurand" [68]. Commonly, accuracy is used to describe systematic errors and measure statistical bias, but the best practice is to provide the median error as an indicator of accuracy [21].
Uncertainty is a "parameter, associated with the result of a measurement, that characterizes the dispersion of the values that could reasonably be attributed to the measurand" [68]. Uncertainty includes the bias and precision errors and can be estimated by the Root Mean Square Deviation (RMSD). Additionally, the linear model fits are used to quantify the goodness of fit. For this purpose, Major Axis Regressions (MARs) were computed instead of Ordinary Least Squares (OLS) because it is specifically formulated to handle errors in both the x and y variables [71]. Other metrics are used, such as the number of samples (N), which is indicative of the power of the validation; or the correlation coefficient (R, estimated as the Pearson coefficient), which indicates the descriptive power of the linear accuracy test. The percentage of pixels within the predefined uncertainty levels (Table 7, Figure 3) is also quantified.

Satellite Product Intercomparison
The intercomparison of satellite products must account for the differences in the spatial and temporal characteristics of the different datasets. Consequently, a spatial support area and temporal support period were defined consistent with the C3S SA characteristics (Table 1). Thus, the MCD43A4 C6 products were re-sampled to Plate Carrée using a 1 km spatial sampling grid. Satellite products were compared using the closest date of their different temporal composites, using the 10 day temporal frequency of C3S SA products.
The completeness and precision of the different products were evaluated and compared, as well as the temporal consistency among the several products, which was qualitatively assessed per biome type. The uncertainty estimates from the intercomparisons provided an overall figure of the spatio-temporal consistency between products. Scatter-plots and the associated metrics between pairs of products were analyzed and complemented with box-plots of the bias per bin albedo value. Our analysis focuses on the consistency during the 6-month overlap period between PROBA-V SA v1.0 and SPOT/VGT SA v1.0, whereas the intercomparison with MCD43A3 C6 was conducted using one year (2014) of data.

Ground-Based Validation
The main steps in the accuracy assessment of albedo products include the generation of blue-sky albedo [30] for a direct comparison with in-situ measurements and the test of the spatial representativeness of the in situ albedometer footprints for the satellite pixel resolution of interest according to in situ measurements standards [31,32]. Consequently, a careful selection of ground points and the characterization of their spatial representativeness are crucial for a meaningful point-to-pixel comparison, as presented in Section 3.2. For the first step, the blue-sky albedo from a satellite is estimated from the retrieved AL-DH-BB and AL-BH-BB satellite EO product under study, and weighted by the fraction of diffuse down-welling shortwave radiation from the ground station [30]. The next step is the careful selection of the homogeneous sites, which are similar to the footprints of the satellite pixel resolution of interest. The accuracy assessment of the SA satellite products was performed against the measurements coming from the 20 selected GBOV homogeneous stations presented in Section 3.
The study here was extended to the period from 2014 to 2018 to have the maximum concomitant samples from the GBOV dataset to PROBA-V observations. The accuracy assessment of MCD43A3 C6 for the same period and sampling was performed for benchmarking. This exercise was carried out at a resolution of 1 km, which is equivalent to one pixel in the case of PROBA-V SA v1.0, and averages of 2 × 2 pixels in the case of MCD43A3 C6. The 10 day frequency of C3S products was used as temporal sampling, and the average values of the daily ground data were computed to compare with satellite estimations during the corresponding temporal composite window of each product (see Table 1). This analysis was performed over the best quality pixels of C3S PROBA-V SA v1.0 and MCD43A3 C6, according to QFLAGs (Table 3).

Product Completeness
The global map of the percentage of gaps (i.e., values with missing data) during one whole year of data (2014) for the PROBA-V SA v1.0 products ( Figure 5) shows poor spatio-temporal continuity over latitudes higher than 45 • north and over the equatorial belt, with a percentage of missing values up to 100% in some pixels over these areas. Note that the information from the quality flags was not considered here in the computation of the gaps (i.e., a gap or missing value corresponds to a fill value in the product). C3S products are not provided over areas out of the region from latitudes of 60 • S to 75 • N.
The temporal evolution of missing values (Figure 6a The global map of the percentage of gaps (i.e., values with missing data) during one whole year of data (2014) for the PROBA-V SA v1.0 products ( Figure 5) shows poor spatio-temporal continuity over latitudes higher than 45º north and over the equatorial belt, with a percentage of missing values up to 100% in some pixels over these areas. Note that the information from the quality flags was not considered here in the computation of the gaps (i.e., a gap or missing value corresponds to a fill value in the product). C3S products are not provided over areas out of the region from latitudes of 60ºS to 75ºN.  The distribution of the temporal length of the missing values ( Figure 6b) allows one to better understand the impact of the gaps for monitoring the temporal variations. PROBA-V SA v1.0 shows that approximately 65% of the gaps are shorter than 30 days and approximately 35% are shorter than 10 days (time resolution of the product). Similarly, typically, approximately 35-40% of the SPOT/VGT SA v1.0 gaps are shorter than 10 days, and approximately 75% of gaps are shorter than 30 days during the overlap period between SPOT/VGT and PROBA-V. MCD43A3 C6 has a shorter length of gaps, with approximately 40% of the gaps corresponding to 5 days.

Temporal Consistency
The temporal profiles of the different SA products (PROBA-V SA v1.0, SPOT/VGT SA v1.0 and MCD43A3 C6) in the three broadband domains (visible, NIR, and total shortwave) were analyzed over the 725 LANDVAL sites for each main biome type (Figure 7). The analytical focus of the transition between SPOT/VGT and PROBA-V and the 2013-2014 period was represented. All the satellite products are displayed at the center of their temporal composite windows (30 days in the case of PROVA-V SA v1.0, 20 days in the case of SPOT/VGT SA v1.0 and 16 days in the case of MCD43A3 C6). Note that the information of the PROBA-V SA v1.0 QFLAG was also displayed in these graphs: Filled dots correspond to pixels flagged as good quality, and unfilled dots correspond to pixels flagged as low-quality (land pixels with bit 6, 10 or 11 to 1) according to Table 3.
A seasonality effect was observed with the sign of the bias between C3S PROBA-V SA v1.0 and SPOT/VGT SA v1.0 during the short overlap period over desert areas for the NIR and total spectrum. PROBA-V SA v1.0 tends to provide slightly higher values than SPOT/VGT SA v1.0 from December 2013 to February 2014, and the opposite trend was found from March 2014 to May 2014. Additionally, the SPOT/VGT SA v1.0 temporal trajectories show some temporal noise over desert sites, which is not observed for the other satellite products (PROBA-V SA v1.0 and MCD43A3 C6).  (Table 3).

Spatio-Temporal Consistency.
The overall spatio-temporal consistency between PROBA-V SA v1.0 and the reference products is assessed over the LANDVAL network sites considering all good quality observations according to In the case of PROBA-V, filled dots correspond to 'good quality' pixels, and unfilled dots correspond to pixels flagged as 'low-quality' according to QFLAG (Table 3).
For evergreen broadleaved forests, some temporal noise was observed in all satellite products. However, PROBA-V SA v1.0 and MCD43A3 C6 seem to provide less noise (i.e., flatter temporal trajectories) than SPOT/VGT SA v1.0, which seems to be more realistic for this biome type. For the rest of the biomes, the PROBA-V SA v1.0 profiles follow the temporal trends of SPOT/VGT SA v1.0 and MCD43A3 C6. The presence of rapid changes due to snow events, the occurrence of stable profiles and the phenological changes were consistent among the three datasets. However, PROBA-V SA v1.0 displays slightly large variability compared to other satellite products in the NIR domain (also affecting the total spectrum). Note that the use of the C3S PROBA-V SA v1.0 Quality Flag in northern latitudes removes valid snow observations in most cases.
A seasonality effect was observed with the sign of the bias between C3S PROBA-V SA v1.0 and SPOT/VGT SA v1.0 during the short overlap period over desert areas for the NIR and total spectrum. PROBA-V SA v1.0 tends to provide slightly higher values than SPOT/VGT SA v1.0 from December 2013 to February 2014, and the opposite trend was found from March 2014 to May 2014. Additionally, the SPOT/VGT SA v1.0 temporal trajectories show some temporal noise over desert sites, which is not observed for the other satellite products (PROBA-V SA v1.0 and MCD43A3 C6).

Spatio-Temporal Consistency
The overall spatio-temporal consistency between PROBA-V SA v1.0 and the reference products is assessed over the LANDVAL network sites considering all good quality observations according to the quality flags (see Table 3).  Table 9), PROBA-V SA v1.0 tends to provide slightly lower values than SPOT/VGT SA v1.0, with small negative mean biases of −2.2% and −2.8% for black-sky and white sky albedos, respectively. Optimal lineal regression relationships from the MAR were found (offset~0 and slope close to 1). Worse results were found in terms of the RMSD (uncertainty), with values of approximately 0.05 (~35%).  For the near-infrared (Figure 8b-e and Table 9), positive biases (PROBA-V SA v1.0 > SPOT/VGT SA v1.0) of 2.3-2.4% were found, with an RMSD of approximately 10% and high correlations (>0.94). PROBA-V SA v1.0 tends to provide higher values than SPOT/VGT SA v1.0 for albedo values lower than 0.5 and the opposite trend for albedo values higher than 0.5 (typically snow cases). In all cases, a median bias close to zero was found, which was within the optimal level of consistency. than 0.5 and the opposite trend for albedo values higher than 0.5 (typically snow cases). In all cases, a median bias close to zero was found, which was within the optimal level of consistency.

AL-DH-VI AL-DH-NI AL-DH-BB AL-BH-VI AL-DH-NI AL-BH-BB
Remarkably low mean biases (<1%) were found for the total shortwave (Figure 8c-f and Table  9), as well as high correlations (>0.93). Total uncertainties (RMSD) of 0.03 (~15%) were found. A systematic positive bias was found for almost all ranges, except for values lower than 0.1 and higher than 0.7, with the median bias typically within the optimal (GCOS) level of consistency (Table 6). Remarkably low mean biases (<1%) were found for the total shortwave (Figure 8c-f and Table 9), as well as high correlations (>0.93). Total uncertainties (RMSD) of 0.03 (~15%) were found. A systematic positive bias was found for almost all ranges, except for values lower than 0.1 and higher than 0.7, with the median bias typically within the optimal (GCOS) level of consistency (Table 6).  Table 10), with overall RMSD of approximately 14-15%. PROBA-V SA v1.0 tends to provide higher values than MCD43A3 C6 for albedo values lower than 0.5, and the opposite trend for albedo values higher than 0.5.
The worse performance in the comparison PROBA-V SA v1.0 versus MCD43A3 C6 was found for the total spectrum (Figure 9c-f and Table 10), with a large positive bias of ~12-13%. As observed for the NIR, PROBA-V tends to provide higher values than MODIS C6 for albedo values lower than 0.5 and the opposite trend for albedo values higher than 0.5 (i.e., snow cases).     The worse performance in the comparison PROBA-V SA v1.0 versus MCD43A3 C6 was found for the total spectrum (Figure 9c-f and Table 10), with a large positive bias of ∼12-13%. As observed for the NIR, PROBA-V tends to provide higher values than MODIS C6 for albedo values lower than 0.5 and the opposite trend for albedo values higher than 0.5 (i.e., snow cases).

•
Compliance with user requirements The compliance matrix of C3S PROBA-V SA v1.0 versus the reference products (SPOT/VGT SA v1.0 and MCD43A3) with predefined uncertainty levels based on different user requirements is presented in Table 11. Table 11. Compliance matrix (percentage of pixels filing the predefined uncertainty levels) of C3S PROBA-V SA v1.0 versus SPOT/VGT SA v1.0 (December 2103 to May 2014 period) and MCD43A3 C6 (2014 year) broadband albedo products. Computation over LANDVAL sites for good quality pixels according to QFLAGs (Table 3).

Intra-Annual Precision
The Probability Density Functions (PDFs) of the intra-annual precision (the so-called smoothness) are analyzed ( Figure 10). The computation was performed over LANDVAL sites during the overlap period between SPOT/VGT, PROBA-V and MODIS products (i.e., December 2013-May 2014 period).


Compliance with user requirements The compliance matrix of C3S PROBA-V SA v1.0 versus the reference products (SPOT/VGT SA v1.0 and MCD43A3) with predefined uncertainty levels based on different user requirements is presented in Table 11. Table 11. Compliance matrix (percentage of pixels filing the predefined uncertainty levels) of C3S PROBA-V SA v1.0 versus SPOT/VGT SA v1.0 (December 2103 to May 2014 period) and MCD43A3 C6 (2014 year) broadband albedo products. Computation over LANDVAL sites for good quality pixels according to QFLAGs (Table 3).

Intra-Annual Precision
The Probability Density Functions (PDFs) of the intra-annual precision (the so-called smoothness) are analyzed ( Figure 10). The computation was performed over LANDVAL sites during the overlap period between SPOT/VGT, PROBA-V and MODIS products (i.e., December 2013-May 2014 period).
The three products present similar distributions of the smoothness (δ). Most of the delta values are below 0.01, which demonstrates the high stability over a short time scale for the albedo products. In all satellite products, worse intra-annual precision (i.e., higher δ values) was found for white-sky albedos compared with equivalent black-sky albedos. The median δ values (Table 12), which are indicative of the intra-annual precision, show improved results (i.e., lower δ values) of PROBA-V SA v1.0 compared to SPOT/VGT SA v1.0 for visible and total shortwave and worse performance for NIR. Both C3S products provide much worse intra-annual precision than MODIS C6. The three products present similar distributions of the smoothness (δ). Most of the delta values are below 0.01, which demonstrates the high stability over a short time scale for the albedo products. In all satellite products, worse intra-annual precision (i.e., higher δ values) was found for white-sky albedos compared with equivalent black-sky albedos. The median δ values (Table 12), which are indicative of the intra-annual precision, show improved results (i.e., lower δ values) of PROBA-V SA v1.0 compared to SPOT/VGT SA v1.0 for visible and total shortwave and worse performance for NIR. Both C3S products provide much worse intra-annual precision than MODIS C6.  . MCD43A3 C6 provides better inter-annual precision, with median absolute deviations lower than 1%. As observed for the intra-annual precision, all products provided worse results for white-sky albedos compared to black-sky albedos (Table 13). Table 13. The median absolute deviation between two consecutive years of the PROBA-V SA v1.0, SPOT/VGT SA v1.0 and MCD43A3 C6 products. Computation over desert calibration sites.

Ground-Based Direct Validation
To investigate the accuracy of C3S PROBA-V SA v1.0 and MCD43A3 C6 satellite albedo products, scatter plots versus field measurements (GBOV RM) were produced for the 2014-2018 period over 20 homogeneous sites (see Section 3) with different vegetation types. Figure 11 and Table 14 show the scatter-plots, and relevant statistics from the direct validation exercise, whilst Table 15 summarizes the compliance of both satellite products with user requirements (predefined in Section 4.1). The relevant statistics per biome type are presented in Tables 16 and 17 Table 14. Relevant statistics of the direct validation of the best quality pixels (Table 3)   Per biome type, PROBA-V SA v1.0 provides a large positive bias for forest (22.4%) and desert (10.4%) sites compared to crop (7.4%) and grassland/shrubland (4.2%) sites. MCD43A3 C6 systematically provides a large negative bias for most biome cases (croplands, grassland/shrublands and desert) except for forests, where a low positive bias was found (2.8%).    (Table 3)  Note that 15.5% of the 1715 PROBA-V SA v1.0 samples achieved the optimal predefined level (i.e., GCOS requirements) and 28.6% of the target level (i.e., C3S KPI), as shown in Table 15. Slightly improved results were found for MCD43A3 C6 (23.7% and 45.5% of the optimal and target levels, respectively).

Discussion
C3S SA v1.0 based on PROBA-V provides continuity to the CDRs of global albedo products in the C3S from December 2013 onwards. Previously, SA products were derived from SPOT/VGT (1998-2014) and NOAA/AVHRR  data. Good spatio-temporal consistency in the transition from SPOT/VGT to PROBA-V for both black-sky and white-sky albedos, with mean biases below ±3% for the overlap period, was found. However, the comparison of C3S PROBA-V SA v1.0 with MCD43A3 C6 SA showed lower spatio-temporal consistency between satellite products with mean biases up to 13%. C3S PROBA-V SA v1.0 (and MCD43A3 C6) displayed more gaps (typically between 10% and 20%) than C3S SPOT/VGT v1.0. Different aspects of the retrieval methodology play important roles in the existing discrepancies between the satellite products under study (C3S PROBA-V SA v1.0, C3 SPOT/VGT SA v1.0 and MCD43A3 C6) starting from the different input data and different atmospheric correction methods; including the BRDF parameterization, temporal compositing and angular integration; and finalizing with the narrow to broadband conversion.
Regarding the input data, each sensor works in different spectral channels (see Table 2). The SPOT/VGT and PROBA-V channels provide very similar spectral characteristics in the Blue, Red and NIR bands, whereas significant differences are found for the SWIR channel. The central wavelengths of the MODIS spectral bands are comparable to the corresponding bands for PROBA-V and VGT, albeit the MODIS spectral bands are narrower. These differences could translate into reflectance discrepancies in regions with high absorption features like in the visible domain, where large uncertainty values were found between pairs of products.
Each satellite processing chain uses its own method for cloud/shadow screening and atmospheric correction according to the spatial, spectral and directional capabilities of each instrument. Cloud or snow contamination is the main reason for missing data in the EO products derived from optical onboard satellite sensors. The conservative PROBA-V cloud detection algorithm [72] is one of the reasons for the higher fraction of missing data of PROBA-V SA v1.0 products compared to SPOT/VGT SA v1.0.
Discrepancies between different albedo estimates can also be attributed to the different BRDF models used [7]. Moreover, the performance of the BRDF model for good clear-sky observations also depends on the number of available looks during the synthesis period and the angular distribution of the sampling. Large BRDF uncertainties are associated with snow targets (as observed for the highest albedo ranges) for which none of these parametric BRDF models were well suited [73]. PROBA-V, SPOT/VGT and MODIS are wide-FOV sensors on polar-orbiting platforms, and a low impact of discrepancies is expected, due to different sun-view configurations. However, the different compositing periods (see Table 1) could play an important role in the differences between satellite products, mainly in rapid SA variations, such as snow events. Different techniques for temporal composite approaches also affect the completeness of EO satellite products. SPOT/VGT SA v1.0 uses a recursive temporal composition approach [9], which is the main reason for the improved completeness compared to the other satellite products (C3S PROBA-V SA v1.0 and MCD43A3 C6) that are computed using classical composite schemes based on predefined temporal windows.
In the last step, the broadband albedos are defined using slightly different spectral regions. The same broadband spectral regions are defined in both C3S PROBA-V SA v1.0 and SPOT/VGT SA v1.0 products, which contribute to a better agreement, whereas the MCD43A3 C6 broadband albedos are slightly differently defined.
In addition, the temporal noise (i.e., large temporal variability) observed for C3S PROBA-V SA v1.0 in the NIR domain through the qualitative inspection of the temporal variations (see Section 5.1.2) has a strong relationship with the precision, since low intra-annual precision was found compared to both reference EO products in this spectral range. In terms of the inter-annual precision, improved results were found for C3S PROBA-V SA v1.0 (∼1%) compared to C3S SA v1.0 based on SPOT/VGT SA data (2%).
Regarding the accuracy with the ground measurements, PROBA-V SA v1.0 provided a similar accuracy (bias of 11.5%) to that found for C3S SPOT/VGT SA v1.0 during the validation exercise [44], where a positive bias (14%) was also reported for a different sampling (i.e., different stations and dates). PROBA-V SA v1.0 also provided a slightly worse accuracy (median error of 18.2%) than MCD43A3 C6 (median error of 11.2%). PROBA-V SA v1.0 tends to overestimate the ground values, whereas MCD43A3 C6 showed the opposite sign of the mean bias. The positive bias of PROBA-V SA v1.0 was mainly observed for forest sites (SA < 0.2) explained in the fact that the Roujean kernel [42] for geometrical scattering may not fit some cover types well, especially dense forest canopies, where it showed a weak hotspot effect [74][75][76]. The negative bias of MCD43A3 C6 is mainly influenced by some outliers detected in Gobabeb (bare soil) and Cabaw (grassland). For the Cabaw case, the lower MCD43A3 C6 values are explained by the persistent cloudiness at the MODIS overpass times [56]. It is important to note that only the satellite retrievals classified as the best quality, according to QFLAGs (Table 3), were used in the direct validation. As observed in the temporal consistency, the use of QFLAGs removes most of the valid snow retrievals in the case of PROBA-V. Then, this exercise is almost equivalent to snow-free conditions, which is more convenient for assessing the uncertainty of satellite EO products, since it is expected that the spatial representativeness of the pixels dropped during the fall and winter months as a consequence of the increased sub-pixel heterogeneity, due to processes, such as non-uniform patterns of snowmelt [77]. The positive bias of C3S PROBA-V SA v1.0 is consistent with previous studies performed on the CGLS, where a positive bias of approximately 22% was found compared to ground measurements over 17 stations [27], and a positive bias of approximately 14% was found compared to National Ecological Observatory Network (NEON) ground data [78]. The accuracy assessment of MCD43A3 C6 was also consistent with that previously reported when comparing the combined TERRA+AQUA albedo product with eight field stations during the spring and summer months of 2003 and 2004 (i.e., equivalent to snow-free conditions), where an accuracy of 0.013 was reported in terms of the absolute bias [79]. It should be noted that the estimation of the uncertainty of ground reference data is around 5% under ideal conditions [58]. Thus, fiducial reference measurements must be characterized with highly calibrated instrumentation at dedicated cal/val sites to better estimate the satellite product uncertainty budget.
The compliance of satellite EO products versus ground data with user requirements showed a low percentage of pixels within GCOS (only <25%) and C3S (<50%) requirements, which indicates the difficulties of achieving these requirements using current products. To further improve the compliance with requirements, it is recommended that EO programs provide the uncertainties associated with the processing chain, mainly related to the sensor calibration and atmospheric correction. In that way, the steps providing higher error can be improved. However, the comparison between satellite C3S SA v1.0 products showed an overall good spatio-temporal consistency in the comparison of SPOT/VGT versus PROBA-V during the 6-month overlapping period (December 2013-May 2014), with typically more than 60% of the samples within the C3S target requirements for the visible domain, and more than 75% of the samples within the requirements for the NIR and total shortwave.

Conclusions
This paper presents the quality assessment results of C3S PROBA-V SA v1.0 products (broadband albedos) through intercomparison with reference satellite products (C3S SPOT/VGT SA v1.0 and MCD43 C6) at the global scale and the direct validation with a representative amount of ground data (1715 samples) across different biomes types (eight stations over forests, six over crops, three over grassland/shrubs, and three over desert). This validation exercise is a novelty in the literature, since very few global and temporally representative SA validation exercises have been published [21], and they are mainly based on MODIS observations [15][16][17]. The validation methodology adopted the guidelines, protocols and metrics defined by the CEOS LPV best practices for the validation of global albedo satellite products [21]. Additional results, such as the comparison of the spectral albedos or the presentation of the satellite products intercomparison per biome type, are not shown in this manuscript for the sake of brevity; however, they can be found in the product quality assessment report [66].
The main conclusions are the following.

•
The good spatio-temporal consistency of C3S PROBA-V and SPOT/VGT SA v1.0 products assures the continuity of the CDRs in terms of the uncertainty between products. However, C3S PROBA-V SA v1.0 (and MCD32A3 C6) provides low numbers of valid retrievals compared to C3S SPOT/VGT SA v1.0. • C3S PROBA-V SA v1.0 shows similar inter-annual precision (~1%) to MCD43A3 C6, improving the results of SPOT/VGT SA v1.0 (2-3%), since they provide some temporal instability over desert calibration targets. Both C3S products provide lower intra-annual precision than MCD43A3 C6, mainly in the NIR domain where some temporal noise was found.

•
The accuracy of the C3S PROBA-V SA v1.0 best quality retrievals with respect to the ground data over a five-year period (2014-2018) showed systematic positive overestimation, which was mainly observed for the lowest albedo ranges (SA < 0.2) over forest sites. Similar uncertainty (RMSD∼ 0.4) was found for MCD43A3 C6 products using the same sampling, showing the opposite sign for the mean bias.

•
Few current satellite EO albedo products comply with the GCOS, C3S and WMO uncertainty requirements.
Additionally, it is important to remark that the use of PROBA-V QFLAGs (bit 6, input status; and bits 10-11, B2-B0 saturation status) removes most of the valid snow retrievals. Therefore, masking out data by means of the QFLAGs is not recommended for snow applications (users should ignore the information of QFLAGs over snow targets for specific applications).
Based on these results, C3S PROBA-V SA v1.0 has reached validation stage 3 in the validation hierarchy of the CEOS LPV [20]. The continuity of the C3S SA CDR time series will be ensured using Sentinel-3 OLCI and SLSTR data, and the algorithm and design of the processing chain are currently being developed. C3S is also developing multi-sensor albedo products combining NOAA/AVHRR, SPOT/VGT and PROBA-V input data. The long-term CDRs, provided by the Copernicus Climate Change Service from 1981 to the present (with the aim of extending to the future), are an added extra compared with the existing EO programs, providing the longest and state-of-the-art albedo products.  Acknowledgments: This study has been undertaken using data from GBOV "Ground Based Observation for Validation" (https://land.copernicus.eu/global/gbov) founded by the European Commission Joint Research Centre FWC932059 as part of the Global Component of the European Union's Copernicus Land Monitoring Service.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: