Anisotropy Parameterization Development and Evaluation for Glacier Surface Albedo Retrieval from Satellite Observations

Abstract

Glacier albedo determines the net shortwave radiation absorbed at the glacier surface and plays a crucial role in glacier energy and mass balance. Remote sensing techniques are efficient means to retrieve glacier surface albedo over large and inaccessible areas and to study its variability. However, corrections of anisotropic reflectance of glacier surface have been established for specific shortwave bands only, such as Landsat 5 Thematic Mapper (L5/TM) band 2 and band 4, which is a major limitation of current retrievals of glacier broadband albedo. In this study, we calibrated and evaluated four anisotropy correction models for glacier snow and ice, applicable to visible, near-infrared and shortwave-infrared wavelengths using airborne datasets of Bidirectional Reflectance Distribution Function (BRDF). We then tested the ability of the best-performing anisotropy correction model, referred to from here on as the ‘updated model’, to retrieve albedo from L5/TM, Landsat 8 Operational Land Imager (L8/OLI) and Moderate Resolution Imaging Spectroradiometer (MODIS) imagery, and evaluated these results with field measurements collected on eight glaciers around the world. Our results show that the updated model: (1) can accurately estimate anisotropic factors of reflectance for snow and ice surfaces; (2) generally performs better than prior approaches for L8/OLI albedo retrieval but is not appropriate for L5/TM; (3) generally retrieves MODIS albedo better than the MODIS standard albedo product (MCD43A3) in both absolute values and glacier albedo temporal evolution, i.e., exhibiting both fewer gaps and better agreement with field observations. As the updated model enables anisotropy correction of a maximum of 10 multispectral bands and is implemented in Google Earth Engine (GEE), it is promising for observing and analyzing glacier albedo at large spatial scales.

1. Introduction

The shortwave surface albedo is the ratio of the hemispheric fluxes of the upwelling and the downwelling shortwave radiation (300–2500 nm). The albedos of snow and ice surfaces are determined by surface characteristics, such as grain size, impurity and water content, and properties such as surface structure and roughness, and vary with the wavelength of radiation and solar zenith angle [1]. Net shortwave radiation is the main driver of the glacier energy balance and largely controls the seasonal and daily variability of glacier melt energy [2,3,4,5,6,7]. Glacier surface albedo directly determines the net shortwave radiation at the surface, thus driving glacier melt, and can therefore be used as a proxy to estimate glacier mass balance [8,9,10]. Previous studies have demonstrated that remote sensing techniques are an efficient means to derive glacier albedo at different spatial and temporal resolutions and can provide an improved understanding of its spatio-temporal evolution over ground observations that are often limited in space and time [11,12,13].

Since satellite sensors receive the reflected radiance from both the land surface and the atmosphere at specific directions and in narrow spectral bands, albedo retrieval from remote sensing data generally requires three steps: (1) atmospheric correction, which removes atmospheric scattering effects from satellite radiance measurements; (2) anisotropy correction, which accounts for anisotropic reflection on the land surface and produces hemispherical reflectance (narrowband albedo); and (3) Narrowband-To-Broadband (NTB) conversion, to estimate the spectrally integrated shortwave albedo from multiple narrowband albedo values [14,15]. For glacier surfaces, steps 1 and 3 have been thoroughly addressed by other studies, e.g., Masek et al. [16] and Vermote et al. [17] for step 1; Knap et al. [18], Liang [19] and Stroeve et al. [20] for step 3. Step 2, however, has been addressed only for a few spectral bands of specific satellites, which therefore limits its large-scale application. Improving this anisotropy correction for glacier snow and ice and applying it to a variety of sensors are the main focuses of this study.

Several Bidirectional Reflectance Distribution Function (BRDF) parameterizations have been developed for albedo retrieval, e.g., Lucht et al. [21] and Reijmer et al. [22]. For Moderate Resolution Imaging Spectroradiometer (MODIS) albedo, the Ross Thick-Li Sparse model, proposed by Lucht et al. [21], is a well-known semi-empirical linear kernel-driven model and is used in one of the current MODIS albedo products (MCD43A3, Schaaf et al. [15]). However, this model requires a sufficient number of cloud-free observations within a 16-day window to construct a BRDF, leading to gaps in albedo time series, even occasionally on cloudless days. Such cloudy conditions are prevalent in mountainous, glacierized regions, limiting the analysis of spatio-temporal albedo variability on glaciers. For satellite images with lower repeat frequency, such as Landsat, several linear parameterization schemes have been used successfully to estimate glacier surface BRDF, for example with the Landsat 5 Thematic Mapper (L5/TM) band 2 (520–600 nm) and band 4 (760–900 nm) [23,24]. However, the relatively few available measurements of reflectance anisotropy over snow and ice surfaces have led past studies to often simply ignore surface anisotropic effects (e.g., Pope and Rees [25] and Pope et al. [26]). This results in a significant underestimation of glacier surface albedo, which can be up to 10% or ~0.1 for specific surface types, spectral bands and sun zenith angles [23,27,28]. One way to deal with this problem of accurately estimating snow and ice hemispherical reflectance is to use spectral parametric models of BRDF. Studies have also developed NTB conversions of L5/TM multispectral data [18,19], with better results achieved with parameterizations using more bands, which requires developing anisotropy correction models for more bands due to the spectral variability of snow and ice BRDF [28,29,30]. For example, Gatebe and King [31] acquired an airborne spectral BRDF dataset over Arctic sea ice and snow covering the whole shortwave spectrum (340–2300 nm) with high angular sampling resolution for this purpose.

Evaluation of remotely sensed albedo generally includes two methods. For high-resolution satellite albedo products (<50 m), such as L5/TM and Landsat 8 Operational Land Imager (L8/OLI) and Sentinel-2 MultiSpectral Instrument (S2/MSI), field observations at weather stations can be used to evaluate albedo estimates directly. For validation of coarser resolution albedo products, such as MODIS (500 m), near-contemporaneous high-resolution satellite products are sometimes used because of scale differences between field observation and the satellite footprint, a problem amplified by the heterogeneity of glacier surfaces. For example, Naegeli et al. [12,28] validated L8/OLI and S2/MSI albedos at two glaciers against field observations and investigated the spatio-temporal albedo variabilities of 39 glaciers in the Alps. Similarly, Wang et al. [24] validated L5/TM and Landsat 7 Enhanced Thematic Mapper Plus (L7/ETM+) against field observations and evaluated the MODIS snow albedo product (MOD10A1) against L5/TM products for five glaciers on the Tibetan Plateau in the 2000–2011 period. Their results showed that the mean absolute biases between high-resolution satellite albedo and field observations were close to ±0.05, and the MOD10A1 albedo was generally underestimated relative to the L7/ETM + albedo. However, many studies only used local ground observations of glacier albedo at a few sites in specific regions (e.g., High Mountain Asia or the Alps) due to the scarcity of such field data, while L5/TM, L8/OLI and MODIS albedos would rather need to be evaluated across a wider range of glacier surfaces for these products to be used effectively at the regional to global scale. Due to the aforementioned spatial scale difference, in the case of MODIS, such an evaluation requires assessments against intermediate spatial resolution sensors to account for local surface variability. Finally, albedo analyses of large areas or for repeated images require an albedo retrieval method that is also computationally suited to processing big datasets.

At present, the anisotropy correction models are applicable to a few spectral bands, restricting the improvement and application of glacier albedo retrieval. The objectives of this study are three-fold: (1) to develop and evaluate the performance of four anisotropy correction models for more bands for glacier surfaces (snow and ice) against the Gatebe and King [31] airborne spectral BRDF dataset; (2) to develop an updated albedo retrieval procedure using the best anisotropy correction model in Google Earth Engine (GEE) and evaluate its performance for L8/OLI and L5/TM data against Automatic Weather Station (AWS) measurements on glaciers around the world; (3) to retrieve MODIS albedo using the updated model and evaluate its accuracy relative to the MCD43A3 albedo product, AWS data and generated L5/TM and L8/OLI albedo. These improvements of albedo retrieval will pave the way for well-constrained glacier albedo datasets at the large scale.

2. Study Sites and Datasets

2.1. Study Sites

Based on available data, we compiled eight on-glacier AWS records to calibrate and evaluate the updated albedo retrieval method (Figure 1, Table 1). These glaciers are located in diverse mountain ranges around the world encompassing different climates from low (<2000 m a.s.l) to high elevation (>5000 m a.s.l). One glacier is located in Alaska (Figure 1a), one is in the Caucasus (Figure 1b), four glaciers are in the inner Tibetan Plateau and eastern Himalayas (Figure 1c–f) and two glaciers are in the Andes (Figure 1g,h). The AWSs are all located in ablation areas, with the exception of the Mera Glacier AWS. The glaciers vary considerably in size, shape, topography and surface heterogeneity (Figure 1, Figure A1 (Appendix A)). This set of AWSs is not comprehensive but is suited to evaluating albedo retrieval across a wide variety of conditions. The AWSs were situated on different glacier surfaces of snow, ice or a mixture of ice and debris, so different surface types were accounted for in our satellite albedo retrieval validation (Figure A1).

2.2. Datasets

Snow and ice BRDF data: High-quality BRDF measurements on mountain glacier surfaces are rare, and most available measurements sample a limited spectral range at specific view angles [23,32,33]. Instead, we analyzed two airborne BRDF measurement datasets made over snow and sea ice [31]. The snow BRDF data are from the Arctic Research of the Composition of the Troposphere from Aircraft and Satellites (ARCTAS) experiment conducted on 6 April 2008. The measurements were taken on the coast of the Arctic Ocean at Elson Lagoon, Alaska and include 10 individual bands, covering visible, near-infrared and shortwave-infrared spectral bands. The sea ice BRDF data are from the Arctic Radiation Measurements in Column Atmosphere-Surface System (ARMCAS) experiment conducted on 8 June 1995 over the Elson Lagoon, Alaska and measurements include 6 individual bands, covering visible and near-infrared spectral ranges [31,34]. The 6 BRDF bands for snow and 3 BRDF bands for ice were selected and applied to the respective parameterization taking into account the specific spectral features of snow and ice, although the spectral coverage of the dataset we applied was not optimal for ice. Details on the spectral coverage of the bands applied in this analysis are given in Table A1, including central wavelength and spectral width. Note that the coast and sea ice at this location are generally covered by a thick, dry snowpack in April and that the sea ice exhibits a mixture of bare ice and shallow ponds in June [35]. The view zenith angle ranged from 0 to 90° in 0.5° intervals, the relative azimuth angle between the sensor and sun ranged from 0 to 360° in 1° intervals, and the solar zenith angle ranged between 61.9° and 70.9° for ARCTAS (snow) and between 49.4° and 57.6° for ARMCAS (sea ice). We used these data to extract and validate anisotropy corrections, assuming that the BRDF of sea ice is similar to that of glacier ice since both are frozen water and have a similar density (~900 kg/m³) and melting point [36]. This assumption ignores differences in salinity, roughness and optical penetration between sea ice and glacier ice [36], and therefore implies some uncertainties in the derived anisotropy correction, which we consider in the evaluation of our results. As a result, the developed anisotropy correction models covered more shortwave bands than existing models limited to L5/TM band 2 and band 4 [23], which enables glacier albedo retrieval using multispectral satellite data.

Field observations of glacier albedo: Incoming and reflected shortwave radiation were measured by Kipp and Zonen CNR1 or CNR4 Net Radiometers at AWS locations on eight glaciers (Table 1, Figure 1 and Figure A1). The albedo was derived directly as the ratio of reflected and incoming shortwave radiation measured by the sensor. Albedo values > 1 or < 0 were set to 1 or 0, respectively. Some typical problems of radiometers, such as sensor tilt, sensor icing or sensors being covered in snow, can lead to large albedo uncertainties. We therefore checked that fresh snow albedo values were regularly above 0.9 to preliminarily ensure that the radiometers were working normally. The field measurements of irradiance and reflected hemispherical radiance were used as provided by the operator of each observatory and no further correction was applied. The basic information of each AWS, including location, observation period and elevation, is shown in Table 1. The AWSs were installed in the middle of the glaciers to ensure radiometers correctly captured the irradiance on and the radiation reflected by the glacier surface, while avoiding as much as possible shadow and other terrain effects (multiple reflection) (Figure A1). In order to ensure temporal consistency between the satellite’s overpass and the albedo field observations, the mean albedo between 11:00 and 13:00 (all times are local solar time) was used to evaluate the satellite-derived albedo.

Satellite data: Land surface reflectance data from L8/OLI, L5/TM (Surface Reflectance Tier 1) and MODIS carried by both Terra and Aqua satellites (MOD09GA and MYD09GA Collection 6) were used to retrieve glacier albedo with spatial resolutions of 30 m and 500 m, respectively. Six spectral bands of each satellite sensor were used to retrieve albedo (Table A1). The local overpass times of L8/OLI and L5/TM are ~10:30, and for MODIS onboard Terra and Aqua satellites are ~10:30 and ~13:30, respectively. Snow and ice albedo, measured at the AWS, are known to exhibit diurnal variations [37,38]. According to the study of Wang and Zender [38] in the Arctic and Antarctica, the albedo at noon best represents the daily albedo. In addition, the albedo in the central hours of the day plays a key role in driving the glacier surface energy balance, because of the high irradiance in this part of the day. In this study, errors (often underestimation) are inevitable but is acceptable when using instantaneous albedo retrieved by satellite data to estimate the daily albedo. MCD43A3 is a standard MODIS albedo product and can provide both white-sky (in the absence of direct radiation) and black-sky (in the absence of diffuse radiation) albedo. In this study, only black-sky albedo was used to compare with our MODIS albedo product because diffuse irradiance is small and negligible at high elevation and under clear-sky conditions [39]. All albedo retrievals were processed in GEE, where the data are available from January 1984 to May 2012 for L5/TM, from March 2013 to present for L8/OLI and from February 2000 to present for MODIS. ALOS World 3D-30 m (AW3D30) version 2.2 is a global Digital Surface Model (DSM) dataset with a horizontal resolution of approximately 30 m generated by Advanced Land Observing Satellite (ALOS) stereo images during 2006–2011 [40]. This DSM was used to derive glacier surface slope and aspect for the topographic correction (Equations (3a) and (3b)).

remotesensing-13-01714-t001_Table 1 Table 1

Basic information of albedo measurements from the AWSs of the eight study glaciers. * indicates that shortwave radiation was only recorded in June–September at the Djankuat Glacier. The numbers of available scenes of MODIS were not presented because MODIS albedos were retrieved by available pixels and we used nearly all daily scenes.

Region	Glacier	Longitude (°)	Latitude (°)	Start Date, End Date	Elevation (m a.s.l.)	Satellite Data (Numbers of Available Scenes)	References
Alaska	McCall	−143.85	69.32	05/01/2007, 31/12/2014	1720	L5/TM (10), MODIS	Troxler et al. [41]
Caucasus	* Djankuat	42.76	43.20	15/06/2007, 01/09/2017	3000	L5/TM (15), L8/OLI (12), MODIS	P Rets et al. [42]
Inner Tibetan Plateau and eastern Himalaya	Zhadang	90.65	30.47	01/01/2011, 31/12/2014	5660	L5/TM (15), MODIS	Zhang et al. [2,3]
	Parlung No.4	96.93	29.25	05/01/2012, 20/09/2018	4600	L8/OLI (36), MODIS	Yang et al. [43]
	Yala	85.62	28.23	08/05/2016, 19/11/2019	5330	L8/OLI (28), MODIS	ICIMOD
	Mera	86.88	27.72	01/01/2017, 12/11/2019	5769	L8/OLI (11), MODIS	GLACIOCLIM
Andes	Artesonraju	−77.64	−8.96	13/03/2004, 13/05/2013	4797	L5/TM (13), MODIS	Winkler et al. [44]
	Zongo	−68.14	−16.28	06/08/2004, 31/08/2019	5050	L5/TM (24), L8/OLI (18), MODIS	GLACIOCLIM

3. Methods

Since we used atmospherically corrected data for L5/TM, L8/OLI and MODIS from the U.S. Geological Survey [16,17], images only needed to be processed for the last two steps of the albedo retrieval process, the anisotropy correction and NTB conversion. We developed the updated retrieval method in GEE and generated 30 × 30 m L8/OLI, L5/TM and 500 m × 500 m MODIS albedos for the eight study glaciers. MODIS daily albedo was estimated using both Terra and Aqua satellites, i.e., averaging the albedo retrieved from Terra and Aqua as the daily albedo if both were available, otherwise using only one.

3.1. Anisotropy Correction of the Glacier Surface

In order to retrieve glacier albedo from satellite multispectral data, a relationship needs to be built between narrowband albedo (${\alpha }_i$), i.e., the hemispherical in-band reflectance, and directional reflectance ($r$) in the same spectral band and at a specific view angle, i.e., the land surface reflectance observed by satellite. Following Greuell and De Ruyter De Wildt [23], we used:

(1)${\alpha }_i=r-\tilde{f}{\alpha }_i{}^{1-m}$

where $\tilde{f}$ depends on the BRDF parameterization (Table 2) and sun zenith angle. Table 2 shows the four selected BRDF parameterizations developed for snow or glacier ice surfaces, which were derived separately by Reijmer et al. [22], Greuell and De Ruyter De Wildt [23], Warren et al. [45] and Knap and Reijmer [46], and all presented in Reijmer et al. [22]. $m$ is a free parameter and needs to be numerically estimated using BRDF measurements. In this study, $m$ is set to 1, which assumes that the BRDF is independent of surface material and is a good compromise between efficient computation and accuracy of the results [23,46]. According to Greuell and De Ruyter De Wildt [23], setting $m$ to 1 results in a reduction in the total variance of 3% for L5/TM band 2, and 2% for L5/TM band 4, which is acceptable to reduce the computational load of data analysis. Thus Equation (1) reduces to:

(2)${\alpha }_i=r-\tilde{f}$

The $\tilde{f}$ can be determined by hemispherical integration of each BRDF parameterization, and more details can be found in Greuell and De Ruyter De Wildt [23]. $\tilde{f}$ was only developed according to BRDF parameterization P2 (Table 2) and validated for glacier ice surface in their study. Here, we extended this approach to the other three BRDF parameterizations and calculated and validated $\tilde{f}$ for each one for both glacier snow and ice surfaces, i.e., P1, P3 and P4 in Table 2.

Another area to consider is that surface slope and aspect can alter the angles in satellite-surface-solar geometry and then affect albedo retrieval [47,48]. Prior to applying an anisotropy correction, it is therefore necessary to correct these angles to the equivalent horizontal surface because glaciers are rarely flat and can exceed 40 degrees surface slopes, e.g., in the accumulation area of Zongo Glacier (Figure A1h). Two topographic corrections depending on surface slope and aspect were adopted for view and sun zenith angles corrections [47,49]:

(3a)$\cos {\theta }_{vc}=\cos a\cos {\theta }_v+\sin a\sin {\theta }_v\cos \left(b-{\phi }_v\right)$

(3b)$\cos {\theta }_{sc}=\cos a\cos {\theta }_s+\sin a\sin {\theta }_s\cos \left(b-{\phi }_s\right)$

where${\theta }_v$ and ${\theta }_{vc}$ are the view zenith angles of satellite before and after correction,${\theta }_s$ and ${\theta }_{sc}$ are the sun zenith angles before and after correction, ${\phi }_v$ and ${\phi }_s$ are the satellite and sun azimuth angles, and$a$ and $b$ are glacier surface slope and aspect generated from the ALOS DSM data.

Glacier albedo can thus be retrieved using satellite observations of land surface reflectance and the $\tilde{f}$ models in Table 2. Since the BRDFs of the snow and ice surfaces are spectral functions, the $c_1$, $c_2$, $c_3$ and ${\theta }_c$ coefficients of $\tilde{f}$ need to be estimated and validated with BRDF measurements for different spectral bands.

To fit and validate these $\tilde{f}$ models, the measured $\tilde{f}$-values (${\tilde{f}}_m$) were directly calculated by airborne BRDF measurements ($r_{BRDF}$) and measured albedo (${\alpha }_m$) by Equation (4):

(4)${\tilde{f}}_m=r_{BRDF}-{\alpha }_m$

We divided the airborne BRDF measurements into two sample sets, each covering a wide range of solar zenith angles. The first group was used to estimate $c_1$, $c_2$, $c_3$ and ${\theta }_c$ by fitting $\tilde{f}$ to the measurements by Least Squares Minimization. The second group was used to validate each model results. We firstly sorted BRDF measurements based on increasing solar zenith angle, and then selected one out of three or four measurements for validation and the rest were used for fitting. A total of 19 and 9 distinct BRDF measurements were used to respectively fit and validate the $\tilde{f}$ models for snow, and 26 and 9 were used to respectively fit and validate the models for ice. We finally used narrowband albedo with the best-performing anisotropy correction model to generate our final L5/TM, L8/OLI and MODIS broadband albedo products.

3.2. Narrowband to Broadband Albedo Conversion

Broadband albedos can be empirically derived from the narrowband albedos. For L5/TM data, we converted narrowband albedo to broadband albedo following two commonly used conversions [18,19] and evaluated their accuracy against field observations. Based on field observations of broadband and L5/TM band 2 and band 4 narrowband albedos on the glacier surface, Knap et al. [18] developed a conversion based on multiple linear regression analysis, hereafter referred to as the Knap method (${\alpha }_{Knap}$):

(5)${\alpha }_{Knap}=0.726b_2-0.322b_2{}^2-0.015b_4+0.581b_4{}^2$

A well-known problem for NTB albedo conversion is the saturation of the L5/TM green band ($b_2$) on fresh snow surfaces, and in this case another conversion is used for L5/TM albedo retrieval following the Knap method [24,48]:

(6)${\alpha }_{Knap}=0.782b_4+0.148b_4{}^2$

Liang [19] used a radiative transfer model to simulate surface reflectance on different land surface types, and then developed conversions for L5/TM and L7/ETM+ data, hereafter referred to as the Liang method (${\alpha }_{Liang}$):

(7)${\alpha }_{Liang}=0.356b_1+0.130b_3+0.373b_4+0.085b_5+0.072b_7-0.0018$

In Equations (5)–(7), $b_i$ (i = 1, 2, 3…) represent narrowband albedo from the i-th L5/TM band (Section 3.1), the band numbers were adjusted accordingly for the L8/OLI. Most previous studies used the Knap method with the parameter values of P1 (snow) from Greuell and De Ruyter De Wildt [23] and P2 (ice) from Reijmer et al. [22] to retrieve narrowband albedos. In this study, we used our updated best-performing anisotropy correction models to retrieve narrowband albedo for the Liang method only, while using the parameter values directly from previous studies for the Knap method, so that we could compare our results with previous studies.

Similarly, Liang [19] also developed NTB conversion for MODIS data by surface in-band reflectance and suggested using it over ice surface (${\alpha }_{MODIS-ice}$):

(8)${\alpha }_{MODIS-ice}=0.160b_1+0.291b_2+0.243b_3+0.116b_4+0.112b_5+0.081b_7-0.0015$

Over snow surfaces, which show high reflectance, we used the conversion developed by Stroeve et al. [20], which showed better accuracy (${\alpha }_{MODIS-snow}$):

(9)${\alpha }_{MODIS-snow}=0.1574b_1+0.2789b_2+0.3829b_3+0.1131b_5+0.0694b_7-0.0093$

In Equations (8) and (9), $b_i$ (i = 1, 2, 3……) represent narrowband albedo from the i-th MODIS band (Section 3.1). Other than the beforementioned green saturation solution for the Knap method, all narrowband albedos were set to 1 when saturation happened. The ice albedo in the shortwave-infrared bands ($b_5$ and $b_7$ in L5, $b_6$ and $b_7$ in L8, $b_7$ in MODIS) was replaced by the bi-conical band reflectance observed by a space-borne imaging radiometer since the airborne BRDF datasets for ice do not cover these wavelengths, i.e., we could not build the anisotropy correction for these bands. We acknowledge that our replacement approach may lead to inaccurate estimates of the ice shortwave-infrared albedo. On the other hand, the impact is limited since both ice reflectance and albedo are approximately equal to zero in the ~1600 nm and ~2100 nm bands [32,50]. For example, the albedos of snow and ice are <0.001 in these two bands [28], which leads to maximum underestimations of 0.000157 in the final broadband albedo. The spectral bands of the airborne BRDF datasets and the L5/TM, L8/OLI and MODIS sensors are shown in Table A1.

3.3. Glacier Surface Classification and Albedo Validation

In order to classify the glacier surface as either snow or ice, we used a Normalized Difference Snow Index (NDSI) threshold method, with thresholds of 0.45 for L5/TM and L8/OLI and 0.4 for MODIS, respectively, following Girona-Mata et al. [51] and Härer et al. [52].

The final albedos generated from NTB conversions were validated at the eight glaciers. We first used the AWS albedo measurements to validate L5/TM and L8/OLI albedo retrieved with the Knap NTB conversion combined with the parameters directly from previous studies and Liang NTB conversion combined with our updated best-performing anisotropy correction models, for all cloud-free overpasses during the AWS observation period. We then selected the best-performing L8/OLI and L5/TM albedo retrievals to validate MODIS albedo retrievals combined with our updated best-performing anisotropy correction models, respectively. In order to reduce the effect of differences in spatial resolution, we compared the MODIS albedos with the average albedo value of 16 × 16 Landsat pixels around the MODIS pixel center. The validations of L5/TM and L8/OLI albedos were performed at the pixel overlapping each AWS position, while the validations of MODIS retrievals were performed over all pixels of each glacier. We performed these two evaluations based on three commonly used performance metrics, i.e., bias, Mean Absolute Error ($\mathrm{MAE}$) and Root Mean Squared Error ($\mathrm{RMSE}$).

4. Results

4.1. Evaluation of the Anisotropy Corrections

To evaluate the anisotropy correction models (Table 2) for snow and ice surfaces, we compared the $\tilde{f}$-values estimated (${\tilde{f}}_e$) according to Table 2, after having determined the $c_1$, $c_2$, $c_3$ and ${\theta }_c$ coefficients by fitting the calibration BRDF data, with the ${\tilde{f}}_m$ calculated directly from the BRDF data by Equation (4) in both calibration (parameter estimation) and validation (evaluation) BRDF data. This comparison was done for all the spectral bands for snow and ice. Three metrics, i.e., Pearson correlation coefficient (R), Mean Difference (MD) and standard deviation of the residuals (Std) between ${\tilde{f}}_e$ and ${\tilde{f}}_m$, were applied. We first calculated these three metrics for all spectral bands for each model, and then averaged them for each model (Table 3). Overall, the performances of the parameterization schemes P1, P2 and P3 are similar, and better than P4 (Table 3). Parameterizations P1–P3 give high R for both snow (0.93) and ice (0.79), with small positive and negative biases for snow and ice, respectively. In order to test the albedo differences retrieved by P1–P4, taking the snow surface on 26 June 2014 and the ice surface on 10 June 2013 on the Parlung No.4 Glacier as examples, we separately retrieved L8/OLI band 2 and band 4 snow and ice albedos using P1–P4 and calculated the difference due to these four parameterizations for each band. The results showed that the narrowband albedos retrieved by P1–P4 were very similar, and the differences were in a magnitude of 10⁻³, which is very small compared to snow albedo (>0.6) and ice albedo (>0.2). As a result, we still selected P1, the best-performing scheme for snow, to retrieve snow albedo from satellite data, which is consistent with the study by Reijmer et al. [22]. Similarly, and also since ARMCAS does not collect the ~550 nm band BRDF (key for ice albedo retrieval), we selected P2 instead of P3 (the best-performing scheme for ice) to perform our ice anisotropy correction, which is consistent with Greuell and De Ruyter De Wildt [23] so that we could directly use their parameters. According to the uncertainty experiments above, our choice has a very limited impact on the results.

Table 4 presents the resulting anisotropy correction models for snow (P1) and ice (P2). Since the spectral settings of the airborne BDRF are different for snow and ice, the counts and spectral ranges of their models are different. Similar with the evaluation for four correction models, we calculated three metrics (R, MD and Std) between ${\tilde{f}}_e$ and ${\tilde{f}}_m$ in each spectral band for P1 (snow) and P2 (ice). In total, 10 and 6 individual narrowband anisotropy correction models were built using airborne BRDF measurements for different satellite bands and separately covered 330–2260 nm for snow and 462–1281 nm for ice (Table 4). For snow, most Rs are above 0.90, while the ranges of MD and Std are 0.0001–0.0093 and 0.012–0.016. For ice, these ranges are 0.70–0.89, −0.0176–0.0133 and 0.043–0.131 for R, MD and Std, respectively. In general, the models for visible and near-infrared bands performed better than those for shortwave-infrared bands. In addition, higher R and lower MD and Std in snow correction models indicated that their accuracies are overall higher than that of sea ice (Table 4). The reason could be related to the BRDF data quality caused by the different surface properties of snow and sea ice. The quality of snow BRDF data were better because the snow surface was uniform and flat on the Elson Lagoon beach in April, while the quality of sea ice BRDF data was worse because the sea ice surface was a mixture of bare ice and shallow ponds and rough in June (see Section 2.2). Therefore, the performance of the anisotropy correction for snow was higher compared with sea ice.

4.2. Accuracy of L8/OLI and L5/TM Albedo

After applying both the BRDF and NTB corrections, we calculated the Bias, MAE and RMSE of L8/OLI and L5/TM albedo with respect to the AWS measurements in order to evaluate their accuracies (Table 5). The validation of Landsat albedo at each AWS site reveals two key findings. First, both the Knap and Liang methods retrieve albedo well for L8/OLI data, but the overestimation applies to the Liang method, while there is no significant difference when applying the Knap method (Figure 2 and Figure 3, Table 5). Indeed, at all sites the mean biases are 0.00 for the Knap method and 0.02 for the Liang method for the L8/OLI albedos, while they are 0.06 and 0.13 for the L5/TM albedos. Second, the Liang method with the updated anisotropy correction generally performs better than the Knap method in retrieving albedo from L8 data, but worse for L5/TM data (Figure 2 and Figure 3, Table 5). For L8/OLI albedo, the mean (range) MAE and RMSE of the Liang NTB estimation are 0.06 (0.03–0.10) and 0.07 (0.03–0.11), while these values are 0.07 (0.03–0.10) and 0.09 (0.06–0.13) for the Knap NTB conversion. For L5/TM albedo, the Liang NTB conversion overestimates albedo at all sites, especially for measured albedos > 0.5 (Figure 3b). The mean bias of the Liang NTB albedo is 0.13 against 0.06 for the Knap approach (Table 5). Since albedo retrieval heavily depends on local conditions on the glaciers, the comparisons between the two methods need to be done on more glaciers to evaluate overall performance.

4.3. Performance of Our MODIS Albedo Product and MCD43A3

Based on the results of Section 4.2, we used the L8/OLI albedo estimated with the Liang NTB conversion and the L5/TM albedo estimated with the Knap NTB conversion to evaluate our MODIS albedo retrievals (anisotropy correction as explained in Table 4, NTB conversions in Equations (8) and (9)), hereafter referred to as the MODIS/Ren and the MCD43A3 albedo products (Figure 4 and Figure 5, Table 6). The results show that the mean biases of both albedo products were −0.10 and −0.16, respectively, compared with L8/OLI and −0.04 and −0.11, respectively, compared with L5/TM. The maximum biases at all glacier sites were −0.17 for MODIS/Ren and −0.29 for MCD43A3, indicating that MODIS-derived albedos are lower than the L5/TM and L8/OLI albedos for most glaciers. Overall, MODIS/Ren is closer to the L5/TM and L8/OLI albedo than MCD43A3, with a smaller mean MAE (0.12 vs. 0.16 for L8/OLI, 0.12 vs. 0.13 for L5/TM) and RMSE (0.14 vs. 0.18 for L8/OLI, 0.14 vs. 0.16 for L5/TM) although the agreement is not uniformly good for either albedo products (Table 6). In addition, the clearest difference was observed in the high albedo range, i.e., snow surface albedo, where MCD43A3 often significantly underestimated albedo compared to MODIS/Ren, such as in the Mera and Zongo glaciers (Figure 4).

The L8/OLI and MODIS albedos also differ in terms of spatial variability, as can be observed, for example, over Parlung No.4 on 16 December 2014 (Figure 6). The three products (L8/OLI, MODIS/Ren and the MCD43A3) demonstrated a similar albedo spatial variability, i.e., albedo increased from low to high elevation. However, MODIS/Ren and the MCD43A3 albedos both displayed lower albedo and less spatial variability, while L8/OLI albedos were higher and had a greater spatial variability also when resampled to the pixel resolution of MODIS (Figure 6b).

We also directly compared daily albedo from AWS measurements and the two MODIS products to evaluate their performances in deriving albedo temporal evolution (Figure 7). Although there is a systematic bias between AWS measurements and the MODIS albedos, both products can capture seasonal albedo variability when enough albedo values are retrieved, for example for some periods at the Parlung No.4 and Mera Glaciers (Figure 7d,f). However, MODIS/Ren can better capture short-term and seasonal albedo fluctuations, such as during the monsoon season (May–August) over the Mera Glacier (Figure 7f), while MCD43A3 is heavily smoothed and seems not able to reproduce the in-situ derived albedo evolution very well. Furthermore, MODIS/Ren has a longer time sampling interval since the albedo can be retrieved as long as there are no clouds, while MCD43A3 exhibits data gaps for some cloudless days due to insufficient valid observations during a 16-day window to calculate BRDF parameters. On Zhadang Glacier, for example, there were almost no observations in the period of October–March in 2012–2014 in MCD43A3, while MODIS/Ren could capture seasonal variability consistent with the AWS observations (Figure 7b). In addition, in order to compare with MCD43A3 at the 16-day interval, taking Parlung No.4 Glacier as an example, we calculated the 16-day moving average of the MODIS/Ren albedos and compared it with field observations and the MCD43A3 albedo (Figure A2). The results showed that the MODIS/Ren albedos were higher than the MCD43A3 albedos, especially for the large albedo fluctuations caused by snowfall, and also agreed better with field observations in terms of albedo evolution during the 2013–2015 period. This is partly attributed to the better representativeness in temporal variability by the MODIS/Ren albedo.

5. Discussion

5.1. Limitations of the Updated Albedo Retrieval Method

The updated albedo retrieval method has a number of limitations due to: (i) differences between the glacier and sea ice BRDFs; (ii) the BRDF measurements being limited to a small range of solar zenith angles; (iii) spectral differences between the airborne BRDF datasets and the satellite bands; (iv) temporal variability of glacier surfaces; (v) band saturation and (vi) NTB conversion for L8/OLI data. Here follows a discussion of these limitations. (i) We used sea ice BRDF measurements from ARMCAS [31] to build the updated parameterization of glacier ice BRDF, as no other dataset was available, and assumed that the two types of ice have similar properties and thus albedo. This assumption ignores variations in salinity and surface properties between sea and glacier ice [36], but our validation scheme shows that despite these physical differences the updated method is able to retrieve albedo values that are consistent with field observations. (ii) An additional limitation of the BRDF measurements used is that they only span a small range of solar zenith angles and surface slopes and aspects, which makes their application at the global scale less straightforward because the updated anisotropy correction models depend on these angles. Even though our results are in reasonable agreement with field observations under these limitations, uncertainties in mountain glacier albedo could likely be reduced using specific glacier ice BRDF experiments over a larger range of solar zenith angle and surface topography. (iii) The spectral coverage of the airborne BRDF data, OLI, TM and MODIS bands is different (Table A1) and these differences can affect the retrieved narrowband albedo. Since the albedos of ice and snow in the visible bands are high and account for large weighting coefficients in the NTB conversion, a small difference in spectral coverage could lead to errors in the final broadband albedo retrieval. (iv) For this study, we used the AW3D30 DSMs derived from satellite images ranging from 2006 to 2011, but for long-term albedo retrieval one also needs to take into account glacier surface topography accuracy, as surface slopes and aspects from DSMs can be inaccurate, which has a direct impact on albedo retrieval. (v) Satellite visible bands can saturate over snow surfaces due to high reflectance. Since we set the narrowband albedo to 1 in the case of saturation, albedo is likely to be overestimated when this occurs, which is evident in our results (Table 5). Compared with L8/OLI, L5/TM visible bands are very easily saturated because of a smaller dynamic range, and the mean fractional abundance of saturated area in the L5/TM and L8/OLI visible bands (blue, green and red) were 25.0% (52.5%, 6.2% and 16.1%) and 2.0% (1.1%, 2.2% and 2.6%) according to the tests on the Djankuat and Zongo glaciers. The experiments on the Zongo Glacier showed that band saturation of L5/TM data can lead to overestimate by ~0.007 for the final broadband albedo. The lower accuracy of L5/TM albedo estimated by the Liang NTB weighing scheme together with the updated anisotropy corrections can also be explained by easier band saturation of the L5/TM bands 1 (blue band) and 3 (red band) over snow surfaces [24,48]. Therefore, we recommend using the Knap method (Equations (5) and (6)) to retrieve albedo from L5/TM data. (vi) Since the spectral settings of L5/TM and L8/OLI are different, especially in the near-infrared band, directly using NTB conversion for L5/TM to retrieve the L8/OLI albedo may cause errors [28].

Two metrics were used to quantify the uncertainties caused by the limitations mentioned above. On the one hand, we used Standard Error ($SE$) for the sample mean to estimate the uncertainties of different albedo retrieval methods:

(10)$SE=STD/\sqrt{N}$

where $STD$ is the standard deviation of the albedo differences between evaluation and reference observations and $N$ is the number of observations. Like albedo validation in Section 4.2 and Section 4.3, AWSs measurements were used to calculate $SE$s of L5/TM and L8/OLI albedos, while aggregated L5/TM and L8/OLI albedos were used to calculate $SE$ of MODIS albedos. We calculated the $SE$ of available observations on observed glaciers for each method. The results showed that the $SE$ were 0.007 for the Liang method and 0.008 for the Knap method with L8/OLI, 0.017 and 0.012 with L5/TM data. Regarding the MODIS albedos, the $SE$ was 0.003 when using L8/OLI albedos and 0.005 for L5/TM albedos, respectively. On the other hand, the histograms of the albedo differences (Figure A3) showed that 36.0% of the L8/OLI albedo values obtained with the Liang method differed by less than 0.04 from the AWS albedos, and 27.9% for those obtained with the Knap method. These values were 21.0% and 35.5% with L5/TM. For MODIS, 44.3% (43.3%) of the albedo values differed by less than 0.10 from the L8/OLI (L5/TM) albedos. In both cases there were smaller uncertainties for L8/OLI data with the Liang method, but larger for L5/TM data since the Liang method did not account for band saturation, which is especially frequent with L5/TM.

According to the aforementioned limitations, the albedo retrieval method can be improved by three aspects in the future: (1) by collecting high-quality BRDF measurements over different glacier surface types, such as snow with varying density, dirty ice and cryoconite, and debris-covered glacier surfaces in a wider range of solar zenith angles and surface slope and aspect, and for more spectral settings [24,26,28]; (2) by generating high spatial resolution and accurate glacier DSM products, which are very important for high-resolution and long-term albedo retrieval; (3) by developing NTB conversions that better account for band saturation and L8/OLI band settings, and by estimating the coefficients of NTB conversions taking into account the spectral distribution of irradiance at Bottom-Of-Atmosphere (BOA) and the effect of surrounding terrain on irradiance at the observed target (pixel) [53].

5.2. Evaluation of the Albedo Products

Overall, there is a reasonable agreement between our L8/OLI, L5/TM and AWS albedos (better for L8/OLI than for L5/TM) and between the MODIS and L8/OLI, L5/TM albedos, but some errors are still apparent. For the L8/OLI and L5/TM albedo, validation against field observations can introduce uncertainty due to heterogeneous and temporally variable glacier surfaces [24,25]. The Landsat albedo values represent the mean value within a 30 m × 30 m pixel, while field observations are point measurements with a small radiometric footprint (<100 m² depending on sensor height) on a glacier surface. Consequently, subpixel effects, especially in the vicinity of the AWS, could lead to apparent disagreement between L8/OLI, L5/TM and AWS albedos, and similarly between MODIS and AWS albedos. Subpixel effects depend on local conditions on the glacier surface, i.e., on a homogeneous surface such as fresh snow, subpixel effects could be small, but large on a heterogeneous surface, such as mixture of snow, ice and debris in the ablation zone. Since the AWS are often installed at lower elevation in the heterogeneous region of ice and debris because of access and ground stability, the AWS albedo is generally lower than that retrieved by satellite data. Indeed, heterogeneity of surface types and properties were apparent at the AWS locations found during field visits, especially for the Parlung No.4 (Figure A1d) and Zongo (Figure A1h) glacier surfaces, which could explain low accuracy of L8 albedo at these sites (Table 5). The results of Wang et al. [24] indicated similar challenges, leading to large errors of L5/TM albedo in the ablation zone of Laohugou No. 12 Glacier in their study.

The evident contrast between MODIS and higher resolution albedo products in our results has also been indicated by previous studies [20,24,26,54]. The spatial heterogeneity of glacier surface facies and local-scale complexity of terrain could explain this pattern. First, glacier surfaces are frequently covered by a mixture of surfaces (snow, dirty ice, ice, debris and melt water), which not only have variable albedos but also different anisotropic properties. In this study, we performed one anisotropy correction in each MODIS pixel, while >200 anisotropy corrections were performed in each MODIS pixel extent for the L5/TM and L8/OLI albedo retrieval, and this effect can be more pronounced on complex and heterogeneous surfaces and lead to albedo difference [26]. Second, since surface topography and roughness can alter solar-surface-satellite geometry and irradiance very locally, Landsat instruments can capture these changes better than MODIS, which inevitably leads to albedo differences [53,55,56,57]. However, unlike the validation of MCD43A3 in Greenland and large ice caps in Iceland, both our study and that of Pope et al. [26] and Wang et al. [24] showed that MODIS underestimated albedo on mountain glaciers. Two reasons probably caused this underestimation. On the one hand, the MODIS land surface reflectance is lower than L5/TM and L8/OLI values, especially in the high reflectance region, i.e., snow surface. For example, the surface reflectance observed with MODIS/Terra and MODIS/Aqua in visible bands were apparently lower than those observed with L8/OLI with differences of approximately −0.10 over Parlung No.4 Glacier on 6 December 2014 (Figure A4, Table A2). Therefore, MODIS results were inevitably lower than the results obtained with L5/TM and L8/OLI. This is also consistent with Roupioz et al. [53] who found that MODIS land surface reflectance was underestimated because of absence of subpixel terrain correction. The effect of subpixel terrain on a snow surface could be larger since snow is often located at a higher elevation within complex terrain, leading to underestimated MODIS snow albedo. Furthermore, subpixel terrain also increases the spatial variability of L8/OLI surface reflectance, therefore leading to a larger spatial variability of L8/OLI albedo than that of MODIS albedo. On the other hand, Liang [19] developed a single equation for the NTB conversion to be applied to MODIS surface in-band reflectance of several land surface types (>9), which may not be suitable because different land surface types have different reflectance spectra characteristics. It would be better to develop specific NTB conversion for ice surface by using ice reflectance data only.

5.3. Potential and Future Applications of the Updated Albedo Retrieval Method

The updated anisotropy correction method we present here to derive albedo from satellite data is promising for several future applications. On the one hand, it provides reasonable results for L8/OLI albedo retrieval, thus offering a more accurate data product to understand the role of albedo in glacier energy balance. On the other hand, this method can also be an option to retrieve albedo with other satellite data, such as Sentinel−2/MSI, AVHRR (Advanced Very High-Resolution Radiometer), ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), and WorldView−3 and −4 [15,19,28] as it covers the primary satellite shortwave bands for earth observation and is promising to retrieve broadband albedos over an extended selection of narrowband albedos. The differences in spectral coverage between the airborne BRDF data and the satellite sensors could result in uncertainty when directly applying our updated method to other satellite data. In the case of large differences in spectral coverage, the BRDF parameterization should be updated with better BRDF data on snow and ice if available. In addition, the NTB conversion should be updated to take into account differences in spectral coverage of different satellite sensors.

The two main advantages of our MODIS albedo product relative to MCD43A3 are: (1) fewer data gaps and improved accuracy of albedo retrievals, enabling retrieval of albedo throughout the entire year, (2) the potential to efficiently upscale data processing using GEE to cover large spatial domains. Since our MODIS albedo can be retrieved for any clear-sky day while MCD43A3 needs enough observations during a 16-day window, our product has fewer gaps during sporadically cloudy periods, such as December–March over the Zhadang Glacier (Figure 7c) and December–April over the Artesonraju Glacier (Figure 7g). Furthermore, our MODIS albedo generally performs better than the MCD43A3 product thanks to the updated BRDF method that is better suited. Indeed, since glacier surface materials frequently change between snow and ice during a 16-day window, the estimated BRDF of MCD43A3 is likely a combination of snow and ice BRDFs. As a result, during the ablation season the MCD43A3 product is unable to accurately retrieve snow albedo and capture albedo fluctuation. A more continuous data series and higher accuracy albedo means that our product can capture sudden albedo increases caused by snowfall better than the MCD43A3 product, for example on 16–17 June 2016 for Parlung No.4 Glacier (Figure 7d). These two advantages are beneficial for the determination of large-scale and long-term glacier albedo trends and can also help us better constrain and understand its seasonality.

6. Conclusions

In this study, we developed anisotropy correction models applicable to glacier snow and ice surface reflectance with high-quality airborne BRDF measurements. These models cover 10 spectral bands (339–2196 nm) for snow and six spectral bands (471–1271 nm) for ice and can estimate their anisotropic reflection factors well.

Albedo retrievals with the proposed anisotropy correction models were validated at eight glaciers across the main global glacier regions using L5/TM, L8/OLI and MODIS data. The mean bias, MAE and RMSE of the L8/OLI albedo retrieved using the Liang method and the updated anisotropy correction models were low, demonstrating a good agreement with field observations and better accuracy than a previously implemented method (the Knap method). However, the retrieved L5/TM albedo was overestimated because of frequent visible band saturation and underestimated by 0.1–0.5 when discarding a saturated band. Therefore, we recommend retrieving L8/OLI albedo with the Liang method [19] and the updated anisotropy correction models, while L5 albedo should be retrieved using the Knap method to avoid the saturation problem in the L5/TM visible band [18].

Furthermore, we developed a new method for MODIS albedo retrieval and validated it with the aggregated values of our L8/OLI and L5/TM albedo retrievals, and also compared it with the MCD43A3 product at the eight study glaciers. Both MODIS albedo products had lower values than L5/TM and L8/OLI albedo for most glaciers. Our MODIS albedo performed better than the MCD43A3 product in both absolute albedo estimation for six glaciers and glacier albedo temporal evolution for eight glaciers (fewer gaps and a better agreement with field observations). Our evaluation shows that our MODIS albedo product is best fitted to analyze long-term spatio-temporal variability of glacier albedo.

For future applications, because of its applicability to a wide range of narrowband sensors, the updated anisotropy correction can be easily applied to other satellite data and provide an opportunity to develop new, higher accuracy NTB albedo conversions. The updated retrieval method can also be easily applied to detect long-term albedo change in large-scale studies using big data processing platforms like Google Earth Engine.

Author Contributions

Conceptualization, S.R., L.J. and M.M.; methodology, S.R., E.S.M. and M.K.; software, S.R.; validation, S.R.; formal analysis, S.R., E.S.M., L.J., M.M., M.K., P.B., M.J.M. and T.E.S.; investigation, S.R. and E.S.M.; resources, L.J., M.M. and F.P.; data curation, S.R., E.S.M., T.E.S. and W.Y.; writing—original draft preparation, S.R.; writing—review and editing, S.R., E.S.M., L.J., M.M., M.K., P.B., M.J.M. and T.E.S. and F.P.; visualization, S.R.; supervision, L.J., M.M. and F.P.; project administration, L.J., M.M. and F.P.; funding acquisition, L.J., M.M. and F.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (grant no. 2019QZKK0103); the National Natural Science Foundation of China project (grant no. 91737205); the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDA19070102); the Chinese Academy of Sciences President’s International Fellowship Initiative (grant no. 2020VTA0001); the MOST High Level Foreign Expert program (grant no. GL20200161002); the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant no. 772751) and the Dragon 5 project (grant no. 32439).

Acknowledgments

Landsat Level-2 Surface Reflectance Science Product was courtesy of the U.S. Geological Survey. We are extremely grateful to the following individuals and groups for collecting and making available high-quality data of albedo by automated weather station records on glacier and station pictures, as this analysis would not have been possible without their contributions: Chuanxi Zhao and ITP-CAS (Parlung No. 4, Yang et al., 2011); Matt Nolan (McCall; https://arcticdata.io/catalog/view/doi%3A10.18739%2FA27S7HS5V, accessed on 18 February 2021); Wolfgang Gurgiser and the Department of Atmospheric and Cryospheric Sciences at the University of Innsbruck (Artesonraju Glacier, https://acinn-data.uibk.ac.at/, accessed on 15 January 2021); ICIMOD (Yala, ICIMOD. (2016). Meteorological data from Yala Base Camp automatic weather station [Data set]. ICIMOD. https://doi.org/10.26066/RDS.26859, accessed on 18 January 2021); IRD GLACIOCLIM Observation Service (Mera and Zongo, http://globalcryospherewatch.org/cryonet/sitepage.php?surveyid=23, accessed on 20 January 2021); Ekaterina Rets and coauthors (Djankuat, P Rets et al., 2019) and Shichang Kang (Zhang et al., 2016, 2013).

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A

Figure A1

(a–h) Pictures of the AWSs of our eight study glaciers. The green star in (h) represents the position of the AWS of Zongo Glacier.

[Figure omitted. See PDF]

Table A1

The band names and wavelength ranges (nm) of the airborne BRDF datasets used for anisotropy correction and corresponding three satellites’ data (L5/TM, L8/OLI and MODIS). ARCTAS data is used for the snow surface, and ARMCAS data for ice surfaces.

GlacierSurface	Airborne BRDF	L5/TM	L8/OLI	MODIS
Snow	339 (330–350)
	382 (370–390)
	480 (450–495)	band 1 (450–515)	band 2 (450–515)	band 3 (459–479)
	677 (650–720)	band 3 (630–690)	band 4 (630–680)	band 1 (620–672)
	873 (835–910)	band 4 (750–900)	band 5 (845–885)	band 2 (841–890)
	1032 (990–1075)
	1222 (1184–1258)			band 5 (1230–1250)
	1275 (1236–1319)
	1649 (1600–1709)	band 5 (1550–1750)	band 6 (1560–1660)
	2196 (2140–2260)	band 7 (2090–2350)	band 7 (2100–2300)	band 7 (2105–2155)
Ice	471 (462–482)	band 1 (450–515)	band 2 (450–515)	band 3 (459–479)
	675 (665–684)	band 3 (630–690)	band 4 (630–680)	band 1 (620–672)
	868 (858–877)	band 4 (750–900)	band 5 (845–885)	band 2 (841–890)
	1037 (1028–1047)
	1219 (1209–1229)			band 5 (1230–1250)
	1271 (1260–1281)
	560 (520–600)			band 4 (545–565)

Figure A2

Time series of 16-day moving average MODIS/Ren, MCD43A3 albedo product and field observations for the Parlung No.4 Glacier during 2013–2015.

[Figure omitted. See PDF]

Figure A3

Histograms of albedo difference between L8/OLI and AWSs (a,b), L5/TM and AWSs (c,d), MODIS and L8/OLI (e), L5/TM (f).

[Figure omitted. See PDF]

Figure A4

Comparison of surface reflectance in the visible bands (a: blue; b: green; c: red) of the 500 m L8/OLI aggregated (left), Terra/MODIS (middle) and Aqua/MODIS (right) over the Parlung No.4 Glacier on 6 December 2014 (same day as in Figure 6 in the main text). The white pixels in the Aqua/MODIS data represent cloud pixels.

[Figure omitted. See PDF]

Table A2

Range and difference of the 500 m L8/OLI aggregated, Terra/MODIS and Aqua/MODIS surface reflectance in the visible bands (blue, green and red) over the Parlung No.4 Glacier on 6 December 2014.

Band Name	Range			Difference
	L8/OLI	Terra/MODIS	Aqua/MODIS	Terra/MODIS-L8/OLI	Aqua/MODIS-L8/OLI
Blue	0.16–0.89	0.17–0.72	0.14–0.70	−0.09	−0.08
Green	0.19–0.90	0.18–0.77	0.16–0.73	−0.08	−0.08
Red	0.21–0.90	0.17–0.75	0.13–0.72	−0.1	−0.09

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Locations of the eight study glaciers with their AWSs. Glacier outlines are taken from the Randolph Glacier Inventory 6.0.

Figure 2. Broadband albedo from AWSs versus broadband albedo from overlapping L8/OLI pixels with the Knap and Liang methods for six of the study glaciers for which observations overlapped with L8/OLI images (Table 1) (a–f).

Figure 3. (a–d) Broadband albedo from AWSs versus broadband albedo from overlapping L5/TM pixels with the Knap and Liang methods for four of the study glaciers for which observations overlapped with L5/TM images (Table 1).

Figure 4. (a–g) Broadband albedos from the whole overlapping MODIS/Ren and MCD43A3 pixels versus L8/OLI broadband albedo for six study glaciers for which observations overlapped with L8/OLI images (Table 1). Results for Parlung No.4 Glacier were split into panels (c) and (d) for clarity given the large number of observations.

Figure 5. (a–e) Broadband albedos from the whole overlapping MODIS/Ren and MCD43A3 pixels versus L5/TM broadband albedo for the four study glaciers for which observations overlapped with L5/TM images (Table 1). Results for McCall Glacier were split into panels (a) and (b) for clarity given the large number of observations.

Figure 6. Comparison of albedo among 30 m L8/OLI (a), 500 m L8/OLI aggregated (b), MODIS/Ren (c) and MCD43A3 (d) albedo products over the Parlung No.4 glacier on 6 December 2014. More gaps were observed in the MCD43A3 product than in the MODIS/Ren product because of the absence of valid observations during 28 November–13 December 2014.

Figure 7. (a–h) Time series of MODIS/Ren, and the MCD43A3 albedo products and field observations for the eight study glaciers. For clarity, only albedo variability for two years is shown.

Figure 7. (a–h) Time series of MODIS/Ren, and the MCD43A3 albedo products and field observations for the eight study glaciers. For clarity, only albedo variability for two years is shown.

Table 2

Four BRDF parameterizations and $\tilde{f}$ models. BRDF parameterizations P1–P4 were derived separately by Reijmer et al. [22], Greuell and De Ruyter De Wildt [23], Warren et al. [45] and Knap and Reijmer [46]. The general form of BRDF parameterization is: $BRDF=a_0+a_1{g}_1+a_2{g}_2+a_3{g}_3$, where $a_i$ are weighting coefficients and $g_i$ are functions of the satellite view zenith angle after terrain correction (${\theta }_{vc}$) and the relative azimuth angle between the satellite and the sun ($\phi$).$c_1$, $c_2$, $c_3$ and ${\theta }_c$ are undetermined coefficients of $\tilde{f}$ and need to be estimated by numerically inverting BRDF data. Note that * indicates that the $\tilde{f}$ in P2 is from De Ruyter De Wildt [23].

Experiment	BRDF Parameterization			$\tilde{\mathit{f}}$
	${\mathit{g}}_1$	${\mathit{g}}_2$	${\mathit{g}}_3$
P1	${\theta }_{vc}{}^2$	${\theta }_{vc}{}^2\cos \phi$	${\theta }_{vc}{}^2{\cos }^2\phi$	$\left[c_1\left(g_1+\frac{1}{2}-\frac{{\pi }^2}{8}\right)+c_2{g}_2+c_3\left(g_3+\frac{1}{4}-\frac{{\pi }^2}{16}\right)\right]e^{\frac{{\theta }_{sc}}{{\theta }_c}}$
* P2	$\cos {\theta }_{vc}$	${\theta }_{vc}{}^2\cos \phi$	${\theta }_{vc}{}^2{\cos }^2\phi$	$\left[c_1\left(g_1-\frac{2}{3}\right)+c_2{g}_2+c_3\left(g_3+\frac{1}{4}-\frac{{\pi }^2}{16}\right)\right]e^{\frac{{\theta }_{sc}}{{\theta }_c}}$
P3	$(1-\cos {\theta }_{vc})$	$(1-\cos {\theta }_{vc})\left(-\cos \phi \right)$	$(1-\cos {\theta }_{vc})\left(2{\cos }^2\phi -1\right)$	$\left[c_1\left(g_1-\frac{1}{3}\right)+c_2{g}_2+c_3{g}_3\right]e^{\frac{{\theta }_{sc}}{{\theta }_c}}$
P4	${\sin }^2{\theta }_{vc}{\sin }^2\phi$	$\sin {\theta }_{vc}\cos \phi$	${\sin }^2{\theta }_{vc}{\cos }^2\phi$	$\left[c_1\left(g_1-\frac{1}{4}\right)+c_2{g}_2+c_3\left(g_3-\frac{1}{4}\right)\right]e^{\frac{{\theta }_{sc}}{{\theta }_c}}$

Table 3

Evaluation of the four anisotropy correction models for snow and ice with airborne BRDF datasets. The acronyms of the models are the same as in Table 2. In bold are the selected models used to retrieve albedo using L5/TM, L8/OLI and MODIS observations.

Glacier Surface	Parameterization Scheme	Calibration			Validation
		R	MD	Std	R	MD	Std
Snow	P1	0.935	0.0034	0.0736	0.939	0.0031	0.0726
	P2	0.935	0.0040	0.0742	0.938	0.0036	0.0732
	P3	0.931	0.0047	0.0768	0.935	0.0044	0.0757
	P4	0.883	0.0091	0.1076	0.887	0.0087	0.1064
Ice	P1	0.796	−0.0008	0.0807	0.815	0.0011	0.0796
	P2	0.798	−0.0009	0.0805	0.817	0.0012	0.0792
	P3	0.799	−0.0008	0.0803	0.817	0.0013	0.0791
	P4	0.794	−0.0003	0.0821	0.811	0.0022	0.0809

Table 4

Spectral band information and parameters of the anisotropy correction models used in this study ($\tilde{f}$) for each spectral band on the glacier surface (snow for P1 and ice for P2) are derived using airborne BRDF datasets. The weighting coefficient acronyms are the same as in Table 2. The meanings of R, MD and Std are the same as in Table 3, but for each band. Note that * corresponds to the weighting coefficients of ice at 560 nm taken from Greuell and De Ruyter De Wildt [23].

GlacierSurface	Airborne BRDF Dataset
	Central (Range)Wavelength (nm)	Weighting Coefficients				Calibration			Validation
		c1	c2	c3	${\mathit{\theta }}_{\mathit{c}}$	R	MD	Std	R	MD	Std
Snow	339 (330–350)	0.00514	0.00494	0.01585	1.57080	0.91	0.0001	0.014	0.93	0.0000	0.012
	382 (370–390)	0.00189	0.01029	0.02096	1.01490	0.91	0.0005	0.019	0.92	0.0005	0.018
	480 (450–495)	0.00000	0.00001	0.00002	0.12131	0.76	0.0078	0.188	0.77	0.0057	0.187
	677 (650–720)	0.00083	0.00384	0.00452	0.34527	0.95	0.0013	0.034	0.95	0.0010	0.033
	873 (835–910)	0.00123	0.00459	0.00521	0.34834	0.96	0.0015	0.038	0.96	0.0013	0.037
	1032 (990–1075)	0.00417	0.00709	0.00736	0.39306	0.97	0.0014	0.042	0.97	0.0010	0.041
	1222 (1184–1258)	0.00663	0.01081	0.01076	0.46132	0.98	0.0021	0.051	0.98	0.0017	0.050
	1275 (1236–1319)	0.00413	0.00954	0.01018	0.46048	0.97	0.0022	0.061	0.97	0.0019	0.059
	1649 (1600–1709)	0.00798	0.01744	0.01680	0.63119	0.96	0.0083	0.156	0.96	0.0091	0.163
	2196 (2140–2260)	0.00622	0.01410	0.01314	0.55261	0.97	0.0093	0.133	0.98	0.0087	0.125
Ice	471 (462–482)	−0.00369	0.00000	0.00007	0.27632	0.70	0.0119	0.045	0.72	0.0133	0.043
	675 (665–684)	−0.00054	0.00002	0.00001	0.17600	0.71	0.0075	0.053	0.74	0.0099	0.051
	868 (858–877)	−0.00924	0.00033	−0.00005	0.31750	0.82	0.0051	0.060	0.85	0.0073	0.057
	1037 (1028–1047)	−0.03533	0.00297	−0.00032	0.54050	0.87	0.0003	0.080	0.89	0.0024	0.077
	1219 (1209–1229)	−0.02388	0.00656	0.00227	0.58473	0.84	−0.0127	0.117	0.85	−0.0091	0.117
	1271 (1260–1281)	−0.02081	0.00683	0.00390	0.57552	0.84	−0.0176	0.128	0.84	−0.0168	0.131
	* 560 (520–600)	−0.02920	−0.00810	0.00462	0.52360	/	/	0.043	/	/	/

Table 5

Evaluation of L8/OLI and L5/TM albedo with AWS measurements. n is the number of observations. The NTB conversions by the Knap Method are Equations (5) and (6), the Liang Method is Equation (7). The BRDF parameterization schemes applied in combination with both methods are P1 for snow and P2 for ice, respectively.

Satellite	Glacier	n	FieldObservation	Knap Method				Liang Method
			Mean	Mean	Bias	MAE	RMSE	Mean	Bias	MAE	RMSE
L8/OLI	Djankuat	12	0.32	0.28	−0.04	0.07	0.09	0.31	−0.01	0.05	0.06
	Zhadang	15	0.75	0.73	−0.02	0.04	0.06	0.72	−0.03	0.04	0.06
	Parlung No.4	36	0.54	0.61	0.07	0.07	0.08	0.64	0.10	0.10	0.11
	Yala	28	0.67	0.69	0.01	0.06	0.07	0.69	0.01	0.05	0.06
	Mera	11	0.61	0.53	−0.08	0.08	0.08	0.61	0.00	0.03	0.03
	Zongo	18	0.42	0.48	0.06	0.10	0.13	0.49	0.06	0.10	0.11
	Average	/	0.55	0.55	0.00	0.07	0.09	0.58	0.02	0.06	0.07
L5/TM	McCall	10	0.46	0.42	−0.04	0.05	0.06	0.52	0.06	0.06	0.07
	Djankuat	15	0.20	0.28	0.09	0.11	0.14	0.33	0.13	0.14	0.20
	Artesonraju	13	0.29	0.39	0.11	0.11	0.17	0.47	0.19	0.19	0.25
	Zongo	24	0.36	0.45	0.09	0.09	0.10	0.50	0.14	0.14	0.17
	Average	/	0.33	0.39	0.06	0.09	0.12	0.46	0.13	0.13	0.17

Table 6

Evaluation of MODIS/Ren albedo and the MCD43A3 albedo against L5/TM and L8/OLI albedo. n refers to the total number of available MODIS pixels (and corresponding degraded L5/TM and L8/OLI albedos) across all scenes evaluated in Table 5. M_Landsat and M_MODIS are the mean albedos retrieved from Landsat and MODIS, respectively.

Satellite	Glacier	n	MODIS/Ren					n	MCD43A3
			MLandsat	MMODIS	Bias	MAE	RMSE		MLandsat	MMODIS	Bias	MAE	RMSE
L8/OLI	Djankuat	34	0.36	0.29	−0.07	0.08	0.09	34	0.36	0.29	−0.07	0.07	0.08
	Zhadang	132	0.71	0.57	−0.14	0.14	0.18	87	0.69	0.58	−0.11	0.11	0.13
	Parlung No.4	1313	0.66	0.55	−0.11	0.13	0.16	671	0.60	0.39	−0.21	0.22	0.24
	Yala	69	0.58	0.56	−0.02	0.09	0.11	105	0.55	0.49	−0.07	0.10	0.12
	Mera	112	0.56	0.47	−0.09	0.10	0.12	189	0.58	0.38	−0.20	0.20	0.22
	Zongo	37	0.58	0.41	−0.17	0.17	0.18	36	0.58	0.29	−0.29	0.29	0.30
	Average	/	0.57	0.47	−0.10	0.12	0.14	/	0.56	0.40	−0.16	0.16	0.18
L5/TM	McCall	589	0.50	0.40	−0.10	0.13	0.16	769	0.49	0.38	−0.12	0.15	0.19
	Djankuat	57	0.41	0.35	−0.06	0.07	0.09	48	0.38	0.32	−0.06	0.08	0.09
	Artesonraju	47	0.42	0.55	0.13	0.14	0.17	83	0.38	0.35	−0.02	0.07	0.09
	Zongo	45	0.50	0.38	−0.12	0.13	0.15	41	0.51	0.29	−0.22	0.22	0.26
	Average	/	0.46	0.42	−0.04	0.12	0.14	/	0.44	0.33	−0.11	0.13	0.16