Assessing the Fitness of Satellite Albedo Products for Monitoring Snow Albedo Trends

Accurate monitoring of albedo trends over snow is essential to evaluate the consequences of the global snow cover retreat on Earth’s energy budget. Satellite observations provide the best way to monitor these trends globally, but their uncertainty increases over snow. Besides, different products sometimes show diverging trends. A better assessment of the fitness of satellite products for monitoring snow albedo trends is needed. We analyze the consistency of black-sky albedo estimates from global long-term products over snow: advanced very-high-resolution radiometer (AVHRR)-based (CLARA-A2.1, GLASS-v4.2), moderate resolution imaging spectroradiometer (MODIS)-based (MCD43C3-v6.1/v6, GLASS-v4.2), multiangle imaging spectro radiometer (MISR)-based (MIL3MLSN-v4), and multisensor (C3S-v1/v2). We use MCD43C3-6.1 as the reference based on a previous comparison against in situ measurements. CLARA-A2.1 is the one most consistent with MCD43C3, but has a low coverage in high latitudes and an artificial albedo decrease since 2015. The study shows the limitations of MIL3MLSN, Global Land Surface Satellite (GLASS), and Copernicus Climate Change Service (C3S) multisensor products over snow. MIL3MLSN has a too-low coverage of albedo over snow. GLASS-AVHRR overestimates albedo in regions with seasonal snow due to delayed snowmelt and underestimates it in permanently snow-covered regions. GLASS-MODIS is more consistent with MCD43C3 at mid-latitudes, and also underestimates albedo in regions with permanent snow and has an increase in missing values after 2011. Both the GLASS datasets are temporally inconsistent with the other products. Despite the improvements from v1 to v2, C3S-v2 has the largest negative bias over snow and discontinuities in the transitions between sensors. The study evidences the difficulties of AVHRR products to provide stable snow albedo estimates in polar regions, particularly before 2000.

vegetation (VGT) cycle, droughts [3], fires [4], and snow and ice changes [5], [6], among other processes. Monitoring the temporal changes in albedo is essential to quantify the impact of these processes on the surface radiation budget. The Intergovernmental Panel on Climate Change (IPCC) differentiates between the direct radiative forcing caused by anthropogenic land cover changes and the feedback triggered by the global retreat of snow and ice cover [7], [8]. The last IPCC Assessment Report (AR6) [9] reported a negative radiative forcing of −0.2 [−0.3 to −0.1] Wm −2 due to land cover changes since 1750 and a positive surface albedo feedback of +0. 35 [0.10-0.60] Wm −2 • C −1 .
Albedo trends can be monitored either with long-term sensors or by combining several sensors. The first approach is used by CLARA-A2.1 (AVHRR) from the Climate Monitoring Satellite Application Facility (CM SAF) [14], [15], Global Land Surface Satellite (GLASS)-AVHRR and GLASS-MODIS from the Beijing Normal University [16], MCD43 products (MODIS) from NASA [17], [18], and the MISR suite from NASA Jet Propulsion Laboratory (JPL) [19]. The second approach is used by the Copernicus Climate Change Service (C3S), which produces a multisensor climate data record from 1981 with observations from NOAA/AVHRR and VGT (SPOT and PROBA-V) sensors [20]. Either combining measurements from the same sensor onboard different satellites (e.g., AVHRR, VGT) or combining observations from different sensors (e.g., C3S multisensor) brings several stability challenges that can lead to the presence of artificial trends or discontinuities in the climate data records. Even when a unique sensor is used, sensor degradation [12], [21] and orbital drifts [15] can also alter the stability of the product. Besides, the different processing algorithms introduce additional sources of uncertainty. The artificial discontinuities and trends introduced by these issues can be larger than the biogeophysical trends analyzed. Assessing the product stability is therefore essential to evaluate its fitness for the specific application being developed. In this sense, GCOS establishes that the temporal stability of albedo satellite products should be better than max(1%, 0.001) per decade [1].
The performance of satellite albedo products is generally assessed by comparing satellite measurements against in situ observations [22], [23], [24], [25], [26]. Most studies analyze short periods typically below five years that are insufficient to assess the product stability. Moreover, most of them are made over snow-free surfaces, whereas the uncertainty of albedo estimations increases over snow due to snow anisotropy and the low satellite viewing angles over snow-covered polar regions [11], [27], [28]. However, several studies are already exploiting satellite products to monitor snow albedo changes and the consequences of the global snow/ice cover retreat, in Greenland [12], [21], [29], [30], [31], Antarctica [32], and globally [13], [33]. Multidecadal satellite observations are also used to monitor the snow albedo feedback by constraining [9] and even replacing [13] global circulation models (GCMs). A better characterization of the uncertainty and temporal stability of satellite albedo products over snow is needed to understand the reliability of these climate monitoring applications.
The goal of this study is to analyze the fitness of long-term satellite products to monitor snow albedo trends globally. For that we intercompare the black-sky albedo retrievals of products covering at least 20 years: AVHRR-based products (GLASS-AVHRR 4.2, CLARA-A2.1) (1981-present), MODIS-based products (MCD43C3 v6 and v6.1, GLASS-MODIS 4.2) (2000-present), the C3S v1 and v2 multisensor surface albedo (1981-present), and the MISR-based MIL3MLSN (Table I). An intercomparison is the only way to assess the spatio-temporal consistency of albedo products globally due to the lack of a dense network of long-term in situ data. However, this intercomparison is based on the results of a previous study [34], in which we evaluated the products against in situ measurements from the Baseline Surface Radiation Network (BSRN) [35] and FLUXNET [36]. Highresolution satellite images were used to select stations with spatially representative measurements of snow albedo. The comparison showed that MCD43C3 v6.1 has the smallest bias over snow (−0.017), the most stable performance with different snow cover conditions, and the most temporally stable performance at the few stations providing long-term measurements. On the contrary, the quality of both GLASS-AVHRR and C3S-v1/v2 albedo decreases over snow, and the bias of both the products varies with the snow cover conditions, underestimating albedo over snow and overestimating snowfree albedo. Based on this, MCD43C3 v6.1 is used as the reference in this study. We extend the previous study spatially, intercomparing the completeness and snow albedo of the products at a global scale, and temporally, analyzing their stability during the MODIS (2001-2019) and the AVHRR era (1982-2019).

II. DATA: SATELLITE PRODUCTS A. GLASS
The GLASS suite is produced by the Beijing Normal University. It provides two types of albedo products using MODIS (500 × 500 m and 0.05 × 0.05 • ) and AVHRR (0.05 × 0.05 • ) observations, respectively. The parameters available are shortwave (300-3000 nm), visible (300-700 nm), and NIR (700-3000 nm) black-and white-sky albedo every eight days. Surface albedo is calculated with the direct-estimation method generating one broadband albedo estimate from each top-of-atmosphere (TOA) reflectance measurement [6]. A POLDER-based bidirectional reflectance distribution function (BRDF) database is coupled with 6S (second simulation of a satellite signal in the solar spectrum) atmospheric radiative transfer simulations to produce the training database that characterizes the relationship between TOA measurements and broadband surface albedo [37]. Then, a linear regression model is fit for each angular bin and land use (VGT, soil, and snow/ice). The fitted models are applied to each satellite observation to obtain the surface broadband albedo estimates. A spatio-temporal filtering algorithm is used to smooth and gap-filled the product [38]. An internal land use classification is used based on the normalized difference VGT index (NDVI) and TOA reflectance in band 490 nm (r 490 ). Both pure snow (r 490 > 0.4) and intermediate class B (0.25 < r 490 < 0.4) use the snow/ice BRDF dataset. The snow mask is not included in the final product. This study evaluates version 4 of GLASS-AVHRR and GLASS-MODIS, both provided at 0.05 × 0.05 • . Version 4 updates the snow/ice BRDF training dataset and includes a new water surface BRDF model for mixed pixels of water/sea ice [39].

B. MCD43C3
MCD43 BRDF/albedo products are produced by NASA from MODIS sensor on-board Terra and Aqua satellites. The maximum spatial resolution offered is 500 m. In this study, we evaluate the MCD43C3 Albedo Product (MODIS Albedo Daily L3 Global 0.05 • CMG) which provides spectral (MODIS bands 1-7) and broadband (visible 300-400 nm, NIR 700-5000 nm, and shortwave 300-5000 nm) whiteand black-sky albedo [17]. MCD43C3 albedo is obtained daily by applying the MCD43 BRDF weighting parameters to a climate modeling grid of 0.05 × 0.05 • . For each daily estimate, 16 days of MODIS observations are used. They are weighted based on their quality, spatial coverage, and the temporal distance from the day of interest. The RossThick-LiSparseReciprocal BRDF [40] is inverted when at least seven observations are available. Otherwise, a lower quality magnitude inversion is made using a database of archetypal BRDF parameters to supplement the observations. The archetypal database is also derived from MODIS observations and contains both snow and snow-free representations. If most of the 16-day window is snow-covered, snow-free observations Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  I   SUMMARY OF THE SATELLITE PRODUCTS EVALUATED IN THE STUDY. ALL PRODUCTS PROVIDE SHORTWAVE ALBEDO  AT LOCAL SOLAR NOON EXCEPT CLARA-A2.1 are discarded, and vice versa. Thus, snow-covered and nonsnow observations are always processed separately. MCD43 products are designed for clear-sky conditions. Therefore, the algorithm may not be specifically robust against conditions of increased haziness and solar zenith angles (SZAs) above 75 • . A quality flag and the percent of snow at each pixel are included in the final products. We evaluate v6 and the new v6.1, which include changes to the response-versus-scan angle approach that affects the reflectance bands for Aqua-Terra/MODIS, corrections to adjust for the optical crosstalk in Terra/MODIS infrared bands, and corrections to the Terra/MODIS forward lookup table update for the period 2012-2017.

C. MIL3MLSN
The MISR suite is produced by NASA JPL exploiting images of the MISR sensor on-board Terra satellite. MISR consists of nine push-broom cameras measuring radiance in four spectral bands at different angles. The suite includes two types of surface albedo products: MISR Level 2 Land surface product (MIL2ASL*) [41] and MISR Level 3 Global Land Surface product (MIL3MLS*) [19]. Level 2 products provide spectral blue-and white-sky surface albedo, as well as hemispherical-directional and bidirectional reflectance factors for each Terra orbit on a Space Oblique Mercator reference grid with 1.1 × 1.1 km. Level 3 products provide shortwave (400-2500 nm) black-sky albedo on a global geographic grid of 0.5 × 0.5 • . In this study, we use the Level 3 product (MIL3MLSN) because our goal is to benchmark end products of shortwave albedo. The retrieval algorithm is based on a modified form of the BRDF model of Rahman et al. [42]. Level 3 shortwave albedo is approximated from Level 2 spectral values using the algorithm described in Weiss et al. [43]. Only retrievals for which smoothed aerosol optical depth (AOD) in the Level 2 Land Surface product is less than 0.3 are used. In Level 3, both Greenland and Antarctica are excluded due to the typically poor quality of MISR aerosol retrievals in these regions. No specific snow treatment is mentioned in the product ATBD [44], and no snow mask is available in the final product.
The processing algorithm was developed by Meteo-France within EUMETSAT LSA-SAF based on MSG/SEVIRI and Metop/AVHRR and later adapted to SPOT and PROBA-V sensors by C3S and Copernicus Global Land Service (CGLS). TOA observations are atmospherically corrected with the Simplified Model for Atmospheric Correction (SMAC) [42] and harmonized to the SPOT/VGT-2 spectral bands. Bands 1 and 2 are used for AVHRR/2 and bands 1, 2, and 3A for AVHRR/3. The RossThick-LiSparse-Reciprocal BRDF is inverted every ten days (10th, 20th, and 30th of each month) using a compositing window of 20 days, considering the previous BRDF inversion with a recursive method. The narrow-to-broadband conversion is made with a linear equation using different coefficients for snow-free and snow-covered pixels [45]. If most observations within the 20-day window are snow-covered, snow-free observations are discarded, and vice versa. The snow and cloud masks are calculated from the bands of each input sensor. No snow mask was used in C3S-v1 AVHRR [20].
In this study, we evaluate the two last versions of the multisensor product (v1 and v2). The main improvements from v1 to v2 [45] are the addition of a snow mask during the AVHRR period, the introduction of a VGT-derived BRDF climatology to reduce data gaps and improve sensor consistency, and a new harmonized pixel identification approach to improve the satellite cross-consistency when dealing with cloud screening and snow detection.

E. CLARA
The Cloud, Albedo, Radiation data record, AVHRR-based, Edition 2.1 (CLARA-A2.1) is produced by the CM SAF using observations from AVHRR-2/3 sensors on-board NOAA and Metop (excluding Metop B) satellites [46]. CLARA-A2.1 provides surface albedo, surface irradiance, and cloud properties globally from 1982 to the present with a spatial resolution of 0.25 × 0.25 • , but the data are processed internally at 5 × 5 km. The albedo parameters provided are the five-day and monthly averages of shortwave black-sky albedo (250-2500 nm). TOA reflectances from AVHRR channels 1 and 2 are atmospherically corrected with SMAC. Surface albedo is calculated separately for snow-free and snow-covered pixels [15]. In snow-free pixels, the anisotropy factor is used to normalize reflectances to a common viewing and illumination geometry (nadir view, SZA = 0). The anisotropy factor is calculated with Wu's equation [47] using semiempirical kernel parameters derived with different NDVI-based coefficients for each land use (barren, cropland, forest, and grassland). Normalized reflectances are integrated hemispherically with Roujean et al. [48] formula and the semiempirical kernel coefficients to produce the black-sky albedo, which corresponds to the SZA at the time of the observation. In contrast, snow albedo is calculated by directly averaging the bidirectional reflectances collected over five days or more. This method assumes that albedo can be estimated as the mean of bidirectional reflectances if the viewing hemisphere is densely sampled during the observation period, which is generally true at polar regions due to the overlap of multiple AVHRR sensors. However, the quality of snow albedo retrievals may decrease at lower latitudes and when the snow conditions change during the window integrated. In all the cases, only observations made at SZA < 70 • are processed. The snow mask is derived with the AVHRR-Polar Platform System (PPS) package [49], which classifies pixels into four categories: clear, cloudy, partly cloudy, and snow-covered. The Ocean and Sea Ice Satellite Application Facility (OSI-SAF) sea ice extent data are used to confirm snow/ice mask in the Arctic and Antarctic. The snow mask is not included in the final product. In this study, we used the monthly averages of black-sky albedo, extending CLARA-A2.1 (January 1982-June 2019) with the Interim Climate Data Record (ICDR) based on CLARA-A2 methods (July 2019-December 2020).

A. ERA5 Snow Cover Mask
An independent snow mask was used to identify snowcovered regions, and within them, to separate historically snow-covered and snow-free months. We used an external mask to homogenize the validation procedure and because not all the products provide their snow masks. The snow mask does not change in time to monitor surface albedo trends caused by snow cover retreat/gain in each region.
For each pixel, the historically snow-covered and snowfree months were defined based on multiyear snow data from ERA5 reanalysis [50]. ERA5 snow time series has several artificial discontinuities and trends due to the addition of new data into the assimilation system [51] but shows good stability and low bias after the assimilation of the Interactive Multisensor Snow and Ice Mapping System (IMS) product [52] in 2004. Therefore, the snow mask was calculated with ERA5 data between 2005 and 2020. Monthly snow cover fraction (SCF) was derived from ERA5 monthly snow depth (SD): min(1, SD[cm]/10). For each month, a pixel was classified as snow-covered if the ERA5 SCF in that pixel was above 0.5 (SCF > 0.5) for at least one year between 2005 and 2020. Both the snow-covered and snow-free months were used in the analysis.

B. MCD12C1 v6: IGBP Land Cover Classification
NASA's MCD12C1 v6 [53] was used to assess the performance of the products over different land cover types. MCD12C1 provides maps of International Geosphere-Biosphere Programme (IGBP), the University of Maryland (UMD), and Leaf Area Index (LAI) classification schemes at 0.05 × 0.05 • from 2001 to 2020. We used the IGBP land cover classification scheme [54], which defines 17 land cover classes. MCD121C1 provides the percent cover of each class per pixel and year. We aggregated the raw files to 0.25 × 0.25 • by summing the area covered by each class. Then, the class with the highest coverage was assigned to each 0.25 × 0.25 • pixel.

C. Preprocessing Satellite Products
The intercomparison is based on the monthly black-sky albedo to include CLARA-A2.1, which only provides blacksky values. Monthly black-sky albedo (α black,m ) was calculated by weighting submonthly black-sky albedo (α black,i ) with the number of days covered by each submonthly estimate (days i ) where days m is the total number of days in the month m, and N m is the number of submonthly albedo values in the month. All the products were aggregated to the coarser spatial grid (MIL3MLSN, 0.5 × 0.5 • ) for completeness analysis. However, the bias and stability were analyzed at 0.25 × 0.25 • after removing MIL3MLSN due to its high number of gaps over snow. In both the cases, products were aggregated by spatially averaging the native resolution into the 0.25 × 0.25 • or 0.5 × 0.5 • grids.
All the products except CLARA provide black-sky albedo at local solar noon. CLARA-A2.1 retrievals represent the solar zenith angle at the time of the satellite observation, discarding images with SZA > 70 • . The difference between averaging albedo observations at any solar elevation angle and averaging solar noon observations was evaluated using ERA5 hourly albedo for one month (March 2000). For each pixel, the monthly mean of all hourly albedo estimations available was compared against the monthly mean of albedo estimations at local solar noon (see Fig. S4). In both the cases, only clear-sky hours with SZA < 70 • were used. The comparison was made with the blue-sky albedo, as it is the parameter provided by ERA5. Fig. S4 shows that monthly averages using all hourly observations are systematically larger than averages based on solar noon observations over the sea due to the specular reflectance of water. However, the difference is negligible over both the snow-free and snow-covered land surfaces. Based on this, we kept CLARA-A2.1 in the intercomparison.

D. Intercomparison Procedure
Completeness was evaluated based on the percentage of valid retrievals in months with SZA < 90 • . Note that throughout this whole article, all the available observations are used independently of their quality flags.
The consistency of snow albedo estimates was analyzed based on the mean bias deviation (MBD) using MCD43C3.v61 as the reference (deviation = product -MCD43C3v6.1). MCD43C3v6.1 was selected as the reference based on the results of a previous comparison of MCD43C3 v6/6.1, GLASS-AVHRR, and C3S-v1/v2 against spatially representative in situ measurements [34]. MCD43C3 v6.1 had the best performance over snow and the most stable performance with different snow cover conditions. The bias was calculated in two periods to evaluate all the C3S sensors: (2001)(2002)(2003)(2004)(2005) including MCD43C3 v6/v6.1, CLARA-A2.1, GLASS-AVHRR/MODIS, C3S-AVHRR v1/v2, and C3S-SPOT v1/v2 (2015-2019) including C3S-PROBA v1/v2. In each period, the monthly bias was calculated separately for the snow-free and snow seasons, and within each season, the bias was analyzed over the different IGBP land cover types.
Stability was assessed based on the decadal trends of the annual bias during: 1) the MODIS era (2001-2019) and 2) the AVHRR era (1982-2019). As in the bias assessment, MCD43C3 v6.1 was used as the reference in the MODIS era. This choice is supported by the good temporal stability of MCD43C3 at the stations from Urraca et al. [34] providing long-term albedo measurements over snow (see the Supplementary Material). In the AVHRR era, we selected CLARA-A2.1 as the reference because it has the snow albedo values and trends most consistent with MCD43C3 v6.1 during the MODIS era (see Section IV). The decadal trends in absolute albedo were also calculated to benchmark the magnitude of the bias trends. Compared with the bias assessment, trends cannot be calculated separately for snow-covered and snowfree months because their magnitude would depend on the snow season duration in each pixel. The percentage of pixels where the bias trends are within the GCOS stability limits, i.e., max(1%, 0.001) per decade, was also calculated. The temporal acceleration or deceleration of the trends was analyzed with a ten-year moving window. C3S albedo trends were calculated by concatenating the sensors as follows: AVHRR (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998), SPOT (1999SPOT ( -2013, and PROBA-V (2014-2019). C3S v1 was excluded from the stability analysis due to the lack of snow information during the AVHRR period.
In each month, bias and stability were evaluated in those pixels where all the products have more than 90% of valid retrievals during the study period. This reduces the number of estimations validated at polar regions because CLARA-A2.1 does not process observations with SZA > 70 • and C3S albedo is systematically missing above 75 • N. However, this procedure removes artificial differences due to the spatio-temporal differences in the number of missing values, guaranteeing that the results are comparable between products.
Completeness, bias, and stability were evaluated both spatially and temporally. The temporal analysis was made by analyzing the annual variation in completeness, bias, and stability metrics over five snow-covered regions: northern hemisphere (NH) above 65 • N (to isolate the effects of low solar elevation and polar night), NH below 65 • N, Southern Hemisphere, Greenland, and Antarctica. The spatial aggregation was made taking into account the area covered by each pixel.

A. Coverage and Completeness
The completeness of most products decreases during months with low solar elevation due to the fewer observations available to retrieve albedo and the higher uncertainty of these retrievals. The number of gaps and the latitude bands affected vary in each product depending on their retrieval algorithm (see Fig. 1). The temporal variation in the number of valid retrievals is analyzed in Fig. 2 to detect potential trends/changes that could affect product stability.
GLASS gap-fills pixels with missing or low-quality satellite measurements, so both GLASS-AVHRR (whole record) and GLASS-MODIS (2000-2010) are gap-free. However, GLASS-MODIS was not gap-filled after 2011 leading to a sudden increase in the missing values at high latitudes (see Fig. 2). Excluding GLASS, MCD43C3 shows the largest number of valid retrievals globally. MCD43C3 has global coverage with a percentage of valid values of around 80% above 65 • N, Greenland, and Antarctica. v6 already had the best completeness, but the new v6.1 slightly extends the percentage of valid retrievals by around 5% above 65 • N.
CLARA-A2.1 also has global coverage but presents more missing values than MCD43C3 because it only processes images with SZA < 70 • , reducing the number of available observations in polar regions. This reduces completeness over 65 • N (∼63%), particularly in DJF and SON, and also below 65 • N in DJF (∼53%). CLARA-A2.1 completeness also decreases in Greenland (∼64%) and Antarctica (∼55%). Besides, the percentage of missing values is not stable over time, particularly in those regions/seasons where observations are made at SZAs close to CLARA's processing limit (see Fig. 2). These instabilities are due to drifts in the satellite daytime pass caused by the lack of active control to maintain a sun-synchronous orbit [15]. Both morning and afternoon NOAA satellites progressively drift away from solar noon worsening the illumination conditions. Stability improves in 2004 after the launch of an additional mid-morning satellite (NOAA-17), and in 2008 after the launch of Metop series, which compared with NOAA's satellites are able to maintain a sun-synchronous orbit. However, the number of missing values is increasing again since 2015 due to the drift of NOAA's afternoon satellites (NOAA-18, NOAA-19).
C3S albedo covers only from 60 • S to 80 • N, but values above 75 • N are generally missing, particularly during the SPOT period. Below 65 • N, completeness is slightly superior to that of MCD43 products. However, C3S completeness at high latitudes is not consistent between satellite generations. C3S v1 completeness above 65 • N differs for each sensor: AVHRR (∼75%), SPOT (∼62%), and PROBA-V (∼85%). C3S v2 improves the multisensor consistency, having similar completeness during both the AVHRR and PROBA periods (∼75%), but completeness decreases again during the SPOT period (∼62%). Besides, some months are fully missing in both the versions during the transition between NOAA-7 and -9, NOAA-9 and -11, and NOAA-11 and -14 due to the quality degradation of AVHRR observations and the lack of calibration coefficients [55].
MIL3MLSN has the largest number of missing values with only 45% of valid retrievals above 65 • N, decreasing below 25% during snow-covered months. Note that Greenland and Antarctica were already excluded due to the typically poor quality of MISR aerosol retrievals in these regions [19]. Compared with other products, MIL3MLSN missing values are not latitude-dependent but increase in particular regions such as the Gulf of Guinea or the Amazon basin, and the Northern part of America and Eurasia. A low coverage over snow was also observed by studies using the Level 2 product [25]. Due to the low coverage of the product over snow-covered pixels, MIL3MLSN was removed from the bias and stability assessment.

B. Bias Assessment
The consistency of the snow albedo products was evaluated using MCD43C3 v6.1 as the reference. The bias maps were calculated separately for snow-covered (see Fig. 3) and snowfree (see Fig. S5) months, and per IGBP land cover type (see Figs. 4 and S6). The bias was further assessed by plotting the temporal variations in albedo during the whole period covered by each product (see Fig. 5).
The differences between the two MCD43C3 versions (v6 and v6.1) are negligible globally both during snow-covered and snow-free months. CLARA-A2.1 has the smallest bias during both snow-covered and snow-free months. CLARA-A2.1 snow albedo is consistent with that of MCD43C3 v6.1 both spatially and per land cover type. It shows a small negative bias globally of −0.019 in most land cover types, particularly over croplands (Great Plains) and snow/ice (Greenland, Antarctica). No systematic difference is observed between CLARA-A2.1 and MCD43C3 v6.1 during snow-free months (see Fig. S5).
GLASS-MODIS has a moderate negative bias of −0.031 during the snow season with respect to MCD43C3 6.1, driven by a negative bias over snow/ice land cover types, mainly in Greenland and Antarctica. The differences between both the products are negligible during the snow-free season. The bias is much larger in GLASS-AVHRR, which has a large positive bias over regions with seasonal snow around +0.1 to +0.2 and a negative bias around −0.1 over permanently snow-covered regions. Both the effects compensate globally leading to an artificially low mean bias of +0.019 during the snow season (see Fig. 3). The positive bias in regions with seasonal snow appears in all the land cover types, exceeding +0.2 in regions such as Eastern Europe, the Kazakh Steppe, and the Great Plains. The positive bias is larger during the snow melt (see Fig. 5): JJA above 65 • N (+0.15) and MAM below 65 • N (+0.1). A positive bias is also observed during the snow-free season globally (+0.056). On the contrary, and similar to GLASS-MODIS, GLASS-AVHRR underestimates albedo over snow/ice by around −0.08.
The quality of C3S snow albedo improves from v1 to v2. The addition of snow information during the AVHRR period mitigates the underestimation of albedo over snow from −0.163 to −0.038. Snow albedo also improves during the PROBA period correcting the strong variations in C3S PROBA v1 bias with land cover type: a positive bias over forests (e.g., Canada's boreal forest) and a negative  latitudes. The underestimation is also larger in low VGT land cover types (savannas, grasslands, croplands) than over forests. Seasonally, the underestimation is more evident in the middle of the snow season: MAM above 65 • N and DJF below 65 • N. C3S also underestimates snow albedo by about −0.2 in the small part of Greenland covered by the product. The Andes is the only snow-covered region where C3S overestimates snow albedo, particularly during SON. On the other hand, C3S v2 systematically overestimates snow-free albedo compared with MCD43C3 v6.1 (+0.021 to +0.024) in all the land cover types and with the three satellite sensors used.
Despite the improved multisensor consistency, Fig. 5 reveals that the different periods of the C3S multisensor dataset are still not fully consistent. Several discontinuities/jumps appear in the transition between sensors. Discontinuities are larger over regions and periods with more snow, particularly MAM above 65 • N and DJF below 65 • N. Above 65 • N, snow albedo increases from AVHRR to SPOT. Below 65 • N, albedo increases from AVHRR to SPOT, and then decreases from SPOT to PROBA. The discontinuity in the transition between AVHRR and SPOT is measurable in the overlapping years (2000)(2001)(2002)(2003)(2004)(2005), with an annual mean difference of +0.024 above 65 • N. Discontinuities also appear during the AVHRR period in the transition between NOAA satellites, with the albedo increasing from NOAA-7 to -9 (1985) and decreasing again from NOAA-11 to -14 (1994).

C. Stability Assessment
Stability is analyzed by comparing the albedo bias trends and the absolute albedo trends of each product (see Figs. 6-8). The temporal change in the trends is evaluated with a ten-year window function (see Fig. 7).
CLARA-A2.1 trends are those most consistent with MCD43C3 v6.1 with 55.6% of the snow-covered pixels within the GCOS stability requirements. CLARA-A2.1 trends are slightly more negative than those of MCD43C3 v6.1 but both the products are very consistent spatially during the past 20 years. Both show negative trends of around [−0.02, −0.03 per decade] in Siberia, Alaska, and Eastern Europe. On the contrary, both show local albedo gains in Kazakhstan (+0.03 per decade), the Chukchi Peninsula, and Eastern Greenland. Trends of both the products are also temporally consistent over regions with seasonal snow (see Fig. 7). Both gather the same trend dynamics in the NH and the progressive decrease in albedo in the Andes. Trends are also consistent over Greenland, with both the products showing an albedo increase during the past decade. However,   (1981-2000), which can be specially observed above 65 • N in MAM (see Fig. 5). GLASS-MODIS has better stability but has trends inconsistent with those of the other products in all the regions. Particularly, both the GLASS products have a strong positive trend in Antarctica and Greenland after 2010 not shown by any other product.
C3S-v2 shows the trends least consistent with MCD43C3 v6.1. The effects of the discontinuities during the transition between different sensors are visible both temporally and spatially. C3S-v2 has only a 24.1% of snow-covered pixels within the stability requirements during 2001-2019. C3S-v2 trends are generally more negative than MCD43C3 ones, particularly over regions with low VGT such as open shrublands (high latitude), grasslands, and croplands (e.g., the Great Plains and the Kazakh Steppe). Temporally, C3S-v2 trends are only consistent with MCD43C3 below 65 • N and during the SPOT-PROBA period (see Fig. 7). However, the sign of C3S v2 trends differs from MCD43C3 during the same period above 65 • N and in the Andes. Besides, temporal inconsistencies aggravate during the AVHRR period in all the regions. The ten-year window trends show very sharp changes driven by the discontinuities in 1985 (NOAA 7-9), 1994 (NOAA [11][12][13][14], and 2000 (NOAA to SPOT).
Trends are also analyzed spatially during the whole AVHRR era (1982-2019) using CLARA-A2.1 as the reference (see Fig. 8). Similar to the 2001-2019 period, GLASS is the closest product to CLARA (GCOS requirement = 52.5%), and C3S-v2 shows the most inconsistent trends again (GCOS requirement = 28.9%). During this period, CLARA-A2.1 trends increase compared with those observed during the MODIS era, with most regions showing a global decrease in albedo, particularly the Greenland coast, Siberia, Kazakhstan, and Alaska. Contrary, an albedo increase is only observed in the Great Plains and Eastern Russia. Note that to assess the real influence of these trends on the surface energy budget, annual trends need to be calculated by weighting monthly albedos with the average shortwave incoming irradiance of each pixel and month. This is not done throughout this article, but we include the trends in annual effective albedo as in the Supplementary Material (see Fig. S7). These trends show generally similar spatial patterns, but their magnitude is generally smaller. A detailed analysis of the effects of albedo trends on the energy budget would also require quantifying the influence of missing values on the annual effective albedo. A cross-consistency check between albedo and snow cover trends would also be very useful, due to the higher number of snow stations available to evaluate stability during 1982-2020.
V. DISCUSSION As discussed in Urraca et al. [34], the differences between products originate from the quality of the optical measurements and multiple factors regarding their retrieval algorithms such as radiance calibration, snow/cloud screening, atmospheric correction, and hemispherical integration, among others. The differences are also explained by the different spatio-temporal resolutions of the products. A high spatial resolution helps represent scattered snow, whereas a high temporal resolution is essential during the snow melting and onset periods. Some differences also exist in the shortwave spectral range covered by each product and the spectral responses of their satellite sensors, but the magnitude of the deviations introduced in the intercomparison by these effects is much smaller than those observed between products. The percentage of the shortwave range covered by each product oscillates from 98.58% (GLASS-AVHRR, 300-4000 nm) to 99.51% (CLARA-A2.1 250-2500 nm) in terms of downwelling irradiance. The differences in terms of albedo, which is a ratio, should be even smaller. Liang [56] proved that the narrowbands from all the sensors used in this study are equally valid predictors of the shortwave albedo, with the differences between them being below 0.001 (interquartile range). Thus, for intercomparison, we can assume that the shortwave black-sky albedo values provided by each product are physically consistent.
A previous comparison of climate albedo products against 11 spatially representative in situ measurements [34] showed that MCD43C3 v6.1 has the smallest bias over snow and the most stable bias with different snow cover conditions. We extended the previous comparison in the Supplementary Material by evaluating the temporal stability of the products at the stations (six out of 11) providing at least ten years of measurements during the MODIS era. We analyze the temporal stability of the annual bias (satellite-in situ) to identify potential artificial trends in the satellite products affecting our intercomparison (see Figs. S1-S3). MCD43C3 is again the product showing the best temporal stability, with trends in the bias below 0.01 per decade in four out of the six stations, and the smallest standard deviation of monthly differences. A statistically significant trend in the bias of +0.03 per decade is only observed at SXF. As mentioned in the introduction, a comprehensive assessment of the temporal stability of satellite products against in situ measurements is hindered by the few sites available with long-term gap-free albedo measurements that are spatially representative over snow. Despite this, the results obtained at the six long-term sites available suggest that MCD43C3 is the most suitable choice as a reference to evaluate the bias and the temporal stability of satellite products in a global intercomparison. Note that throughout this article, terms such as bias or under/overestimation always refer to differences with respect to MCD43C3. We acknowledge that MCD43C3 has its own bias and trends in the bias, which could be relevant in some regions where MCD43C3 quality degrades. However, the comparison against in situ data shows that the magnitude of these effects is generally low and much lower than the differences between satellite products or the annual albedo trends. Indeed, the results obtained by C3S-v2 and GLASS using MCD43C3 v6.1 as a reference are similar to those obtained in the comparison against in situ data, supporting the selection of MCD43C3 v6.1 as the reference for this intercomparison.
The good performance of MCD43 albedo products over snow was also observed in other studies [25], [27], [28], [57]. The higher temporal resolution of MCD43 albedo products (one day) compared with CLARA-A2.1 (five days), GLASS (eight days), or C3S-v1/v2 (∼ten days) may be one of the main factors explaining its superior performance over snow. As in the comparison against in situ data, no significant differences were found between v6.1 and v6 in terms of bias and stability, but v6.1 extends the percentage of valid retrievals by around 5% in polar regions.
CLARA-A2.1 snow albedo values and trends are those most consistent with MCD43C3 v6.1 globally over different land cover types. Their only discrepancies are a small negative bias in CLARA of around −0.021 globally over snow, and a progressive decrease in CLARA-A2.1 snow albedo globally since 2015. Karlsson et al. [58] reported a similar anomalous decrease after 2016, arguing that it could be caused either by natural reasons or by unaccounted drifts in AVHRR calibration. The dataset used in this study is composed of the original CLARA-A2 (January 1982-December 2015), the A2.1 extension (January 2016-June 2019), and the ICDR extension (July 2019-December 2020). The anomalous negative trend starts with the A2.1 extension and is not observed in any other product, so the most likely explanation is the presence of an artificial trend during the A2.1 extension. While all the CLARA datasets share the same algorithm and preprocessing, AVHRR calibrations were not fully recalculated for the A2.1 extension but temporally extended from 2016 onward (Aku Riihelä, personal communication, 2022). Orbital drifts in NOAA afternoon satellites are also progressively increasing the SZA of CLARA albedo retrievals. The cutoff of SZA = 70 • mitigates the effects of orbital drifts at polar regions, at the expense of reducing the product completeness. However, both the orbital drifts and the extension of AVHRR calibrations since 2015 are the most likely explanation for the anomalous negative trend during 2016-2020.
The main limitation of CLARA is its large number of missing values over snow due to CLARA's cutoff limit at SZA = 70 • . Besides, the number of missing values changes during the period when only one NOAA satellite is available due to the orbital drifts. All the pixels affected by these temporal variations must be removed to monitor the albedo trends, reducing, even more, the coverage of the product at high latitudes. The effects of these gaps in applications related to the energy radiation budget should be low, due to the low incoming irradiance during those months, but they have to be quantified.
Among the two GLASS versions, GLASS-MODIS is the one most consistent with MCD43C3 6.1, particularly during the snow-free season. Its main limitation is a sudden increase in missing values in 2010, going from gap-free to having a similar number of missing values to CLARA-A2.1. Besides, both the GLASS versions underestimate albedo compared with the other products in regions with permanent snow, particularly Antarctica and Greenland. In addition, GLASS-AVHRR bias varies strongly with the snow cover conditions. It underestimates albedo over fully snow covered regions, as also observed in the comparison against in situ data at Antarctic stations and in most stations during fully snowcovered months. On the contrary, it has the largest positive bias during the snow season in regions with seasonal snow caused by an overestimation of albedo above +0.2 during the melting season: NH in JJA above 65 • N and MAM below 65 • N (see Fig. 5). This issue was also observed in the in situ comparison and it is likely caused by an artificial delay of snowmelt by the GLASS snow mask.
GLASS trends from both the products are neither spatially nor temporally consistent with those of CLARA-A2.1 and MCD43C3 v6.1. GLASS-AVHRR has larger trends and more abrupt changes. Larger albedo trends over snow could be due to the larger number of pixels classified as snow-covered. Abrupt temporal changes are more frequent from 1982 to 2004, when fewer AVHRR observations are available, and in polar regions, which are those most affected by AVHRR orbital drifts, so they are likely caused by trying to gap-fill missing AVHRR observations. GLASS-MODIS does not show these oscillations but its trends are inconsistent with MCD43C3 and CLARA, especially in Antarctica and Greenland.
C3S multisensor product has improved its performance over snow from v1 to v2, particularly during the AVHRR period, with the addition of snow information, and during the PROBA-V period. However, as observed in the in situ comparison, the performance during the three periods (AVHRR, SPOT, PROBA) strongly varies with the snow cover conditions. C3S v2 has the largest negative albedo bias over the snow season (−0.029 to −0.040) and a positive overestimation in the snow-free season (+0.021 to +0.024). Besides, despite the improved multisensor consistency, the product stability is insufficient to monitor accurately the snow albedo trends. Some instabilities were already observed in the comparison against in situ data during the AVHRR-SPOT transition, but the intercomparison confirms that these discontinuities occur globally and in most sensor transitions, being more evident during the AVHRR period in the transition between different NOAA sensors, and in the transition from AVHRR to SPOT. Sanchez-Zapero et al. [55] reported a positive bias of 17.7%, −2.5%, −4.4%, +3.9%, +3.4%, and +1.2% for NOAA-7, 9, 11, 14, 16, and 17 AVHRR sensors with respect to SPOT-VGT, respectively, over 20 desert calibration sites. The effects are larger over snow, due to the larger albedo values and additional inconsistencies in the snow mask. Sanchez-Zapero et al. [55] acknowledged that the magnitude of the snow-covered area is not consistent between sensors, with SPOT having the largest snow-covered area followed by PROBA-V, NOAA- [16/17], and NOAA [7][8][9][10][11][12][13][14]. The latter period had a significantly smaller snow area due to the lack of band 3a in AVHRR-2 sensors [55].
AVHRR provides the longest dataset for monitoring surface albedo globally. However, the study evidences the challenges of monitoring snow albedo trends with AVHRR data. Only CLARA-A2.1 is able to produce acceptable estimations in terms of bias, while all the products face stability issues, particularly the C3S dataset, leading to inconsistent snow albedo trends between products temporally and spatially. AVHRR instabilities are driven by orbital drifts in NOAA satellites that progressively move the overpass time away from the local solar noon, worsening the illumination conditions for retrieving albedo. This also affects the calibration of visible bands, which currently relies on vicarious calibration methods. Besides, the number and type (morning, mid-morning, afternoon) of satellites available are not constant since 1982.
The largest uncertainty appears before 1992 when only one afternoon satellite is available. Morning and afternoon satellites were simultaneously available during 1992-2004, and a mid-morning satellite was added in 2004. The most stable period comes after the launch of Metop satellites in 2008, which are able to maintain a sun-synchronous orbit. There are also small differences in the spectral bands and responses of AVHRR-2 and AVHRR-3, with the latter including an additional visible channel (1580-1640 nm) during daytime that is exploited by some products.
The stability of AVHRR-based products depends on the satellites and bands used, their calibration, and the retrieval algorithm. The best stability is obtained by CLARA using channels 1 and 2 from all the satellites simultaneously available. AVHRR radiances are calibrated with MODIS data during 2000-2014 and with a vicarious method before 2000. Therefore, uncertainty increases before 2000, particularly before 1992 when only one satellite is simultaneously available. CLARA-A2.1 validation report [58] shows the increased uncertainty before 2004 using Greenlandic stations as the reference, but no systematic trend was observed. The anomalous negative trend since 2015 observed in this study is most likely due to extending AVHRR calibrations during the A2.1 extension. In contrast, C3S uses only one satellite at a time (discarding NOAA-12, 15 morning satellites) and exploits the additional band 3a of AVHRR/3 when available. This less consistent approach introduces discontinuities in the C3S AVHRR record. The accuracy of C3S snow albedo is also lower despite using a more complex BRDF-based algorithm. CLARA has better accuracy by just averaging TOA bidirectional reflectances and by discarding all the observations with SZA > 70 • . This suggests that when exploiting AVHRR data, consistency should always be the first priority, as this may not only lead to better stability but also smaller bias. Otherwise, complex retrieval algorithms become useless. The comparison also shows that discarding observations made at extreme solar angles (CLARA) instead of using them (C3S) or filling them (GLASS) leads to more stable estimates because the orbital drifts in NOAA satellites complicate even more the albedo retrieval under these conditions.

VI. CONCLUSION
We evaluate the consistency of snow albedo values and trends of MCD43C3 v6/v6.1, CLARA-A2.1, GLASS-AVHRR, GLASS-MODIS, and C3S v1/v2 multisensor product globally. MCD43C3 v6.1 was used as the reference due to its superior performance over snow when compared with in situ measurements in a previous study. CLARA-A2.1 is the only product exploiting the 40-year AVHRR dataset (1982-2020) with low bias and acceptable stability over snow. A large number of missing values of MIL3MLSN albedo, and the low accuracy and stability of both GLASS (particularly GLASS-AVHRR) and C3S-v2 albedo over snow, made them inappropriate to monitor snow albedo trends. The bias of both GLASS-AVHRR and C3S-v2 changes with the snow cover conditions, underestimating and overestimating albedo over snow-covered and snow-free surfaces, respectively. GLASS-AVHRR overestimation particularly increases during snow melt likely due to an artificially extended snow season. Besides, GLASS-AVHRR trends are larger and more unstable than those of MCD43C3. GLASS-MODIS has a good consistency with MCD43C3 during the snow-free period and snow-covered mid-latitudes, but it has a sudden increase in missing values after 2011. Besides, both the GLASS products have a negative bias over regions with permanent snow and trends inconsistent with those from other satellite products.
The C3S multisensor dataset has improved over snow from v1 to v2, but the improvements made are still insufficient. C3S v2 systematically underestimates albedo over snow-covered regions. Besides, discontinuities in the transition between satellites and sensors, aggravated by discontinuities in their snow masks, lead to spatially and temporally unstable snow albedo trends.
The study evidences the challenges faced by AVHRR-based products to provide stable albedo products over snow. Snow albedo values and trends are inconsistent between the three AVHRR-based products evaluated and only CLARA has low bias and acceptable stability. On the satellite community side, consistency should always be a priority when developing AVHRR products. On the user side, assessing the consequences of these instabilities in their applications becomes mandatory. The uncertainty of AVHRR-based products increases even more before 2000, due to the fewer satellites available, and also the fewer stations and products available to validate or intercompare AVHRR-based products. Assessing the cross-consistency with other variables (e.g., snow cover extent) could help in understanding better their stability during this period. Ruben Urraca received the Ph.D. degree in data science from the University of La Rioja, Logroño, Spain, in 2018.
He worked on reducing the uncertainty of solar radiation measurements (in situ, satellite, reanalysis) for PV applications, releasing a new satellite-based quality control method for solar radiation measurements (www.bqcmethod.com). Since 2019, he has been with the Joint Research Centre, European Commission, Ispra, Italy. He worked in the PVGIS service, using satellite data to estimate the PV potential, and in the Covenant of Mayors, evaluating the effectiveness of local action plans to mitigate and adapt to climate change. He currently works for the Copernicus Program, improving the validation protocols of satellite products by investigating topics, such as the spatial representativeness of in situ measurements, the temporal stability of climate data records, and the uncertainty budget closure. His main research interests include the study of surface energy balance.
Christian Lanconelli received the Ph.D. degree in atmospheric physics from the University of Ferrara, Ferrara, Italy, in 2007, exploring the effects of Earth surface reflectance properties on aerosol direct radiative forcing by combining modeling and measurement approaches.
He worked at the Institute of Atmospheric Sciences and Climate (ISAC), Italian National Research Council, Rome, Italy, from 2002 to 2015, developing experience in radiative transfer modeling and solar radiation measurement, aimed to assess the effects of aerosols and clouds on Earth's energy balance. Since 2015, he has been with the Joint Research Centre, European Commission, Ispra, Italy, where he worked on the Quality Assurance for Essential Climate Variables and Copernicus supported projects, by developing competencies in Monte Carlo 3-D radiative transfer modeling, mainly to assist the fitness for purpose evaluation of satellite-derived biogeophysical variables and surface reflectance. He has also been the Project Manager of the Baseline Surface Radiation Network since 2019.
Fabrizio Cappucci received the master's degree in physics from the Physics Department "Aldo Pontremoli," University of Milan, Milan, Italy, in 2009.
He is currently with the Directorate-D, Nature Conservation and Observation Unit, Joint Research Centre, European Commission, Ispra, Italy. His main research interests include development of advance techniques for assessing the fitness for the purpose of remote sensing data, especially within the Copernicus Program, to implement a rigorous control chain to assess Earth observation products' quality and reliability.
She was a Research Fellow at the European Space Agency (ESA) from 1997 to 1998 and a Grant-Holder at the Space Applications Institute, Joint Research Centre, European Commission, Ispra, Italy, from 1998 to 2000. She is the Project Leader at the Directorate Sustainable Resources Unit D.6 Nature Conservation and Observations, Ispra. She works for Copernicus Land Service and supports the Knowledge Center on Earth Observation (KCEO). She has authored or coauthored more than 100 refereed publications. Her main research interests include research on the theory of radiation transfer in plant canopies and, more generally, the development of tools to quantitatively interpret satellite remote sensing data in the optical spectral domain via inverse methods.
Dr. Gobron is a member of the Copernicus Sentinel-3 OLCI Mission Advisory Group. She received the Group Achievement Award by NASA in 2001 and the Best Young JRC Scientific in 2004.