Assessment of the EUMETSAT Multi Decadal Land Surface Albedo Data Record from Meteosat Observations

: Surface albedo, deﬁned as the ratio of the surface-reﬂected irradiance to the incident irradiance, is one of the parameters driving the Earth energy budget and it is for this reason an essential variable in climate studies. Instruments on geostationary satellites provide suitable observations allowing long-term monitoring of surface albedo from space. In 2012, EUMETSAT published Release 1 of the Meteosat Surface Albedo (MSA) data record. The main limitation effecting the quality of this release was non-removed clouds by the incorporated cloud screening procedure that caused too high albedo values, in particular for regions with permanent cloud coverage. For the generation of Release 2, the MSA algorithm has been replaced with the Geostationary Surface Albedo (GSA) one, able to process imagery from any geostationary imager. The GSA algorithm exploits a new, improved, cloud mask allowing better cloud screening, and thus ﬁxing the major limitation of Release 1. Furthermore, the data record has an extended temporal and spatial coverage compared to the previous release. Both Black-Sky Albedo (BSA) and White-Sky Albedo (WSA) are estimated, together with their associated uncertainties. A direct comparison between Release 1 and Release 2 clearly shows that the quality of the retrieval improved signiﬁcantly with the new cloud mask. For Release 2 the decadal trend is less than 1% over stable desert sites. The validation against Moderate Resolution Imaging Spectroradiometer (MODIS) and the Southern African Regional Science Initiative (SAFARI) surface albedo shows a good agreement for bright desert sites and a slightly worse agreement for urban and rain forest locations. In conclusion, compared with MSA Release 1, GSA Release 2 provides the users with a signiﬁcantly more longer time range, reliable and robust surface albedo data record.


Introduction
Land surface albedo is one of the Essential Climate Variables (ECVs) defined by the Global Climate Observing System (GCOS) [1]. Surface albedo refers to the ratio of the radiant flux reflected from the surface into the atmosphere to the incident radiant flux [2], and depends on both the anisotropy of the surface and the atmosphere (e.g., [3]). The key elements for albedo retrieval are (i) a robust and reliable Bidirectional Reflectance Distribution Function (BRDF) describing the angular distribution of the radiation reflected at the surface, (ii) the availability of multi-angular surface measurements and (iii) a reliable knowledge of the state of the atmosphere. Finally, a sound assessment of the uncertainty is important for the generation of Climate Data Record (CDR) [4] for climate monitoring. Surface albedo has been retrieved exploiting imagery acquired by instruments on board several platforms, both polar and geostationary satellites [5], such as for instance MODIS [6,7], MISR [8], AVHRR [9,10] and Meteosat Second Generation (MSG) [11]. Observations acquired by geostationary satellites have the advantages of offering both a long-term dataset and an angular sampling of the surface, as well as providing diurnal sampling of key parameters influencing the retrieval, such as cloud cover and aerosol load. At EUMETSAT geostationary land surface albedo is derived using the Geostationary Surface Albedo (GSA) algorithm that exploits imagery acquired from the Meteosat Visible and Infrared Imager (MVIRI) on board Meteosat First Generation (MFG) and from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on board Meteosat Second Generation (MSG). The pixel size for the MVIRI visible channel measurements is 2.5 km at the sub satellite point that is defined as the point where a straight line drawn from a satellite to the center of the Earth intersects the Earth's surface. The SEVEIRI channel spectrally similar to the MVIRI visible channel (Figure 1) is the High Resolution Visible (HRV) channel that has a sub satellite point pixel size of 3 km. The output surface albedo, including associated uncertainties [12], is provided as a 10-day composite in order to minimize the gaps due to cloud coverage. The surface albedo is derived over the full disk, with maximum extensions of 65 • N-65 • S and 65 • W-65 • E around the nominal Sub Satellite Point (SSP) (Figure 2). This algorithm is an improved and extended version of the original Meteosat Surface Albedo (MSA) algorithm developed at the Space Applications Institute of the Joint Research Centre (JRC) of the European Commission in cooperation with EUMETSAT [13]. In GSA the BRDF is described using the Rahman-Pinty-Verstraete (RPV) model [14]. Other algorithms such as the one for MODIS use a semi empirical model based on kernels [15]. Both approaches have pro and contras but are also equally robust and reliable in characterize the BRDF [16]. The first release of the surface albedo data record (called MSA Release 1) was made available in 2012. The validation showed [17] that the MSA Release 1 data record, under various viewing conditions, agreed well with corresponding albedo values derived from other satellites or obtained from in-situ measurements. However, quality issues due to undetected clouds were found. The authors [17] concluded that these quality issues needed to be considered for specific applications and addressed in the context of future improvements. In this paper, the second release of the surface albedo data record (called GSA Release 2) is presented and its quality evaluated.

Retrieval Algorithm
The GSA algorithm retrieves both Black-Sky Albedo (BSA) or Directional Hemispherical Reflectance (DHR) and White-Sky Albedo (WSA) or Bi-Hemispherical Reflectance under perfect isotropic illumination conditions (BHRiso) [3,15]. In GSA, the DHR is calculated assuming a sun zenith angle of 30 • (DHR30). The Global Climate Observing System (GCOS) implementation plan specifies that the DHR is the product required for climate change analysis [1] and for it reason it is the variable analyzed in this paper. The algorithm performs the inversion of a fast Radiative Transfer Model (RTM) ingesting vicarious calibrated visible reflectance [18,19], cloud mask [20,21], Total Column Water Vapour (TCWV) and Total Column Ozone (TCO3) estimates from model reanalysis data produced at the European Centre for Medium-Range Weather Forecasts (ECMWF) [22]. The retrieval algorithm is based on a method proposed by Pinty et al. [13]. The system land-atmosphere is divided in three distinct layers: (i) an anisotropic surface described with the Rahman-Pinty-Verstraete (RPV) model [14], (ii) a lower atmospheric layer where the scattering processes take place, and (iii) a higher atmospheric layer where the gaseous absorption takes place. The RPV models the BRDF as a product of angular functions and a reflectance level ρ 0 [13]. These angular functions depend on the satellite and sun position and three model parameters, i.e., the empirical surface parameter ρ c to characterize the hot-spot, and k, Θ to account for surface anisotropy controlling the shape of the surface bidirectional reflectance factors (BRF) [14]. The model inversion is performed by accumulating cloud free top of atmosphere BRF values over one day at different viewing angles. For each pixel and each day, a minimum of six cloud free daylight observations are needed to retrieve the albedo variables. According to the reciprocity principle, measurements taken at different viewing angles are exchangeable with measurements taken at different sun illumination conditions [23]. Assuming that the geophysical properties controlling the radiance field emerging from a given pixel do not evolve much over a day, the acquisition of radiance data over such a period corresponds to an angular sampling of the same radiance field for various solar geometries. The usage of Meteosat single visible channel data only cannot guarantee a robust and accurate retrieval. In order to constrain the retrieval to the surface-aerosol system only, the atmospheric state is further characterized ingesting estimates of Total Column Ozone (TCO3) and Total Column Water Vapour (TCWV) from the ERA-Interim reanalysis produced at the ECMWF. The approach followed in solving this surface-aerosol scattering problem is an extension of the Multi-angle Imaging SpectroRadiometer (MISR) algorithm for retrieving aerosol optical depth values over dark surfaces [24]. The GSA algorithm also estimates the retrieval parameter uncertainty [12]. Radiometric noise and the lack of spectral information limit the possibility to distinguish unequivocally among the various solutions that could fit the measurement vector (top of atmosphere cloud free BRF values) at a given level of confidence. This level of confidence, which depends on the size of the measurement vector, can change from pixel to pixel according to cloud and illumination conditions and to possible gaps in the input satellite images. Hence, each solution has a probability, depending on the number of degrees of freedom, i.e., the number of cloud free input BRFs. The most critical variables of such a system are then the aerosol optical depth (τ), the surface brightness and anisotropy. It is assumed that actual atmospheric situations fall within the range of standard aerosol discrete optical depth values between 0 < τ < 1 .
The surface brightness is simultaneously estimated during the retrieval process from a set of predefined solutions (parameters k, Θ), which describe the anisotropic properties of typical surfaces and the reflectance level ρ 0 , a non-constraint parameter estimated from the model inversion. A 10-day temporal compositing technique is applied to maximise the spatial coverage of cloud-free pixels. Finally, the retrieved surface state variables are used to derive the surface albedo quantities, together with their uncertainties. A probability, depending on the number of input cloud free measurements and uncertainty is also estimated. Due to the huge amount of calculations needed, in particular for the generation of a data record that spans over several decades, the number of possible states of the coupled surface-aerosol system is limited to a fixed number. The exploitation of precalculated Look-Up Tables (LUT) speeds up the RTM inversions. This retrieval scheme has been successfully applied to geostationary satellites from different operational agencies within the international Sustained, Coordinated Processing of Environmental Satellite Data for Climate Monitoring (SCOPE-CM) framework [25].

Release 2 Improvements
Compared to the previous release, the algorithm used for Release 2 allows for (i) ingestion of an external cloud mask data from a new cloud detection module based on the method presented by [20], (ii) allows for retrievals from both MVIRI and SEVIRI imagery, and (iii) allows retrievals over expanded spatial and temporal scales.

Cloud Removal
The most relevant difference between the two releases is the approach for cloud masking. In Release 2 the information on cloud free observations is obtained in two steps. Firstly, from the ingested cloud mask data generated with the method developed at the EUMETSAT Climate Satellite Application Facility (CM-SAF) [20,21]. Secondly, from the internal cloud masking method based on an analysis of the diurnal cycle of the top-of-atmosphere BRF for observations identified as cloud free in the first step [26]. In Release 1 no external cloud masking is exploited. The MVIRI input to the external cloud masking algorithms comprised of recalibrated infrared radiances [27], and vicarious calibrated visible reflectances [19]. The SEVIRI input to these algorithms comprised of operational calibrated infrared and vicarious calibrated visible reflectances [18]. The retrieved DHR30 in the instrument native VIS channel only using the internal cloud mask procedure for Release 1 was compared to the results of using both the internal cloud-masking procedures for Release 2. This comparison revealed that the new cloudmasking procedure detects much more clouds and, as a result, residual cloud contamination in Release 2 has drastically decreased. To verify the impact of removing clouds more effectively, two regions were analysed: one north (CAFR) ( Figure 3) and one south (SAFR) (Figure 4) of the Equator. The precipitation regime in these two regions is well known and is depending on the shift of the Intertropical Convergence Zone (ITCZ) and the African Monsoon Source: Climate and Ocean: Variability, Predictability and Change (CLIVAR) http://www.clivar.org/african-monsoon (accessed on 18 May 2021). More precipitation is expected from June to September north of Equator and from January-March south of the Equator. The presence of precipitation is associated with persistent cloud coverage. The impact of including or excluding the use of an external cloud mask on the retrieved albedo should have a clear signature matching in time and space the precipitation regime.
The VIS DHR30 retrieved over those two regions of Release 1 and Release 2 have been compared for the complete year 2001. As expected, the introduction of an external cloud mask removed almost entirely the clouds signature in the product. The remaining clouds are due to the limitation in the external cloud mask. The impact of removing clouds more effectively is clear in the period June-September north of the Equator (Figure 3) and in the period January-March south of the Equator (Figure 4). The number of retrievals decrease between the 20% and 60% in the season of permanent cloud coverage, while it remains almost constant for periods with a lower presence of clouds and in region such as the Sahara where the presence of clouds is minimal. The remaining retrievals are of higher quality, more reliable and not sohwing any artefact due to remaining cloud contamination.

Conversion to Broadband Albedo
In order to enable a comparison among different Meteosat satellites and among data records derived from different satellite instruments, the albedo quantities retrieved with GSA must be converted into a larger common spectral broadband interval (0.3-3.0 µm). The broadband albedo conversion is performed from the native instrument visible channel of MVIRI (Meteosat-2 to Meteosat-7) or the High-Resolution Visible (HRV) channel of SEVIRI (Meteosat-8 to Meteosat-10). Figure 1 shows the differences in the spectral bands of the used channels and this will result in a different spectral albedo value retrieval. As already mentioned, in this paper the DHR is used for analysis, but for the sake of completeness, the conversion is here extended to the BHR iso The conversion to broadband is done for MVIRI (Meteosat-2 to Meteosat-4) following the method described in [28] and for MVIRI (Meteosat-5 to Meteosat-7) and SEVIRI (Meteosat-8 to Meteosat-10) following the method described in [29]. The need for two methods is due to an issue in the nominal Sensor Spectral Response (SSR) that was leading to some temporal discrepancies in the observations from Meteosat-2 to Meteosat-4. Loew and Govaerts adopted an empirical method to derive the spectral coefficients. However, the relation between spectral and broadband albedo is the same for both methods and is given by a third-order polynomial (Equation (1)): 3 (1) where VAR is either DHR or BHR iso , each using a different set of conversion coefficients a-d for the DHR (Table 1) and BHRiso (Table 2) according to the Meteosat satellite number.

Albedo Data Record
The data record was generated exploiting MVIRI and SEVIRI imagery from MFG and MSG (see Table 3

Validation Datasets
This section presents the sites and the reference data used for the validation of the GSA Release 2 data record. For the temporal stability analysis six sites, including bare soil, urban, shrubland, and vegetation, were selected from the list of available Surface Albedo Validation Sites (SAVS) [33]. The selection criteria include coverage of different parts of the Meteosat disk and representativeness of the most common surface types. The locations of the selected sites are given in Table 4.  [34].
Where the blue-sky albedo refers to the albedo calculated under real-world conditions with a combination of both diffuse and direct lighting based on atmospheric and view-geometry conditions [35].
It is important to note that the compared quantities are not the same. The SAFARI comprises BHR (blue-sky albedo) retrievals, whereas GSA and MODIS comprise DHR30 (black-sky albedo) retrievals. The BHR contains a direct (85% under clear sky conditions) and diffuse (15% under clear sky conditions, depending on aerosol optical depth and Rayleigh scattering) component of albedo, while DHR30 only contains the direct component of albedo. The diffuse component is strongly depending on the aerosol content, which is expected to be most relevant over an urban site such as MONGU. However, this comparison provides a good qualitative indication of retrieval quality, considering that the ratio between the diffuse and direct radiation does not change drastically in time and space [36].

Results
In this section the assessment outcome for the GSA Release 2 is presented. It includes an evaluation of the retrieval consistency, an analysis on the temporal stability and finally a comparison against reference dataset.

Retrieval Consistency Assessment
In order to assess the retrieval consistency among the different missions, a comparison of the broadband DHR30 was done for an overlapping transect. A period in April 2005 was is shown, but similar results were obtained for other periods. In order to minimize the potential discrepancies due to a different viewing geometry, the transect was chosen at the longitude of 31.5 • . This longitude is exactly in the middle between the nominal locations of the 0 • mission (both MFG and MSG) and the IODC mission at 63 • . Figure 5 shows the very good agreement between the retrievals from the first (Meteosat-5 and -7) and second generation (Meteosat 8) of imagers. The DHR30 correlation is always higher than 0.98. Table 5 summarizes the statistics (correlation, the Mean Absolute Error (MAE) and Root Mean Square Error (RMSE)).

Temporal Stability Analysis
In order to assess the data record stability over time, DHR30 time series for the selected sites (see Table 3) were extracted. Surface albedo is supposed to be relatively stable with time over desert sites and decadal trends (linear increase or decrease over a period of 10 years) are expected to be lower than 1%. Significant decadal trends (>1%) can be expected for other sites, such as urban sites. In order to assess the decadal trend, the DHR30 broadband time series has been resampled averaging over the months. A trend-seasonal-residual decomposition [37] procedure, as shown in Equation (2), has been applied. The method includes an outliers handling, trying to limit their impact on the decomposition.
where Y t is the observed parameter (DHR30), T t is the trend component, S t is the additive seasonal component, and R t is the remainder component, for each observation at time t (t = 1. . . N number of observations). The Mean Trend Value is the mean value of T t over the full period. A linear regression (linear fit of T t ) on the trend component allowed estimating the decadal trend (Figures 6 and 7, upper plot). In order to confirm that any seasonality from the input time series has been removed, an Auto Correlation Function (ACF) plot has also been computed up to 24 lags (2 years, 24 months). The objective is to test the hypothesis that the time series is auto correlated due to a seasonal component. This component is present, as expected, for the MONGU urban site (Figure 7, bottom plot) in the input broadband DHR30 (red dots) array, but it disappears when looking at the trend component (green diamonds). The shaded areas indicates the upper and lower bounds of the 95% confidence level (only lags inside this area are statistically significant). The same analysis for the LIBIA_00001 desert site ( Figure 6, bottom plot) shows, as expected, no clear seasonality.
The above trend-seasonal-residual decomposition procedure has been performed for all selected sites (see Table 3). The decadal trend over more than thirty-six years and nine different platforms is close to 1% (LIBIA_00001, EGYPT_ONE, SOV, see Table 6) for all bright desert sites, the scrubland site (SKUKUZA) and forest site (BELMANIP_00025) as listed in Table 3. The MONGU site shows a decrease of about 2% per decade. For the MONGU site, the seasonal variations (see Figure 7) are due to surface dynamics that are linked to the precipitation regime and/or to minimal residual cloud contamination in the input cloud mask. Such variations are not present for the LIBIA_00001 site (see Figure 6), where the issue of potential minimal undetected clouds is much less relevant due to the much lower average annual cloud coverage of that site.  The corresponding 95% confidence regions are also displayed. Only lags inside those areas are statistically significant.

Validation against Reference Dataset
This section presents the results of the black sky albedo (DHR30) comparison over the six selected SAVS sites. In the time series, like the example in Figure 8, the outliers are found using a filter based on the median [38]. The RMSE and MAE are calculated excluding the outliers. Such outliers for GSA are most likely due to clouds undetected by the cloud mask. This hypothesis is supported by the fact that outliers have values higher than the yearly variation range. PERC is the ratio between the RMSE and the DHR30 average (AVG) expressed as percentage and calculated over the number of decadal (10 days) values available during a year (SAMP value with a maximum of 37, i.e., the maximum number of GSA values in a year). Each value represent the average over a 5 × 5 pixels around the location and the error bar is the standard deviation calculated over the same area. The location of the site on the Meteosat disk is shown in the plot. The results presented in Figure 8 and Table 7 confirm the good agreement between GSA and MODIS, with PERC values for DHR30 ranging between 2% and 7% for bright desert sites, and with value between 8% and 9% for the urban (MONGU), vegetation (BELMANIP_00025), and shrubland (SKUKUZA) sites. The comparison versus SAFARI (see Figure 9 and Table 8) shows a worse agreement, with PERC values of 21%. Although, the comparisons of SAFARI (BHR) versus MODIS and GSA (DHR30) show similar RMSE and PERC values, the MAE value of the comparison versus GSA is higher than the corresponding value for MODIS. The latter is an indication of the presence of residual clouds not screened out by the cloud mask, which is confirmed by the significantly higher AVG values of GSA than of MODIS (see Table 7). For the SKUKUZA site the differences between the GSA vs SAFARI and MODIS vs SAFARI are much smaller, with slightly higher RSME, PERC, and MAE values for GSA. The AVG values for GSA and MODIS are almost identical (see Table 8).

Discussion
Black sky and white sky albedo data records, together with associated uncertainties, have been generated for both MVIRI on the MFG and SEVIRI on the MSG geostationary satellites of EUMETSAT over a period of 36 years, from 1982 to 2017 for the main mission at 0 degrees longitude. In addition, a 20 years data record (1998-2017), covering the Indian Ocean has been generated from IODC data. The retrieval algorithm applied is the same for both platforms and for all three longitudinal orbit positions (0°, 57°E, 63°E ) [12][13][14]. The study has clearly demonstrated that for historical satellites, also for those where less calibration information is available, the GSA data record shows very good temporal stability and a good agreement with reference data. The temporal analysis of the data records confirms the homogeneity and stability of the Meteosat surface albedo time-series, with decadal trends over the desert sites (LIBIA_00001, EGYPT_ONE, SOV) and forest site (BELMANIP_00025) below 1%. In addition, the decadal trend over the shrubland site (SKUKUZA) is well below 1%. The urban site (MONGU), in southern Africa, reveals a trend larger than 2% with clear seasonal variations. However, this site by definition (urban site) is not supposed to show a stable behaviour over such long time range. The analysis of the spatial overlap ( Figure 5 and Table 6) shows very good consistency among the three data records (MFG 0 • , MSG 0 • and MFG IODC 63 • ), with correlation is close to 0.99 and the RMSE values around 7%. The validation against reference data revealed that the GSA agree well with MODIS and and slightly worse with SAFARI (Tables 7 and 8). The RMSE of GSA versus MODIS is less than 6% for the desert sites. The agreement is slightly worse for the urban (8%), shrubland (8.1%) and rain forest (8.9%) sites. The use of an external cloud mask, besides cloud screening based on an analysis of the daily cycle of the top-of-atmosphere BRF, significantly improved the quality of the GSA Release 2 data record compared to the MSA Release 1 data record. The cloud detection method is more complicated at high satellite view zenith angles (longer path through the atmosphere and high pixel distortion), in particular over region with permanent cloud overcast. Therefore, a post processing through a threshold comparison with a climatological background may help to remove the remaining noise in the surface albedo data due to undetected clouds. This could, for example, be achieved by using a climatological seasonal albedo data record, created from the GSA data record itself, in order to preserve the overall system energy [17]. Although such a post processing could remove good retrievals as well, it is preferable for a climate data record to have a limited amount of retrievals with a high quality, instead of a high amount of retrievals with a low quality. Some known limitations from the previous release remain. For instance the limited range of τ in the look-up tables might lead to an albedo overestimation due to the saturation of the aerosol contribution (τ > 1) in the coupled surface-scattering system. This effect has been confirmed by the analysis done for Release 1 [17] and it is also valid for Release 2 because this part of the algorithm has not been changed. There is potential for further improvement of the quality of the data record presented in this paper by replacing the input of operational calibrated visible imagery with input of recalibrated visible imagery using reconstructed spectral response functions [39][40][41]. In conclusion, GSA Release 2 contains both a spatial and temporal extensions and a clear improvement in quality compared to the previous release. The overall quality of the product is good and Release 2 offers a longer and a more reliable data record for usage in climate monitoring. Data Availability Statement: The data are accessible via the EUMETSAT Data Centre and can be easily found using the data record Digital Object Identifier (DOI) number. To access the data from EUMETSAT, a user need to register with the EUMETSAT Data Centre. For enquiries or feedback the user can contact the EUMETSAT User Service Helpdesk by email: ops@eumetsat.int.