{'fraction_non_alphanumeric': 0.03677731026959639, 'fraction_numerical': 0.02337540906965872, 'mean_word_length': 4.777940348902645, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anisotropy Parameterization Development and Evaluation for Glacier Surface Albedo Retrieval from Satellite Observations

: 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 efﬁcient means to retrieve glacier surface albedo over large and inaccessible areas and to study its variability. However, corrections of anisotropic reﬂectance of glacier surface have been established for speciﬁc 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 Reﬂectance 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 ﬁeld measurements collected on eight glaciers around the world. Our results show that the updated model: (1) can accurately estimate anisotropic factors of reﬂectance 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 ﬁeld 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.


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 highresolution 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.

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).

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 parameter-ization 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 3 ) 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, Figures 1 and 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].  [40]. This DSM was used to derive glacier surface slope and aspect for the topographic correction (Equations (3a) and (3b)).
Sens. 2021, 13, x FOR PEER REVIEW 4 igure 1. Locations of the eight study glaciers with their AWSs. Glacier outlines are taken from the Randolph Glacier ventory 6.0.

Datasets
Snow and ice BRDF data: High-quality BRDF measurements on mountain gl surfaces are rare, and most available measurements sample a limited spectral rang

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.

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 (α 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: where 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: The 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]. 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 Remote Sens. 2021, 13, 1714 7 of 28 parameterizations and calculated and validated f for each one for both glacier snow and ice surfaces, i.e., P1, P3 and P4 in Table 2.  [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 (θ vc ) and the relative azimuth angle between the satellite and the sun (ϕ). c 1 , c 2 , c 3 and θ c are undetermined coefficients of f and need to be estimated by numerically inverting BRDF data. Note that * indicates that the f in P2 is from De Ruyter De Wildt [23].

Experiment
BRDF Parameterizationf 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]: where θ v and θ vc are the view zenith angles of satellite before and after correction, θ s and θ sc are the sun zenith angles before and after correction, ϕ v and ϕ 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 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 θ c coefficients of f need to be estimated and validated with BRDF measurements for different spectral bands.
To fit and validate these f models, the measured f -values ( f m ) were directly calculated by airborne BRDF measurements (r BRDF ) and measured albedo (α m ) by Equation (4): 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 θ c by fitting 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 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.

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 (α Knap ): 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]: 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 (α Liang ): 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 (α MODIS−ice ): (8) Over snow surfaces, which show high reflectance, we used the conversion developed by Stroeve et al. [20], which showed better accuracy (α MODIS−snow ): (9) In Equations (8) and (9), b i (i = 1, 2, 3 . . . . . . ) represent narrowband albedo from the ith 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.

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 (MAE) and Root Mean Squared Error (RMSE).

Evaluation of the Anisotropy Corrections
To evaluate the anisotropy correction models ( Table 2) for snow and ice surfaces, we compared the f -values estimated ( f e ) according to Table 2, after having determined the c 1 , c 2 , c 3 and θ c coefficients by fitting the calibration BRDF data, with the 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 f e and 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 −3 , 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 f e and 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 nearinfrared 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. 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.  Table 4. Spectral band information and parameters of the anisotropy correction models used in this study ( 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 (Figures 2 and 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 (Figures 2 and 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. 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.

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 (Figures 4 and 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).   (Table 1).   (Table 1).   (Table 1). Results for McCall Glacier were split into panels (a) and (b) for clarity given the large number of observations. 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). 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.

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: 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 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].

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 2 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.

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.

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.  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). ARC-TAS data is used for the snow surface, and ARMCAS data for ice surfaces.   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 ARM-CAS data for ice surfaces.    Figure 6 in the main text). The white pixels in the Aqua/MODIS data represent cloud pixels.   Figure 6 in the main text). The white pixels in the Aqua/MODIS data represent cloud pixels.