Retrieval and Evaluation of Global Surface Albedo Based on AVHRR GAC Data of the Last 40 Years

Abstract

In this study, the global land surface albedo namely GAC43 was retrieved for the years 1979 to 2020 using Advanced Very High Resolution Radiometer (AVHRR) global area coverage (GAC) data onboard National Oceanic and Atmospheric Administration (NOAA) and Meteorological Operational (MetOp) satellites. We provide a comprehensive retrieval process of the GAC43 albedo, followed by a comprehensive assessment against in situ measurements and three widely used satellite-based albedo products, the third edition of the CM SAF cLoud, Albedo and surface RAdiation (CLARA-A3), the Copernicus Climate Change Service (C3S) albedo product, and MODIS BRDF/albedo product (MCD43). Our quantitative evaluations indicate that GAC43 demonstrates the best stability, with a linear trend of ±0.002 per decade at nearly all pseudo invariant calibration sites (PICS) from 1982 to 2020. In contrast, CLARA-A3 exhibits significant noise before the 2000s due to the limited availability of observations, while C3S shows substantial biases during the same period due to imperfect sensors intercalibrations. Extensive validation at globally distributed homogeneous sites shows that GAC43 has comparable accuracy to C3S, with an overall RMSE of approximately 0.03, but a smaller positive bias of 0.012. Comparatively, MCD43C3 shows the lowest RMSE (~0.023) and minimal bias, while CLARA-A3 displays the highest RMSE (~0.042) and bias (0.02). Furthermore, GAC43, CLARA-A3, and C3S exhibit overestimation in forests, with positive biases exceeding 0.023 and RMSEs of at least 0.028. In contrast, MCD43C3 shows negligible bias and a smaller RMSE of 0.015. For grasslands and shrublands, GAC43 and MCD43C3 demonstrate comparable estimation uncertainties of approximately 0.023, with close positive biases near 0.09, whereas C3S and CLARA-A3 exhibit higher RMSEs and biases exceeding 0.032 and 0.022, respectively. All four albedo products show significant RMSEs around 0.035 over croplands but achieve the highest estimation accuracy better than 0.020 over deserts. It is worth noting that significant biases are typically attributed to insufficient spatial representativeness of the measurement sites. Globally, GAC43 and C3S exhibit similar spatial distribution patterns across most land surface conditions, including an overestimation compared to MCD43C3 and an underestimation compared to CLARA-A3 in forested areas. In addition, GAC43, C3S, and CLARA-A3 estimate higher albedo values than MCD43C3 in low-vegetation regions, such as croplands, grasslands, savannas, and woody savannas. Besides the fact that the new GAC43 product shows the best stability covering the last 40 years, one has to consider the higher proportion of backup inversions before 2000. Overall, GAC43 offers a promising long-term and consistent albedo with good accuracy for future studies such as global climate change, energy balance, and land management policy.

1. Introduction

Land surface albedo, defined as the ratio of reflected shortwave radiation to incoming solar radiation, represents the integrated surface hemispherical reflectivity across the solar spectrum. It is a critical variable linking the land surface to the climate system by regulating shortwave energy exchanges [1]. In particular, seasonal vegetation phenology can significantly alter surface albedo, which in turn influences the global energy balance and potential albedo-related warming or cooling feedback. Therefore, long-term, accurate albedo data are essential for both regional and global climate models. These models require surface albedo measurements with an absolute accuracy of 0.02–0.05 units for snow-free and snow-covered land across various spatial scales (ranging from 10s of meters to 5–30 km) and temporal scales (from daily to monthly) [2,3]. Consequently, the Global Climate Observing System (GCOS) has designated albedo as an essential climate variable [4,5]. However, the surface albedo varies spatially and temporally due to both natural processes (e.g., solar illumination, snow fall, and vegetation growth) and human activities (e.g., the clearing and replanting of forests, the sowing and harvesting of crops, and the burning and grazing of rangelands). Thus, a consistent, accurate, and global multi-decadal surface albedo dataset is crucial for climate dynamics [6,7], surface energy budget [8], and possible carbon sequestrations [9].

Currently, global land surface albedo products have been extensively developed using polar-orbiting satellite observations, benefitting from their global coverage, regular revisit schedules, and long-term maintenance [10]. These include the Terra and Aqua/Moderate Resolution Imaging Spectroradiometer (MODIS) [11], Suomi NPP/Visible Infrared Imaging Radiometer Suite (VIIRS) [12], National Oceanic and Atmospheric Administration (NOAA) and Meteorological Operational (MetOp)/Advanced Very High Resolution Radiometer (AVHRR) [13], SPOT/VEGETATION (VGT) [14], PROBA/VGT [15], Sentinel-3/Ocean and Land Colour Instrument (OLCI), and Sea and Land Surface Temperature Radiometer (SLSTR) [16]. These albedo products provide data spanning different periods, with spatial resolutions ranging from 500 m to 25 km and temporal resolutions from daily to monthly intervals. Although MODIS surface albedo products are widely used, the Terra and Aqua satellites have been operating well beyond their expected lifespans (launched in 1999 and 2002, respectively), limiting long-term analysis to approximately 20 years. However, climate studies require multidecadal observations [13]. So far, only two global satellite-based albedo products offer 40 years of data. One is developed under the Copernicus Climate Change Service (C3S) by the European Union, offering a 10 day temporal resolution and various spatial resolutions, which depend on the sensors used. This dataset includes data from NOAA/AVHRR, SPOT/VGT, PROBA-V/VGT, and potentially data from Sentinel-3/OLCI and SLSTR in the future. A critical spectral harmonization step is performed to homogenize data from these sensors by accounting for differences in spectral response functions [15]. While several studies have evaluated the performance and consistency of the C3S albedo based on SPOT/VGT, PROBA-V/VGT, and Sentinal-3/OLCI and SLSTR, few studies have yet evaluated the consistency of the AVHRR-based albedo compared to other albedo products within the C3S framework. Such consistency is crucial for long-term research dating back to the 1980s. The second long-term product is the third edition of the Satellite Application Facility on Climate Monitoring (CM SAF) CLoud, Albedo and surface RAdiation (CLARA-A3) dataset [13,17], which offers temporal resolutions at pentad or monthly intervals, with a spatial resolution of 25 km. CLARA-A3 is expected to offer better data consistency than C3S, as it relies solely on AVHRR data from the NOAA and MetOp satellite series from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), despite the discontinuation of the AVHRR/1 sensor in approximately 1991 and the transition from AVHRR/2 to AVHRR/3 sensors around 2000. However, CLARA-A3 does not provide albedo values under fixed illumination and observation geometry; instead, it offers average albedo values under mean instantaneous observation conditions within a pentad or month. This limitation complicates long-term surface albedo monitoring due to surface anisotropy caused by the bidirectional reflectance distribution function (BRDF [18]).

The BRDF describes the surface anisotropy and the corresponding BRDF parameters can be retrieved from sufficiently well-distributed multi-angle satellite observations [11]. Land surface albedo is then derived by integrating the BRDF parameters over all hemispherical viewing and illumination angles. Among the BRDF models, the most widely used kernel-based model known as RossThick-LiSparseR (RTLSR) [19,20] is frequently adopted for developing albedo products from various multi-angular space-borne sensors. This model consists of three kernels that form a linear combination to represent surface reflectance with corresponding illumination and observation geometries: isotropic kernel, volumetric kernel, and geometric kernel. The isotropic kernel accounts for isotropic scattering, the volumetric kernel represents the scattering properties of a turbid medium, and the geometric kernel captures the shadowing effects of sparse vegetation [21]. It is worth mentioning that the success of BRDF models also depends significantly on the revisit frequency and variability of angular sampling provided by the sensors. Increasing the number of observations and enhancing angular sampling diversity through multi-sensor integration significantly improves the robustness of BRDF retrievals compared to using data from a single sensor [22,23]. This is because angular observations from individual sensors are constrained by their scanning configurations and the orbital characteristics of their platforms [24]. However, robust and accurate BRDF retrievals under snow-covered conditions remain challenging with the aforementioned kernels due to the wide variety of snow types, each with distinct reflectance anisotropy characteristics [25]. Therefore, this study focuses on global snow-free surface albedo estimates and evaluations.

The AVHRR sensor, a multispectral imaging sensor aboard the polar-orbiting NOAA and MetOp satellites, offers a unique global remote sensing dataset spanning from the 1980s to the present. Several platforms operate simultaneously, capturing observations under varying geometric conditions, which provide a wide range of viewing and solar angle combinations. This variety enhances the effectiveness of the BRDF model and improves the quality of BRDF retrievals. Consequently, this study utilizes the standard MODIS BRDF/albedo algorithm (RTLSR) to AVHRR global area coverage (GAC) data, aiming to generate a global BRDF/albedo climate data record from 1979 to 2020.

The primary objective of this study is to develop a novel framework for retrieving global surface albedo from AVHRR GAC data, focusing primarily on snow-free conditions over the past four decades. Additionally, the framework is applied to process a long-term data series, generating a global albedo product, which is evaluated using in situ measurements and comparisons with other satellite-based products. The paper is organized as follows; Section 2 introduces the data used for retrieving surface albedo, as well as the datasets validated at in situ and intercomparison with other satellite-based albedo products; Section 3 presents the algorithm for albedo retrieval; Section 4 gives the results and preliminary evaluations; and Section 5 and Section 6 provide discussion and conclusions, respectively.

2. Data and Pre-Processing

2.1. AVHRR GAC

In this study, AVHRR GAC data from NOAA and MetOp polar-orbiting satellites, spanning the period from 1979 to 2020, were utilized to generate global surface albedo. The AVHRR instrument, first launched in 1978 aboard the TIROS-N satellite from NOAA, has been continuously deployed on all subsequent NOAA satellites, from NOAA-6 to NOAA-19, as well as on the three MetOp satellites, thus providing a more than 40 year data record. The nominal 3 km × 4.4 km GAC data represent reduced-resolution imagery created onboard the satellite by sampling one line for every three and averaging four of every five adjacent samples along each scan line [26]. However, in this study, we utilize its resampled climate model grid size of 5 km × 5 km. This resampling process further reduces the original local area coverage (LAC) data with a nominal spatial resolution of 1.1 km × 1.1 km.

Notably, all GAC data were processed using the PyGAC V2 software [27], ensuring complete calibration for reflectance and radiance, inter-sensor consistent calibration [28], and accurate geometric information [26], such as solar zenith angle (SZA), view zenith angle (VZA), and relative azimuth angle (RAA, defined by the absolute difference between the solar azimuth angle and the view azimuth angle). Each satellite transmits GAC data twice daily, corresponding to ascending (crossing the equator from the south) and descending (crossing the equator from the north) nodes. It is worth mentioning that the pixel with the lowest VZA is selected during resampling when swaths overlap.

Figure 1 provides an overview of the satellites equipped with the AVHRR instrument used to generate the AVHRR GAC albedo (GAC43, following the naming convention of MODIS albedo product, namely MCD43) data throughout the study period. This figure shows that almost all AVHRR sensors experienced significant orbit drift, as indicated by local time (Figure 1a) and SZA changes (Figure 1b) over their lifespan, except for those onboard the MetOp satellites. Additionally, data availability varied depending on the number of operational platforms in orbit. For instance, during the 1980s and 1990s, typically only one or two sensors were active at any given time. However, since 2005, at least three NOAA or MetOp satellites equipped with AVHRR sensors have been operational simultaneously, with even four satellites providing data for albedo retrieval in 2007, 2009, 2010, and 2013. Furthermore, the CLARA-A3 dataset includes cloudiness probabilities for each imaged AVHRR pixel, obtained through Bayesian classification [29], which will be used for masking out the cloudy top-of-atmosphere (TOA) reflectance during the atmospheric correction process. The snow flags were also included to identify snow-free pixels as part of the cloud processing [17].

2.2. Ancillary Data

In addition to AVHRR GAC TOA reflectance, the inputs for generating the albedo product primarily include data from the atmospheric correction process, specifically total column water vapor, ozone content, surface pressure, and aerosol optical depth (AOD) at a wavelength of 550 nm. A detailed description of this process is provided in Section 3.2, atmospheric correction.

2.2.1. ERA5 Meteorological Data

ERA5 is the fifth-generation atmospheric reanalysis product of the European Centre for Medium-Range Weather Forecasts (ECMWF), providing global climate data from 1940 to the present, with a temporal resolution of 1 h and a spatial resolution of 0.25° [30]. It is an updated version of the ERA-interim dataset, benefitting from the continuous development of numerical meteorological models and improved methods [31]. ERA5 has been widely acknowledged for its reliable performance across various parameters and is frequently used for atmospheric correction in various sensors [13,32,33]. In this study, ERA5 hourly data, including total column water vapor, surface pressure, and total column ozone, were used for atmospheric correction. These inputs are consistent with those used in the CLARA-A3 albedo product [13]. The needed ERA5 meteorological data were bilinearly interpolated to match the 0.05° × 0.05° spatial resolution of the AVHRR GAC data.

2.2.2. MERRA-2 AOD Data

The Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2), developed by National Aeronautics and Space Administration (NASA), provides a consistent and continuous record of weather and climate variables from the 1980s to the present using the Goddard Earth Observing System (GEOS) model with advanced data assimilation methods [34]. The original MERRA-2 dataset, M2I3NXGAS, with 3 h temporal resolution and a spatial resolution of 0.625° × 0.5° [35], was aggregated into a daily product and interpolated to a 0.05° × 0.05° spatial resolution to align with the AVHRR GAC data. MERRA-2 AOD data have shown good agreement with worldwide AERONET measurements [36,37,38] and are used as the unique AOD source for atmospheric correction in this study. Noted is that the AOD data for 1979 was substituted with MERRA-2 AOD data from the corresponding day in 1980.

2.3. Global BRDF/Albedo Satellite Products

In order to comprehensively evaluate our generated albedo product namely GAC43, two additional satellite-based albedo products, CLARA-A3 and C3S, are used for comparisons of the entire study period. The widely used MCD43 BRDF/albedo product will be used as the ideal reference for intercomparison after the year 2000.

2.3.1. CLARA-A3 Albedo Data

The CLARA-A3 albedo product, developed by CM SAF/EUMETSAT, spans the period from 1979 to the present. It is based on intercalibrated AVHRR data and provides black-sky albedo (BSA), white-sky albedo (WSA), and blue-sky albedo at a spatial resolution of 25 × 25 km, with temporal resolutions of pentads and months. BSA is defined as the albedo in the absence of diffuse component and varies with the SZA. WSA is defined as the albedo in the absence of direct component when the diffuse component is isotropic [39]. Blue-sky albedo represents the actual albedo and is calculated as a weighted average of BSA and WSA, where the weight corresponds to the fraction of diffuse skylight determined by the instantaneous SZA and AOD [11,40]. Blue-sky albedo is primarily used to evaluate with in situ measurements, whereas BSA is utilized for spatial comparisons of satellite-derived albedo products. Atmospheric correction is performed using the Simplified Method for Atmospheric Correction (SMAC) algorithm, which relies on ERA5 meteorological data for total column water, surface pressure, and total column ozone, but excludes AOD [13]. The AOD data are instead derived from TOMS and Ozone Monitoring Instrument (OMI) Aerosol Index (AI) observations [41]. For BRDF correction, atmospherically corrected surface reflectance, excluding snow-cover areas, are used with the Roujean kernel [19], which accounts for volumetric and geometric scattering. To account for seasonal BRDF variations driven by vegetation phenology, the TOA normalized difference vegetation index (NDVI) is applied. This is derived separately for specific land cover types, including forest, barren land, grassland, and cropland, based on their distinct physical characteristics [42]. It means the kernel coefficients are empirically determined as functions of TOA NDVI at various land cover types. However, inaccuracies may arise if land use classification does not align with the actual physical characteristics of the surface. Importantly, the CLARA-A3 albedo estimates represent the mean SZA at corresponding compositing period, including pentad and monthly, rather than fixed geometric conditions. Therefore, it should be normalized to fixed geometric conditions, specifically SZA at local solar time and nadir VZA for intercomparison purposes [43]. Lastly, to maintain spatial consistency, the original 25 km spatial resolution of CLARA-A3 is bilinearly interpolated to match the 5 km spatial resolution of GAC43.

2.3.2. C3S Data

The C3S of the European Union is hosted by VITO, which provides global BSA at local solar noon and WSA with a 10 day temporal resolution and three different spatial resolutions. The C3S albedo product is derived from satellite observations, including NOAA AVHRR (1981–2005, ~4 km), SPOT-VGT (1998–2014, ~1 km [14]), PROBA-V (2014–2020, ~1 km [15,44]), and the newest member, Sentinel-3 OLCI/SLSTR (from 2018 onwards, ~300 m [16]). To ensure consistency across the entire period, spectral harmonization is applied to account for differences in spectral response functions [15], except for Sentinel-3. However, some inconsistencies remain in the atmospheric correction and BRDF inversion processes when compared to other albedo products. For instance, both AVHRR-based and SPOT VGT-based albedo products in C3S apply MERRA-2 data for total column ozone, total water vapor column, surface pressure, and AOD, while the PROBA-V-based albedo product uses monthly AOD from the MODIS product, with the other meteorological inputs remaining the same. The Sentinel-3 product employs daily Copernicus Atmosphere Monitoring Service (CAMS) near-real-time (NRT) products due to the NRT requirement. Regarding the BRDF model, all four albedo products within the C3S framework apply the RossThick kernel [19] for volumetric scattering and the LiSparseR kernel [45,46] for geometrical scattering, with the exception of PORBA-V, which uses the Roujean kernel for geometric scattering. Additionally, the compositing periods differ: 30 days for PROBA-V and 20 days for AVHRR, SPOT, and Sentinel-3. It is important to note that the Sentinel-3 albedo is considered pre-operational and is therefore excluded from intercomparison within the C3S framework. Furthermore, C3S albedo data with an associated uncertainty (ERR) greater than 0.2 or a mean age (AGE, in number of days) greater than 20 are discarded, as they are considered lower-quality data [15]. Lastly, C3S is interpolated to the spatial resolution of GAC43 using the same method as CLARA-A3.

2.3.3. MODIS BRDF/Albedo Data

The latest MCD43 collection V6 products are generated daily using a 16 day moving window of observations with predefined weights from MODIS data onboard the Terra and Aqua satellites. These products provide BSA at local solar noon and nadir VZA, and WSA, spanning from 2000 to the present. As in the previous collection V5, the second simulation of the satellite signal in the solar spectrum (6S) model [47] is employed for atmospheric correction. The primary inputs for the model include MCD04 AOD, MCD05 water vapor, MCD07 ozone, and MCD35 cloud mask, all derived from MODIS products, while surface pressure is obtained from the Data Assimilation Office. The widely used RTLSR model is applied for BRDF correction in both snow-free and snow-covered conditions, depending on the dominant surface condition during the retrieval period. Additionally, a major improvement in collection V6 is the backup inversion method, which uses pixel-based updates from the latest full inversion, rather than relying on a land cover-based database as in collection V5. In this study, the aggregated MCD43C3 product with a 0.05° spatial resolution is used as the global albedo reference product after 2000, and the original 500 m spatial resolution product (MCD43A3) is used to distinguish between full and backup inversions. As with the C3S albedo product, only high-quality data indicated by quality control (QC) flags from the MCD43C3 product are used for site-scale validation and spatial intercomparison.

2.4. In Situ Measurement

The sites selected for validating satellite-based albedo estimates were primarily chosen from the SURFRAD [48] and FLUXNET [49] networks, spanning between 1996 and 2020, as illustrated in Figure 2. The SURFRAD sites are all located within the United States, while the FLUXNET sites are mainly distributed across North America, Europe, and Australia, representing diverse land cover types. These sites encompass a range of land surface categories, including savannas (BSV), croplands (CRO), deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), evergreen needleleaf forest (ENF), grasslands (GRA), mixed forest (MF), open shrubland (OSH), permanent wetlands (WET), and woody savannas (WSA), as defined by the International Geosphere-Biosphere Programme (IGBP) in the MODIS collection 6 land cover data (MCD12C1) [50,51]. Both SURFRAD and FLUXNET networks are widely used for surface albedo validation due to their established protocols, long-term operations, and robust infrastructure [52,53,54,55]. FLUXNET provides only downward and upward shortwave radiation at a half-hour resolution, from which nominal blue-sky albedo is calculated. In contrast, SURFRAD provides more detailed measurements, including direct and diffuse downward shortwave radiation along with upward shortwave radiation, at higher resolutions (1 or 3 min). This allows for more accurate downward radiation measurements compared to that from FLUXNET. Daily AOD values at 1 km spatial resolution were extracted from the MCD19A2 [56] product to calculate diffuse skylight using the 6S radiative transfer model when diffuse radiation data were unavailable. If MCD19A2 AOD values were missing, a constant AOD value of 0.1 was used. Additionally, SURFRAD measurements utilized total sky imagers to mask cloud-contaminated values, while FLUXNET measurements used the ratio of measured to potential downward shortwave radiation to classify observations as cloudy when the ratio was less than 0.6. For snow-covered pixels, daily global snow cover data from the European Space Agency (ESA) Climate Change Initiative Snow (Snow_cci) project was used to mask snow-covered pixels. Since this study focuses on albedo at local solar noon, the measured local noon blue-sky albedo was calculated by averaging observations between 11:00 and 13:00 local time, reducing the influence of potential outliers. Ultimately, cloud-free and snow-free ground albedo measurements at local solar noon were generated for validations. For temporal matching between in situ measurements and satellite-based albedo products, all valid measurements at each site within the individual satellite compositing period were averaged, except for the daily MCD43 albedo product. Theoretically, the albedo measured from a tower covers a circular footprint that depends on the tower height in combination with the beam width of the measurement device that should ideally correspond to the pixel size of satellite estimates. However, satellite footprints are often much larger than tower footprints. Consequently, the representativeness of the surrounding landscape at each site is evaluated using geostatistical indices, such as NDVI and digital elevation model (DEM) data, based on the semi-variogram model [57] as shown in

Table 1. The representativeness indices at selected homogeneous sites, where the ST_NDVImax, ST_NDVImin, and ST_DEM represent the overall measures of the spatial representativeness using three indicators, including maximum NDVI, minimum NDVI, and digital elevation model (DEM), respectively.

Site	Lat (°N)	Lon (°S)	Network	IGBP	STNDVImax	STNDVImin	STDEM
AU-Tum	−35.66	148.15	FLUXNET	EBF	2.39	5.53	1.72
AU-Wac	−37.43	145.19	FLUXNET	EBF	1.30	1.17	1.69
CA-NS6	55.92	−98.96	FLUXNET	OSH	2.11	2.62	2.64
CA-Oas	53.63	−106.20	FLUXNET	DBF	13.60	1.93	4.35
CA-SF3	54.09	−106.01	FLUXNET	ENF	1.65	2.18	8.33
DE-Geb	51.10	10.91	FLUXNET	CRO	15.40	\	2.60
DE-Hai	51.08	10.45	FLUXNET	DBF	5.65	5.87	1.05
DE-Kli	50.89	13.52	FLUXNET	CRO	1.47	1.62	0.0023
FI-Hyy	61.85	24.30	FLUXNET	ENF	\	2.87	1.85
FR-Pue	43.74	3.60	FLUXNET	EBF	\	\	2.06
IT-Col	41.85	13.59	FLUXNET	DBF	\	7.80	2.53
RU-Che	68.61	161.34	FLUXNET	MF	0.53	\	6.11
SF_DRA	36.63	−116.02	SURFRAD	BSV	0.62	1.31	1.77
SF_FPK	48.31	−105.10	SURFRAD	GRA	10.70	2.79	4.05
SF_GCM	34.25	−89.87	SURFRAD	GRA	4.09	\	0.00043
SF_PSU	40.72	−77.93	SURFRAD	CRO	2.91	2.28	7.89
SF_SXF	43.73	−96.62	SURFRAD	OSH	5.55	\	2.13
SF_TBL	40.13	−105.24	SURFRAD	GRA	\	7.78 × 10−6	2.03
US-ARM	36.61	−97.49	FLUXNET	CRO	2.08	1.88	2.28
US-Ivo	68.49	−155.75	FLUXNET	WET	\	1.03	2.42
US-MMS	39.32	−86.41	FLUXNET	DBF	3.41	1.37	2.40
US-Me2	44.45	−121.56	FLUXNET	ENF	\	2.96	2.24
US-Ne2	41.16	−96.47	FLUXNET	CRO	1.65	1.47	1.76
US-SRM	31.82	−110.87	FLUXNET	WSA	3.27	0.16	3.06
US-WCr	45.81	−90.08	FLUXNET	DBF	1.51	3.51	0.00027

. We consider that sites are heterogeneous or not spatially representative when no indicator is larger than 1.2 (see Table 1), and sites with higher indicator values are likely to be more representative of their surroundings than those with lower values. More detailed information about the representativeness of the measurement sites is provided by https://savs.eumetsat.int/ (accessed on 28 November 2024) [58]. It is worthy to note that satellite pixels are not necessarily internally homogeneous. Additionally, 19 pseudo-invariant calibration sites (PICS) are used to evaluate the stability of satellite-based albedo products [52], as these sites are generally recognized for their high stability.

3. Algorithm Description

3.1. Method Overview

AVHRR channels 1 and 2, with wavelengths of 0.58–0.68 μm (CH1) and 0.725–1.10 μm (CH2), respectively, serve as the radiance source for GAC43 albedo generation. The workflow diagram of this study is illustrated in Figure 3. It is divided into two main sections: the product generation stage outlined by light grey frames, and the product evaluation stage represented by a dark grey frame. The data processing begins by excluding cloud-contaminated GAC pixels from the EUMETSAT AVHRR fundamental data record (FDR). Next, GAC TOA reflectance data collected under favorable geometric conditions are subjected to atmospheric correction to produce top-of-canopy (TOC) or atmospherically corrected surface reflectance. Subsequently, a BRDF model is used to normalize surface anisotropy for TOC reflectance. The resulting BRDF parameters are then integrated over the illumination and viewing hemisphere to derive the spectral albedo. Finally, a narrow-to-broadband conversion was applied to transform the spectral albedo into broadband albedo. For the evaluation stage, in situ measurements and selected satellite-based albedo products are used for comparison at the site scale and spatial distributions, respectively.

3.2. Atmospheric Correction

For the atmospheric correction of the AVHRR GAC data, the SMAC [59] updated to the 6S [47] radiative transfer method was employed, where all the pertinent radiative quantities are parameterized as a function of auxiliary data as mentioned in Section 2.2 Ancillary data. The SMAC approximation to the 6S method is reported to achieve accuracy within 1% for most scenarios under the assumption of a Lambertian (i.e., isotropic) surface. However, the SMAC sacrifices precision for processing speed through utilizing simplified and parameterized equations for radiative transfer. Consequently, observations with unfavorable illumination or observation geometry (i.e., SZA > 70° and VZA > 60°), cloud probability greater than 20%, or AOD values exceeding 1.0 were discarded [60]. The inputs to SMAC include TOA reflectance, geometry information, pixel geolocation, and ancillary meteorological data, while it outputs TOC directional reflectance. It is noteworthy that AOD is one of the most significant factors that affects atmospheric correction. Therefore, the SMAC algorithm provides an option to compute corrections especially for desert areas, characterized by dry conditions with a high AOD content. The land cover data used to identify desert and non-desert regions, essential for applying different SMAC coefficients, are derived from the yearly MCD12C1 land cover data. When MCD12C1 data are unavailable, the most recent year’s MCD12C1 data are used instead. Mathematically, the SMAC algorithm originally solves the surface reflectance following an atmospheric scattering equation based on an inverse solution, which could be expressed as follows:

$$\rho \left({\theta }_s,{\theta }_v,\Delta \phi \right)={\displaystyle \frac{{\rho }_{TOA}-t_g{\rho }_a\left({\theta }_s,{\theta }_v,\Delta \phi \right)}{t_{g}T\left({\theta }_s\right)T\left({\theta }_v\right)+S{\rho }_{TOA}-t_g{\rho }_a\left({\theta }_s,{\theta }_v,\Delta \phi \right)S}}$$

(1)

where ${\rho }_{TOA}$, $t_g$, ${\rho }_a$, $T\left({\theta }_s\right)$, $T\left({\theta }_v\right)$ and S represent TOA reflectance, total gaseous transmission, atmospheric reflectance, downward scattering transmission, upward transmission, and spherical albedo of the atmosphere as viewed from the ground, respectively, which all are provided by SMAC coefficients modeled by 6S previously. For more details about the SMAC algorithm, refer to the algorithm theoretical basis document (ATBD) of the CLARA-A3 albedo.

3.3. BRDF Inversion

GAC43 BRDF/albedo is computed using of the MODIS BRDF/albedo algorithm corresponding to the MCD43 collections V5. The BRDF retrieval algorithm operates on an 8 day overlapping procedure and two daily rolling modes, with the center of the moving 16 day input window aligned with the date of each retrieval. This 16 day compositing period strikes an appropriate balance between the availability of sufficient angular samples and the temporal stability of the land surface, although this assumption of stability may fail during periods of significant phenological change, such as vegetation green-up, senescence, or harvesting. The BRDF model used in this study includes the RossThick component to account for volumetric scattering, and the LiSparseR [45,61] component to represent geometric scattering. The model is described by the following expression:

$$R=f_0+f_1{k}_1({\theta }_s,{\theta }_v,\phi )+f_2{k}_2({\theta }_s,{\theta }_v,\phi )$$

(2)

where the R, subscript 0, subscript 1, and subscript 2 denote observed reflectance, nadir reflectance, volume scattering term, and geometrical scattering term, respectively. Additionally, the terms of k and f represent the model kernel functions and weighting parameters, respectively. Notably, the isotropic kernel function is always set as a constant value of 1.0. For more details and complete mathematical expressions regarding the three kernels models, we refer to [39]. The determination of the kernel coefficients is performed using a linear least-square fit, which carries more statistical significance than using TOA NDVI to determine the kernel coefficients in the BRDF inversion of the CLARA-A3 albedo. Once both kernel coefficients and functions are determined, surface observations under any illumination and viewing conditions can be simulated.

Additionally, the GAC43 albedo provides a quality index to identify whether it is retrieved through a high-quality full inversion or low-quality backup inversion. Three criteria must be met for the full inversion of GAC43; first, at least seven clear-sky observations must be available within 16 days; second, the root mean square errors (RMSEs) of the kernel model fit must be less than 0.08 for both spectral channels; and third, the weight of determinations (WODs) needs to be less than 1.65 for nadir BRDF-adjusted reflectance at 45° SZA and 2.5 for WOD–WSA. It is worth mentioning that the above uncertainty indicators are constructed during the MCD43 albedo inversion process [62]. To simplify the uncertainty comparison between GAC43 and MCD43, this study applied consistent metrics with MCD43. If any of these conditions are not met, a backup inversion is performed. For backup inversion, the mean BRDF climatology from high-quality full inversions on the corresponding date between 2005 and 2020 is used as priori information to improve performance. This approach assumes that it can provide more high-quality albedo inversions, thereby more stable BRDF climatology, as at least three AVHRR constellations are available during the given period. Furthermore, if more than three available observations are available within 16 days for backup inversion, the mean BRDF climatology is used as the initial value of kernel coefficients and serves as the inputs for subsequent linear least-square fitting. If no climatology is available, the initial value of kernel coefficients is set to [0.1, 0, 0] for the three kernels, followed by the fitting process. Otherwise, if fewer than three observations can be collected within the compositing period, the GAC43 pixel would be filled with the BRDF climatology. The backup inversion values should be used with caution, as they are usually considered lower quality inversions and are flagged accordingly in the QC field.

3.4. Albedo Computation

After determining the BRDF kernel coefficients and functions, the spectral albedo is obtained by integrating over all viewing and illumination directions, resulting in the BSA at desired SZA and WSA for the entire hemisphere. Since the instruments only observe radiances in discrete spectral bands, climatological applications require deriving the total shortwave broadband albedo. This is achieved by extending the estimated spectral albedo to the full broadband spectral range through a process known as narrow-to-broadband conversion, as described in [63]:

$$\alpha =-0.3376 \times {\alpha }_{RED}^2-0.2707 \times {\alpha }_{NIR}^2+0.7074 \times {\alpha }_{RED} \times {\alpha }_{NIR}+0.2915\times {\alpha }_{RED}+0.5256 \times {\alpha }_{NIR}+0.0035$$

(3)

However, only the blue-sky albedo is directly comparable to field observations. To address this, we converted the BSA and WSA to the blue-sky albedo using the function below:

$${\alpha }_{blue-sky}=f_{dir} \times {\alpha }_{BSA}+(1-f_{dir}) \times {\alpha }_{WSA}$$

(4)

here ${\alpha }_{BSA}$ and ${\alpha }_{WSA}$ are computed by Equation (3) after determining spectral albedo, while $f_{dir}$ represents the fraction of direct skylight, which can be measured at sites. However, some sites used in this study do not provide direct radiation measurements. In such cases, AOD values from the daily MCD19A2 product at 1 km spatial resolution are used to calculate the $f_{dir}$ using the 6S radiative transfer model, which utilizes a previously computed look-up-table (LUT) [64]. If AOD data are unavailable, a constant AOD value of 0.1 is used for the calculation.

3.5. Evaluation Metrics

The final albedo estimates obtained from the GAC43 albedo product were evaluated in terms of direct validations against homogeneous in situ measurements, spatial intercomparisons with three other satellite-based products including CLARA-A3, C3S, and MCD43C3, and assessments of inversion quality consisting of temporal stability and the percentage of full inversions. Three commonly used statistical metrics as below would be applied to quantify accuracy such as the coefficient of determination ($R^2$), RMSE, and bias. Additionally, linear regression analysis was conducted to evaluate trend dynamics at PICS sites.

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(e_i-o_i)}^2}{{\sum }_{i=1}^n{(o_i-\overline{o})}^2}}$$

(5)

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(o_i-e_i)}^2}$$

(6)

$$\mathrm{Bias}={\displaystyle \frac{1}{n}}\sum _{i=1}^n(e_i-o_i)$$

(7)

where $e_i$ is the estimated albedo, $o_i$ is the corresponding ground measurement, and $\overline{o}$ is the average value of ground measurements.

4. Results

4.1. Validation at Site Scale

The performance of full inversion and full plus backup inversions based on newly developed GAC43 albedo products was first evaluated across various land cover types, as shown in Figure 4. It is worth mentioning that approximately 85% of the matching samples come from full inversions; therefore, the performance of the combined full and backup inversions is expected to be primarily determined by the retrievals from the full inversions. The results indicate comparable performance across most land surface types, except for the ENF type, where the full plus backup inversion performs slightly worse than the full inversion. Conversely, for the OSH type, the full inversion performs worse than full plus backup inversion. This discrepancy is probably due to an unrepresentative statistical bias caused by the limited number of OSH samples, fewer than 100. Although the backup inversions demonstrate reasonable accuracy, we recommend using GAC43 albedo values derived from it with caution when conducting related research, as suggested in previous studies [11,23,65]. Additionally, more detailed performance of the GAC43 albedo with full inversions across diverse land surface types is presented in Figure 5. The evaluation of GAC43 retrievals against in situ measurements is limited to sites with the highest homogeneity to minimize the effect of scale mismatch. Figure 5 shows that GAC43 tends to overestimate albedo in forest areas at GAC (5 km) spatial resolution, particularly in EBF, followed by DBF, while ENF is the least affected, with biases of 0.028, 0.024, and 0.014, respectively. This is likely due to non-negligible surface heterogeneity within an AVHRR GAC footprint. It is reasonable to expect that forest cover is not internally homogeneous across the AVHRR GAC pixel, and the presence of other low-density vegetation types is likely to increase the resulting albedo as shown in Figure 6a. Furthermore, some land surfaces, such as DNF, experience more significant seasonal heterogeneity due to phenological changes. For BSV areas, GAC43 tends to underestimate albedo as shown in Figure 5a, possibly due to the influence of other scattered vegetation types, roads, or roofs captured in satellite observations as shown in Figure 6b. However, the GAC43 albedo does not perform very well in CRO and GRA, as shown in Figure 5b,f, likely due to human activities or complex pixel patterns, as shown in Figure 6c,d, which are often associated with higher heterogeneity. Overall, these biases are probably due to the fine-scale spatial variability of the plant cover and the insufficient representativeness of tower footprint in the AVHRR GAC pixel.

To provide a more comprehensive evaluation of the GAC43 albedo product, Figure 7 shows its accuracy with three other satellite-based products: CLARA-A3, C3S, and MCD43C3. Figure 7a–d presents the overall accuracy of the four albedo products using all available sample data matched with site measurements, while Figure 7e–h shows their cross-accuracy using the same samples. It is worthy to note that these satellite-based albedo products are generated based on various composing periods, which may affect accuracy, particularly during rapid surface albedo changes, such as snow or fire events. To minimize these effects, the comparisons were made using the closed centers of their respective temporal composites. It is evident that the widely used MCD43C3 product provides the most accurate albedo estimates, with an RMSE of 0.023 and a minimal bias as shown in Figure 7h. GAC43 and C3S follow, both exhibiting similar accuracy with an RMSE of approximately 0.03, although both show a tendency to overestimate albedo. CLARA-A3 is the least accurate among the four albedo products, with the highest RMSE (~0.041) and a notable positive bias of 0.02 as shown in Figure 7f. Figure 8 presents the accuracy statistics for the four albedo products across classical biome types, using RMSE (Figure 8a) and bias (Figure 8b) as evaluation metrics. For forested areas, GAC43 and C3S show similar levels of uncertainty (RMSE) and a positive bias, with values around 0.029 and 0.023, respectively. In contrast, CLARA-A3 significantly overestimates the albedo in forests, showing a positive bias of approximately 0.035 and the highest estimation uncertainty, around 0.045. Meanwhile, MCD43C3 demonstrates minimal bias and achieves the lowest RMSE of 0.015 for forest surfaces. On grass/shrub surfaces, GAC43 and MCD43C3 exhibit similar uncertainties and positive biases, approximately 0.024 and 0.010, respectively, while CLARA-A3 and C3S have higher uncertainties (0.033–0.038) and noticeable positive biases (0.022–0.026). As expected, the uncertainties of all four albedo products are relatively high on crops surfaces, ranging from 0.034 to 0.038, with deviations generally within ±0.016. Finally, desert surfaces show the lowest estimation uncertainties and biases for all four products, although GAC43 displays a slight underestimation.

In addition to direct validation with site measurements, Figure 9 illustrates the temporal performance of four albedo products across five typical vegetated land surface types at different observation times. It is important to note that this study does not specially evaluate albedo products during snow seasons. Instead, the temporal validation results highlight the products’ capability to capture annual albedo dynamics, including seasonal snow and snowmelt periods, as shown in Figure 9, despite some biases. For instance, all four albedo products tend to overestimate albedo values for the EBF type (Figure 9a), with CLARA-A3 showing the largest overestimation, followed by GAC43 and C3S, while MCD43C3 is closest to the in situ measurements, consistent with the results in Figure 8. A similar pattern appears for the ENF type (Figure 9b). For more complex forest types, such as DBF (Figure 9c), GAC43, CLARA-A3, and MCD43C3 display two distinct biases: overestimations during the snowy seasons and underestimations during non-snowy seasons. However, GAC43 closely aligns with in situ measurements from August to October when there are more full inversions. In contrast, C3S aligns more closely with measurements for DBF from May to July but overestimates from August to October. For the GRA surface type (Figure 9d), all four products perform well with minimal bias. However, for CRO (Figure 9e), all products significantly overestimate albedo during the growing season, while estimates are closer to measured values in the non-growing season. Lastly, discrepancies can arise when snow appears, sometimes resulting in lower-than-expected albedo values and increasing uncertainty, particularly when snow only partially covers pixels. In summary, while all four satellite products effectively capture surface albedo changes, they tend to overestimate albedo over forested surfaces, likely due to mismatches between site representativeness and satellite footprint and the RSLSR model’s limitations under snow cover conditions. However, these limitations are beyond the scope of this study.

4.2. Spatial Performance

In addition to the direct validations between satellite-derived albedo and in situ measurements, spatial distribution differences among albedo products need to be evaluated. The monthly BSA for July 2013 is used as an example because there is less snow-covered surface, illustrating the spatial distributions of GAC43 and the corresponding differences with CLARA-A3, C3S, and MCD43C3 at same month, as shown in Figure 10a–d. Appendix A Figure A1, Figure A2, Figure A3 and Figure A4 present additional comparisons of monthly albedo for July in the years 2004, 2008, 2012, and 2016. The GAC43 albedo demonstrates satisfactory global patterns, with lower albedo over forested areas, higher albedo over deserts, and intermediate values elsewhere as shown in Figure 10a. Furthermore, comparisons between GAC43 and the other three products show that GAC43 and C3S exhibit the closest overall agreement, with approximately 81% of absolute differences within 0.02 (Figure 11b), followed by CLARA-A3 at 67% (Figure 11a). In contrast, only 38% of absolute differences between GAC43 and MCD43C3 fall within 0.02 (Figure 11c), with GAC43 showing higher values than MCD43C3 in nearly 90% of the regions globally. Spatially, the higher GAC43 values compared to CLARA-A3 are concentrated in the northern African desert, southern South America, and eastern coastal areas of China, while lower values are primarily found in southwestern China and eastern Russian. These discrepancies may be attributed to differences in AOD data sources, particularly in areas with high AOD. Compared to C3S in Figure 10c, GAC43 shows slightly higher albedo values primarily in the Northern Hemisphere, including eastern and southern China, northern India, and eastern North America. In contrast to the MCD43C3 albedo products, GAC43 shows higher values across most regions, except for southern Africa, central Australia, and south-central South America. It is worth noting that in regions with dense forest cover, such as the Amazon, central Africa, and the European plains, the albedo values from GAC43, CLARA-A3, and C3S are relatively consistent and tend to overestimate compared to MCD43C3.

In addition to the spatial comparison of satellite albedo products, Figure 12 further presents the GAC43 monthly BSA alongside the corresponding CLARA-A3, C3S, and MCD43C3 products, evaluated on a grid-cell basis for July 2013. These scattering results closely align with the above results as shown in Figure 10 and Figure 11, demonstrating that the albedo estimates from GAC43 and C3S are the closest with negligible bias and the smallest RMSE of 0.017. Although CLARA-A3 also shows minimal overall bias compared to GAC43, their distributions are more dispersed, leading to a higher RMSE of 0.027. As expected, GAC43 BSA has the largest positive bias of 0.025 and the highest RMSE of 0.032 compared to MCD43C3 BSA. Additionally, Figure 13 further illustrates the BSA distributions among the four satellite albedo products across various land cover types. It reveals that GAC43 and C3S show very similar albedo distributions across almost all land cover types, regardless of the degree of vegetation cover. While GAC43 exhibits higher albedo values than MCD43C3 across nearly all land cover conditions, including diverse forest types, CRA, CRO, SAV, and WSA, this is consistent with Figure 10. When comparing GAC43 and CLARA-A3, the albedo differences are relatively small under land cover types such as CRO, GRA, SAV, and WSA. However, for almost all other forest cover types besides MF, CLARA-A3 tends to retrieve higher albedo values than GAC43.

The similarity between the albedo estimates derived from GAC43 and C3S is primarily attributed to both products utilizing the SMAC algorithm for atmospheric correction and the RTLSR model for BRDF correction, despite differences in the data used for atmospheric correction and the spectral response functions of sensors. GAC43 and CLARA-A3 are nearly identical in their atmospheric correction stage, with the only difference being the selection of AOD data sources. However, the kernel models used by the two products differ slightly, particularly in their geometric scattering functions, and their kernel coefficients are entirely distinct. Additionally, CLARA-A3 does not provide BSA at local solar time, which may introduce bias for end users, although normalizations were applied in this study to mitigate this issue. In contrast, the MCD43 BRDF/albedo product uses the more accurate 6S model for atmospheric correction, relying almost entirely on internally sourced data. While the applied BRDF model is the same, the MODIS sensor has a more detailed band setting, contributing to more precise broadband albedo retrieval.

4.3. GAC43 BRDF/Albedo Inversion Quality

After comprehensive evaluation against in situ measurements and spatial intercomparison with three other satellite-based albedo products (CLARA-A3, C3S and MCD43), the inversion quality of the GAC43 BRDF/albedo is further assessed based on the temporal stability at PICS sites and the frequency of full or backup inversions globally.

Given that this study focuses on albedo dynamics over nearly 40 years, it is crucial for long-term albedo products to exhibit stability. PICS are typically used to assess stability, as they are expected to experience minimal temporal variations. Figure 14 displays examples of monthly average albedo profiles during 1982–2020 at three randomly selected desert calibration sites. Additionally, Figure 15 presents the decadal trends in monthly BSA across four periods (1982–1990, 1991–2000, 2001–2010, and 2011–2020) for four satellite-based albedo products. The GAC43 product (shown in green scatters) demonstrates a stable temporal performance throughout almost the entire period at each site, as shown in Figure 14, with most slopes remaining within ±0.002 per decade and showing no significant fluctuations in Figure 15a–d. In contrast, the C3S product exhibits notable anomalies during transition periods, particularly between NOAA-11 and NOAA-14 in 1995 and during the NOAA-7 era (1982–1984), likely due to imperfect sensor calibration as shown in Figure 14. However, its stability improves significantly with the use of SPOT from 1998, with most slopes remaining within ±0.00125 per decade after the 2000s as shown in Figure 15c,d. The CLARA-A3 albedo temporal trajectories show significant temporal noise, especially before 2002 as shown in Figure 14, which is largely attributed to the limited availability of observations for anisotropy correction. Notably, both GAC43 and CLARA-A3 exhibit a slight decreasing trend during 2011–2020, particularly for CLARA-A3 after 2013, as indicated by the lower monthly albedo values in Figure 14. Furthermore, the MCD43C3 product shows strong stability after 2000, with a slight increasing trend within ±0.000625 per decade as shown in Figure 15c,d. Overall, the GAC43 albedo product demonstrates a more consistent stability than CLARA-A3 and C3S over the entire study period, despite a minor tendency to decrease since the 2010s. Therefore, GAC43 is expected to be a highly promising data source for long-term albedo research with overall stability satisfactory.

Another important consideration is the frequency of BRDF/albedo full inversions, because full inversions without prior information are usually considered to be of better performance than backup inversions. The spatial distributions of global full inversion percentages between GAC43 (top) and MCD43A3 (bottom) for the years 2004, 2008, 2012, and 2016 were shown in Figure 16 as examples. Even though the MCD43 BRDF/albedo product provides daily product, the median day of that within the 16 day synthesis period at an 8 day temporal resolution was selected for an easy comparison with the GAC43 inversion quality. The result indicates that GAC43 exhibits a significantly higher frequency of full inversions particularly below 50°N, although a lower frequency of full inversion at latitudes above 50°N, compared to MCD43A3. This discrepancy arises because GAC data selects the pixel with the smallest VZA when multiple observations occur on the same day, which increases the likelihood of high-quality BRDF inversion failure at higher latitudes due to fewer pieces of diverse geometry information. However, GAC43 provides significantly more full inversions than MCD43 even in cloudy areas, such as the Amazon, central Africa, southeast Asia, and southeast China, potentially but not necessarily supporting more reliable analyses. In general, the BRDF inversion quality is primarily influenced by cloud cover and geometry observations during the synthesis period. For example, both albedo products show a high percentage of full inversions across north Africa, central Australia, western United States, and some other areas. Overall, it successfully demonstrates that more available satellite observations could effectively improve BRDF/albedo inversion quality.

The percentage of full inversions was used further to assess the inversion quality of GAC43 data throughout the study period from 1979 to 2020, as illustrated in Figure A5. The result indicates that the percentage of full inversions was significantly lower in areas with frequent rain events before 2000s due to fewer available satellite observations. In contrast, dry regions, such as north Africa, south Africa, central Australia, and the Middle East, maintained relatively high percentages of full inversion. Although the GAC43 albedo provided backup inversions to address this limitation, the scarcity of observations likely reduced its effectiveness during this period. However, starting in the 2000s, particularly after 2002, GAC43 produced a markedly higher number of full inversions globally, facilitated by increased satellite observations. For example, almost all land areas below 50° had more than 70% full inversions in 2013, including parts of the Amazon and central Africa, though southeast Asia remained slightly lower, at around 50%. Figure 17 further breaks down the percentage of full inversions over non-ice/snow-covered surfaces by continents. It shows that before the 2000s, only Oceania and Africa had a full inversion percentage of around 40%, due to the presence of large deserts and arid environments. In contrast, the percentages of full inversions on other continents were below 20% during same time period. After 2002, full inversion percentages rose significantly across all continents, with Europe, North America, and Asia stabilizing between 30% and 50% (Europe averaging 30%, North America around 40%, and Asia reaching 50%). Meanwhile, South America, Africa, and Oceania recorded higher percentages, ranging from 70% to 90%. These statistics demonstrate that since the 2000s, GAC43 albedo data have been more stable and of high quality, compared to the larger uncertainties of the 1980s and 1990s when more backup inversions were required. Nevertheless, GAC43 still provides critical albedo data through backup inversions before the 2000s, despite these uncertainties.

5. Discussion

This study introduces a newly developed global land surface albedo dataset namely GAC43 spanning from 1979 to 2020, which includes both complete atmospheric correction and BRDF correction. Through robust sensors calibration and a rigorous inversion process, GAC43 demonstrates satisfactory stability at PICS sites, including periods before the 2000s. In contrast, the CLARA-A3 and C3S datasets display limited stability, especially before 2000, with significant noise or biases observed, respectively. Another key feature of GAC43 is its use of all available AVHRR sensors for anisotropic correction, which enhances the variety of ground-viewing geometries, allowing for more successful high-quality inversions. Comprehensive evaluations against in situ measurements and three other satellite-based albedo products reveal that GAC43 performs well in direct validation, temporal variation, and spatial distribution assessments. This product is a valuable albedo data source for long-term studies related to global climate change, energy balance, and land management policy. However, some limitations and uncertainties remain that should be mentioned.

Firstly, GAC43 relies on only two bands from AVHRR sensors to retrieve broadband surface albedo, which increases systematic inversion deviation due to the coarse spectral settings of the AVHRR sensors. In contrast, the MODIS sensor provides seven shortwave bands to generate the MCD43 albedo product, potentially allowing for more accurate albedo retrievals. Secondly, the AVHRR sensor was not originally designed for quantitative retrievals based on the visible/infrared channels and has inherent limitations [66,67]. For example, because the AVHRR sensor lacks sufficient spectral information to retrieve atmospheric components as effectively as MODIS, it must rely on reanalysis products like ERA5 or MERRA-2 for atmospheric correction. Thirdly, GAC43 employs a simple model as SMAC for atmospheric correction, which is more efficient but less accurate. Although SMAC claims to be within 1% difference of the 6S model under favorable conditions, it is less effective in challenging atmospheric conditions such as high AOD content or observation angles. Additionally, the eight day temporal resolution of GAC43 limits its ability to capture rapid albedo changes in vegetation, especially during growth and senescence periods. The relatively coarse 5 km spatial resolution of GAC data restricts its applicability for regional albedo studies, as it cannot provide finer spatial details. Lastly, due to the limited availability of in situ measurements, comprehensive assessments could not be conducted before the 2000s.

Discrepancies between satellite-based albedo estimates can be contributed to several factors, including differences in spectral response functions, sensor footprints, data sources or models for atmospheric correction, the kernel models used in BRDF correction, and varying compositing periods. For instance, although the same RTLSR model is used for BRDF correction in GAC43, C3S, and MCD43, slight differences remain. C3S and MCD43, for example, apply semi-Gaussian distribution [15] and Laplace distribution [11], respectively, to assign temporal weights to observations within the compositing period, emphasizing the day of interest. In contrast, GAC43 applies uniform weighting over the retrieval period. Additionally, mismatches between site representativeness and satellite footprint significantly affect site-scale validations. The sites measurements only represent a limited surface area, dependent on their tower heights and observation angles, making it challenging to find a homogeneous surface within the 5 km spatial range of the GAC pixel.

6. Conclusions

This study developed a global land surface albedo product, the GAC43 albedo, which provides black-sky albedo (BSA), white-sky albedo (WSA), and potentially blue-sky albedo in the future, using AVHRR GAC data on board NOAA and MetOp satellites spanning from 1979 to 2020. The study employs the less accurate but more efficient SMAC algorithm for atmospheric correction and classical RTLSR kernels for BRDF correction. By leveraging all data from multiple AVHRR sensors, the GAC43 albedo generates a higher number of full inversions than the widely used MCD43 BRDF/albedo products, even in rainy regions like the Amazon, central Africa, and southeast Asia, thus enabling more reliable analyses in tropical rainforest areas. Furthermore, the GAC43 albedo product exhibits a slight tendency within ±0.002 per decade at nearly all PICS, although it shows a slight decreasing trend between 2011 and 2020. In contrast, CLARA-A3 displays significant noises prior to the 2000s due to fewer satellite observations, and C3S shows significant biases at same time period because of imperfect intercalibrations among multiple AVHRR sensors. Therefore, the GAC43 albedo provides a promising long-term and consistent albedo product for future studies.

Comprehensive validation using globally distributed homogeneous sites across four satellite-based albedo products, including GAC43, CLARA-A3, C3S, and MCD43C3, reveals that GAC43 and C3S demonstrate overall similar accuracy, both with an RMSE of approximately 0.03, though GAC43 has a lower positive bias with 0.012. As expected, MCD43C3 performs best against in situ measurements, with an RMSE of ~0.023 and minimal bias. However, CLARA-A3 exhibits the highest RMSE and bias, at 0.042 and 0.02, respectively, consistent with precious study [13]. More accuracy analyses were conducted across various land surface types. Our quantitative assessments show that the GAC43 albedo generally aligns well with in situ measurements, both in terms of direct validation and temporal variations, although slight deviations do occur. For example, GAC43 tends to overestimate in forested areas with positive biases ranging from 0.014 to 0.028, such as EBF, ENF, and DBF. Moreover, while GAC43 shows no significant tendency in GRA and CRO land cover types, it exhibits notable uncertainty for these regions, with RMSEs 0.030 and 0.034, respectively. In terms of spatial performance, GAC43 and C3S exhibit similar behavior across most land cover types. Additionally, CLARA-A3 generally exhibits slightly higher albedo values than GAC43 and C3S under most surface conditions, except for GRA, MF and WSA, where the estimated values are comparable. However, GAC43, CLARA-A3, and C3S consistently show significantly higher albedo values than MCD43C3 across almost all surface conditions. Notably, discrepancies in representativeness between in situ measurements and satellite-based albedo estimates have a substantial impact. Exploring the performance of the longest historical GAC data in albedo inversion is crucial for the potential application of AVHRR data. The proposed GAC43 albedo, which demonstrates consistent stability and satisfactory accuracy, is a valuable resource for future studies. However, caution is advised when using GAC43 albedo estimates before the 2000s due to the limited availability of observations during that period.