Spatial and temporal variations of the surface albedo and other factors influencing Urumqi Glacier No. 1 in Tien Shan, China

Abstract

Glacier albedo controls the surface energy budget, the variability of which affects the glacier surface melt rate and, in turn, impacts the mass balance of the glacier. During 2013 and 2014, spatial and temporal variations of albedo were investigated using 18 Landsat images of Urumqi Glacier No. 1. Factors influencing these spatiotemporal profiles were analyzed. An established retrieval process, including geolocation, radiometric calibration, atmospheric, topographic, and anisotropic correction and narrow- to broadband conversion, was applied for the first time to Landsat-8 images. Differences between Landsat image derived albedo values and albedo values measured using a handheld spectroradiometer ranged from −0.024 to 0.049. Spatial and temporal variations of surface albedo were significant, especially in the ablation area. The variability of the values of ice albedo ranged from 0.06 to 0.44 due to topographic effects and light-absorbing impurities. The results suggest that this retrieval method can be used to investigate the spatial and temporal variability of surface albedo from Landsat-8 images on mountain glaciers. Moreover, as constant albedo values for ice and snow cannot be assumed, the distribution of albedo was not completely dependent on altitude under conditions of more intense ablation, and by reason of light-absorbing impurities during the melt season.

1. INTRODUCTION

Surface albedo is generally considered to be the broadband hemispherical averaged reflectance over the shortwave spectrum from 0.35 to 2.8 µm (Brock and others, 2000; Schaepman-Strub and others, 2006). It directly governs the net shortwave radiation flux at the surface of a glacier and largely determines the surface energy balance over the entire glacier (Dumont and others, 2011, 2012). Furthermore, net shortwave radiation is the primary energy source for the process of glacial melt (Van de Wal and others, 2005; Giesen and others, 2008; Six and others, 2009). Therefore, surface albedo plays an important role in controlling the rate of glacial melt. Glacier surface melt rates are highly sensitive to albedo (Oerlemans, 1992; Tedesco and others, 2011). Small changes in surface albedo can cause marked variations in the absorption of incoming solar radiation, and consequently influence the glacier surface melt rate, which in turn affects the mass balance, water availability and regional ecology (Winther and others, 2002; Dumont and others, 2014; Sun and other, 2014; Skiles and Painter, 2016).

Glacier albedo is determined by factors related to surface layer characteristics such as grain size, water content, impurity content, surface roughness, crystal orientation and structure, as well as by factors related to incident shortwave radiation, such as cloud cover and solar incidence angle (Kylling and others, 2000; Jonsell and others, 2003; Gardner and Sharp, 2010). Therefore, enormous spatial and temporal variations in surface albedo are observed, ranging from >0.9 on fresh dry snow to <0.1 on debris-rich ice surfaces (Cuffey and Paterson, 2010). Hence, as the crucial parameter in glacial energy-balance models and hydrological models, knowledge of albedo variations is vital for the accurate calculation of glacier surface melt by energy- and mass- balance models (Brock, 2004).

In general, an automatic weather station (AWS) mounted radiometer or portable spectrograph is used to obtain surface albedo in the field. (e.g., Oerlemans and Knap, 1998; Jonsell and others, 2003; Brock, 2004; Oerlemans and others, 2009; Jiang and others, 2010; Azzoni and others, 2016). Although field measurement provides reliable and high temporal resolution datasets, it is limited in its ability to reflect changes in albedo over the entire glacier for long periods. Moreover, the remoteness and hazardous nature of most glaciers pose a significant challenge to field data acquisition. Therefore, at the glacier-wide scale, spatiotemporal variations of surface albedo are normally extrapolated using modeled albedo parameterizations or data retrieved from satellites. Albedo parameterizations generally use ‘proxy’ variables, such as time and accumulated daily temperature, in order to account for the factors that physically control changes in the reflectance properties of the surface material, such as snow metamorphism and increases in the concentration of light-absorbing impurities in snow or on ice (Brock, 2004; Hock, 2005; Gardner and Sharp, 2010). However, the enormous spatial and temporal variability in glacier surface albedo, combined with the complex interplay of different factors determining it, lead to inaccurate estimations (Hock, 2005). This inherent shortcoming also becomes the major limiting factor for the accuracy of energy- and mass-balance models (Fitzgerald and others, 2012). To date, most studies have focused on the retrieval of surface albedo from satellite data to investigate the spatiotemporal variation of albedo. For instance, Landsat 5 Thematic Mapper (TM) or Landsat 7 Enhance Thematic Mapper Plus (ETM+) data, with a high spatial resolution of 30 m, have been used to research the spatial distribution of the surface albedo of a single glacier (Koelemeijer and others, 1993; Knap and others, 1999a; Klok and others, 2003; Paul and others, 2005; Fugazza and others, 2016). Moderate Resolution Imaging Spectroradiometer (MODIS) albedo products with a high temporal resolution have been used to monitor temporal variation in large-scale areas (Dumont and others, 2012; Brun and others, 2015; Moustafa and others, 2015; Wang and others, 2015; Wu and others, 2015). Other sensors, such as AVHRR, SPOT HRV, and AVIRIS, have also been selected to analyze the surface albedo, with a variety of purposes (Jacobsen and others, 1993; Stroeve, 2001; De Ruyter De Wildt and others, 2002; Dumont and others, 2011). Although the satellite products provide reasonable albedo estimates, some limiting factors include restriction to cloudless conditions and retrieval problems arising from topographic effects on reflected radiation.

Previous studies have suggested that a large-scale decline in albedo has appeared (Stroeve and others, 2013; Wang and others, 2014a; Ming and others, 2015). Glacier albedo variations occur on many timescales (e.g., seasonal, as the surface type changes from snow to ice and vice versa; short-term, due to melting of a thin snow layer and the accumulation of meltwater on the surface; within minutes, as a result of atmospheric conditions such as cloud cover and solar incidence angle) (Jonsell and others, 2003; Brock, 2004; Seo and others, 2016). Reduced albedo is linked to stronger warming trends in some regions, early melt onset, a lack of wintertime accumulation, expansion of bare ice area, high concentrations of impurities (cryoconite, dust and soot), melting of outcropped ice layers enriched with mineral content, and enhanced meltwater production and runoff (Moustafa and others, 2015).

Urumqi Glacier No. 1 located in the central region of Central Asia, is surrounded by vast areas of desert: the Taklimakan Desert and Qaidam basin desert to the south, the Gurbantunggut Desert to the north, the Gobi Desert to the east and the Peski Muyunkum and Peski Sary-Ishikotrau deserts to the west (Fig. 1a). As a reference glacier in the World Glacier Monitoring Service (WGMS) glacier monitoring network, it provides the longest and best observed glaciological and climatological monitoring record in China, since 1959. The spatial variation of surface albedo was quite obvious; the mean ranged from 0.12 to 0.81 during the summer ablation period (Bai and Yu, 1985). Furthermore, Li and others (2011) demonstrated that the reduction of glacial surface albedo due to the continued increase of the ablation area was one of the most important reasons for the rapid and accelerated shrinkage of Urumqi Glacier No. 1 over the past 30 years. However, previous studies primarily focused on monitoring changes based on several field points during the ablation period. For example, Takeuchi and Li (2008) measured the surface albedo values at six sites from 2 to 5 August 2006; the results suggested that the albedo was significantly reduced by surface dust. Ming and others (2016) conducted a similar study during the 2013 ablation period. It seems notable that albedo was not distributed over the entire glacier for a long period, and no research has been conducted using satellite data with high spatial resolution to investigate the spatiotemporal variation of surface albedo at the glacier-wide scale. Therefore, this study aimed to assess the spatiotemporal variation of surface albedo over entire glacier during one mass-balance year, and to analyze the impact factors of spatiotemporal profiles and feedback mechanisms of albedo, which will improve our comprehension of surface albedo over an entire glacier and provide references for energy- and mass-balance models of continental glaciers.

Fig. 1. Geographical location (a) and map (b) of Urumqi Glacier No. 1 in Tien Shan, China. Black circles mark sites of ground-based measurements taken with an Analytical Spectral Device on 18 August 2015 and 4 August 2016.

2. STUDY SITE AND DATA

2.1. Site description

Urumqi Glacier No. 1 (43°06′N, 86°49′E) is a northeast-facing valley glacier with two branches, located on the northern slope of Tianger Summit II (4484 m a.s.l.) in eastern Tien Shan, central Asia. The glacier comprises the East Branch (4267 m to 3743 m a.s.l.) and the West Branch (4484 to 3845 m a.s.l.), and covered an area of 1.646 km² in 2009 (Fig. 1b and Fig. 2) (Wang and others, 2014b). The glacier has been shrinking since 1959, and the shrinkage has accelerated since the 1980s (Li and others, 2011). Consequently, the glacier was separated into two branches in 1993. During the past 50 years, the glacier equilibrium line has averaged an altitude of ~4050 m a.s.l.

Fig. 2. Landscape photographs of the eastern (left) and western (right) branches of Urumqi Glacier No. 1. The pictures taken on 5 August 2016 (top) represent the bare ice surface during the ablation period, and the pictures taken on 1 September 2016 (bottom) represent the after-snow appearance.

The region is dominated by three dynamic elements: the westerly jet stream in the upper troposphere, the Siberian anticyclonic circulation and the cyclonic disturbances of the west wind circulation (Farinotti and others, 2015). According to records from the Daxigou Meteorological Station (3593 m a.s.l.), located 3 km southeast of Urumqi Glacier No. 1, the annual mean air temperature was ~−5.0°C, and the annual mean precipitation was 460 mm, from 1959 to 2015. The precipitation mainly originates from moisture carried by westerly winds; 78% of the annual total precipitation occurs from May to August (Fig. 3), predominantly as solid precipitation, such as snow, hail and sleet, at higher elevations (Yang and others, 1992). During winter, the temperature is controlled by the Siberian anticyclonic circulation (Mölg and others, 2014), causing very low temperatures and little precipitation. The long-term climate and glacier mass-balance observation data for 1959–2015 suggest that ablation of the glacier mainly occurs from May to August, with accumulation occurring simultaneously at higher elevations; a small amount of accumulation was observed during other months due to little snowfall and lower air temperatures.

Fig. 3. Mean air temperature and precipitation recorded by Daxigou Meteorological Station (located 3 km southeast of Urumqi Glacier No. 1 at 3539 m a.s.l.) from 1959 to 2015.

2.2. Data

In this research, we selected 18 Landsat-8 Operational Land Imager (OLI) and Landsat-7 ETM+ images with a spatial resolution of 30 m during cloudless days from June 2013 to October 2014 (http://glovis.usgs.gov/) to find the spatial and temporal variations of surface albedo on the Urumqi Glacier No. 1. Detailed information for these satellite images is listed in Table 1, including observation date, sensor type, sun Azimuth and sun elevation. In addition, photographs taken on 18 August 2015 and 4 August 2016 were used to estimate the accuracy of the surface albedo retrieved from satellite images.

Table 1. Landsat OLI and ETM+ data used to retrieve the surface albedo of the glacier; satellite images were downloaded from the US Geological Survey (USGS).

Furthermore, topographic correction and anisotropic correction demand accurate terrain geometry. Therefore, we used a topographic map with a scale of 1:50 000 and a contour interval of 20 m, measured using an aerial photograph from 1981. The topographic map was digitized using the Arc/Info geographic information system software. A DEM was created using a triangular irregular network data model and then transferred into gridded data (Li and others, 2002). The cell size of the DEM was 30 m, and the map projection used was Universal Transverse Mercator Projection (UTM) in World Geodetic System-1984 Coordinate System (WGS-84).

In order to estimate the accuracy of the surface albedo retrieved from satellite images, a field survey was conducted to acquire ice and snow spectral reflectance measurements with an Analytical Spectral Device (ASD) Field Spec Pro Spectroradiometer Handheld 2 in the East Branch (Fig. 1b) on 18 August 2015 and 4 August 2016 when the satellite passed overhead. The optical measurements followed the method used by Moustafa and others (2015). The instrument collects electromagnetic radiation ranging from visible to near-infrared (325–1075 nm) wavelengths with a spectral resolution of 3 nm and an instrument accuracy of <4%. In a standard top-of-atmosphere solar irradiance reference spectrum, the 325–1075 nm range comprises 81% of the total irradiance in the 300–3000 nm range. Therefore, the dominance of reflectance in the ASD visible and near-infrared wavelengths in determining broadband albedo can be used quantitatively to provide the spatiotemporal distribution (Moustafa and others, 2015).

The measured reflected radiation was converted into reflectance units, normalizing for incident radiation collected using a calibrated Lambertian Spectralon panel. Spectral data were collected in clear-sky conditions within 2 h of local solar noon. Spectral measurements were performed using a bare optical fiber with a 25° field of view. The distance between the optical fiber and the surface of glacier was 40 cm. The instrument was calibrated before measurement at each site using the surface of a white reference panel that is ~100% reflective and diffusive. Albedo measurements over a sloping surface must be made in a surface parallel plane to reduce geometrical errors. Three reflective spectra were recorded with 5 s intervals at each site to calculate the mean value. Broadband albedo values were calculated as a weighted average based on their spectral response curve and the amount of incoming solar radiation over the entire spectral range at each site (Wright and others, 2014; Moustafa and others, 2015).

To ensure a high-quality albedo dataset, quality was rigorously controlled. Apparent outliers were identified, defined as physically unrealistic spectral albedo values (>1.0) and raw spectra that were markedly different from other spectra across the entire spectral range (visible and near-infrared wavelengths) taken for the same sample. For 18 August 2015, four spectra were deemed outliers (total spectra collected = 24), and for 18 August 2016, six spectra were deemed outliers (total spectra collected = 48). The uncertainty of the field data was ±0.01. This was likely due to differences in incident shortwave radiation, such as cloud cover.

3. RETRIEVAL METHOD

We adopt the methodology in Klok and others (2003) to retrieve the glacier surface albedo from Landsat images. This is the only method that considers all important processes influencing the relationship between the satellite signal and surface albedo. Moreover, it has been widely applied for a number of mountain glaciers (e.g., Fugazza and others, 2016), and the good agreement between Landsat–derived albedo and measured albedo suggests that the retrieval methodology generally produces accurate estimates of glacier albedo. The steps include geolocation, radiometric calibration, atmospheric correction, topographic correction, anisotropic correction and narrow- to broadband conversion. However, there are a few differences with the methodology of Klok and others (2003), including atmospheric correction and topographic correction. These will be described in detail in the following sections.

3.1. Geolocation

We apply geolocation to each image by comparing the image to the topographic map. The resolution of the topographic map is greater than the pixel size of the satellite images, which is 30 m. First, we select the ground control points (GCP) on the topographic map. These GCPs must be easily recognizable, obvious and constant features such as a steep ridge, sharp mountain peak, or road junction (Guo and others, 2014). GCPs were focused around, but did not include, glacier margins owing to their annual variability. Then, we locate these feature points on most images within an accuracy of one pixel. Furthermore, when at least seven GCPs are obtained, a uniform distribution is achieved throughout the sub-images.

3.2. Radiometric calibration

In order to perform atmospheric corrections, the radiometric calibration must be computed differently for Landsat OLI and ETM+ data using the following equations.

For Landsat OLI data:

(1) $$L_i = \alpha DN + \beta. $$

Here, i is the Landsat OLI specific band, L is the radiance at the sensor's aperture [W/(m² μm sr)], DN is the satellite image digital number, β is the offset of the sensor, and α is the sensor gain [W/(m² μm sr)].

For Landsat ETM+ data:

(2) $$L = \displaystyle{{\left( {L_{\max} - L_{\min}} \right) \times \left( {DN - DN_{\min}} \right)} \over {\left( {DN_{\max} - DN_{\min}} \right)}} + L_{\min}. $$

Here, L _min and L _max are the minimum and maximum spectral radiance band limit, respectively. DN _max is the maximum satellite image digital number (255), and DN _min is the minimum satellite image digital number.

3.3. Atmospheric correction

Of the varied kinds of atmospheric correction models, 6S, MODTRAN, ATOCOR and LOWTRAN are widely applied. According to Zhang (2012), the error of snow albedo retrieval using FLAASH is less than that of 6S; hence, we applied FLAASH to perform atmospheric correction in the ENVI 5.3 software. FLAASH, based on the radiative transfer model MODTRAN4+, is a software package to support hyperspectral and multispectral sensors for shortwave infrared, visible and ultraviolet detection. FLAASH quantifies atmospheric properties, including the water vapor column, aerosol and cloud optical depths, as well as surface properties such as altitude and temperature. These parameters, required as inputs to FLAASH for atmospheric correction, are shown in Table 2, taking 20150818 OLI data as an example.

Table 2. Parameters required as input to FLAASH for atmospheric correction (using OLI data from 18 August 2015 as an example)

3.4. Topographic correction

In the discussion above, it was implicitly assumed that the surface was horizontal. While Urumqi Glacier No. 1 is relatively steep and surrounded by high mountains, we must take into account the topographic effects on the derived surface reflectance. In this research, we calculate the surface albedo with an improved topographic correction approach based on the C-factor correction. Compared with the C-factor correction, this method is proposed without linear regression. Experiments using a Landsat TM image and DEM in Three Gorges proved that the results of the improved method were better than those from C-factor and cosine correction (Huang and others, 2005). The improved expression is as follows:

(3) $$L_{\rm H} = (L_{\rm T} - L_{\min} ) \times \displaystyle{{[\cos \beta - \cos \alpha _{\min} ]} \over {[\cos \alpha - \cos \alpha _{\min} ]}} + L_{\min}. $$

Here, L _T is the reflectance of the inclined surface, L _H is the reflectance of the horizontal surface, L _min is the minimum reflectance of the satellite image, α is the local solar incident angle and β is the solar zenith angle.

3.5. Anisotropic correction

Snow and ice reflect solar radiation anisotropically, so the satellite-derived albedo depends on the view geometry. The function that describes the reflection of solar radiation as a function of the view geometry is called a bidirectional reflectance distribution (BRDF). It varies with wavelength, solar elevation and (sub) surface properties such as grain-size and impurity content. To account for this and calculate the surface band albedo, we divide the surface reflectance from Eqn (4) by the anisotropic reflection factor, f:

(4) $$a = \displaystyle{r \over {\,f\,(\theta _r, \theta _v, \varphi )}}.$$

Here, φ is the relative view azimuth angle, and θ _v is the satellite zenith angle with respect to the inclined surface. The anisotropic reflection factor differs for ice and snow owing to their different surface reflection patterns. We calculated f from BRDFs derived for ice by Greuell and de Ruyter de Wildt (1999) and for snow by Reijmer and others (2001), shown as Eqns (5) and (6).

For ice:

(5) $$f_{{\rm ice}} = a_0 + a_1 {\ast} \cos \theta + a_2 {\ast} \theta ^2 {\ast} \cos \varphi + a_3 {\ast} \theta ^2 {\ast} (\cos \varphi )^2. $$

For snow:

(6) $$f_{{\rm snow}} = b_0 + b_1 {\ast} \theta ^2 + b_2 {\ast} \theta ^2 {\ast} \cos \varphi + b_3 {\ast} \theta ^2 {\ast} (\cos \varphi )^2. $$

In the two equations above, f is the anisotropic reflection factor and the variables a _i and b _i are coefficients used in two previous studies (Greuell and de Ruyter de Wildt, 1999; Reijmer and others, 2001). The view zenith angle and relative view azimuth angle are denoted by θ and φ, respectively.

3.6. Narrow- to broadband conversion

Finally, we calculated the broadband albedo from the narrow albedo in bands 3 and 5 for OLI, and bands 2 and 4 for ETM+, following the parameterization in Knap and others (1999b):

(7) $$\eqalign{\alpha &= 0.726{\ast}\alpha _{2{\rm or}3} - 0.322{\ast}\alpha _{2{\rm or}3} ^2 - 0.051{\ast}\alpha _{4{\rm or}5}\cr &\quad + 0.581{\ast}\alpha _{4{\rm or}5} ^2.}$$

Satellite images in band 2 easily saturate, because fresh snow has high reflectivity at visible wavelengths. When this occurs, another formula is used:

(8) $$\alpha ^{\prime} = 0.782{\ast}\alpha _{4{\rm or}5} + 0.148{\ast}\alpha _{4{\rm or}5} ^2. $$

This yields an accuracy with residual standard deviations of only 0.009 and 0.014.

As mentioned, the BRDF parameterizations differ for ice and snow, but whether an area is ice or snow is not known beforehand. We solved this problem following the method used by Klok and others (2003). Specifically, we calculated broadband albedo with both BRDF parameterizations. Then, we assumed that areas having a calculated surface albedo <0.5 with both BRDF parameterizations represented ice and that pixels having a calculated surface albedo >0.5 with both BRDF parameterizations pertained to snow. We omitted the remaining pixels unless the difference between the two values did not exceed 0.1, in which case we used the mean value.

4. RESULTS

4.1. Comparison of satellite-derived albedo with ground measurements

All processing steps introduce some uncertainty in the surface albedo. However, the uncertainties resulting from most processing steps are small and contribute to providing a better explanation of the spatial and temporal patterns shown by the satellite-derived albedo of the glacier. Hence, to validate the accuracy of the above methodology, we compared the albedo measured from ASD with albedo retrieved from satellite images at the sample site pixels, and a scatterplot of satellite-derived albedo against ground-based albedo for the time when the satellite passed over is shown in Figure 4. Despite different spatial resolutions, a good agreement can be found, with differences from −0.024 to 0.049 at all sample sites; the average absolute error is 0.022, and the R ² is 0.994, which is comparable with that determined by Klok and others (2003) and admissible within the 0.05 error ranges of the climate model used (Wang and Gao, 2004).

Fig. 4. Comparison between the albedo retrieved from Landsat data and the albedo measured by an Analytical Spectral Device on 18 August 2015 and 4 August 2016.

Next, we analyzed the possible error sources, which include: (1) the error may be a result of the atmospheric correction process. To account for the lack of actual atmospheric parameters when the Landsat images were acquired, we used a standard atmospheric model proposed by FLAASH code (the middle latitude summer/winter model) to estimate atmospheric effects. In this regard, there is a slight disconformity with reality; (2) the error is probably due to a contrast of spatial resolution between the ground-based and satellite-derived data. Surface albedo can have high spatial variability over distances of only a few meters while the area covered by ground measurements is relatively small, with an area of ~1 m × 1 m, compared with the OLI pixel size. Thus, it is very probable that any particular sample point had a significantly lower or higher albedo than the surrounding 30 m × 30 m area (Knap and others, 1999a); (3) the error is likely introduced by anisotropic correction of ice and snow. The ranges of BRDF parameterizations of ice and snow are 26.4–74.2° and 15.91–65.51°, respectively (Greuell and de Ruyter de Wildt, 1999; Reijmer and others, 2001). The application of the BRDF parameterizations to albedo may be outside these ranges in some areas. Besides, we also applied the BRDF parameterization of ice to zones covered with debris, for which the BRDF parameterization actually does not apply; (4) the topographic map was measured in 1981 but satellite images were obtained in 2015. In this study, topographic correction and anisotropic correction depend largely on the topography calculated for each grid cell. Therefore, the time lag between the topographic map and the satellite images is also an important factor that must be considered when discussing error. As the results of Wang and others (2014b) show, the average thinning of Urumqi Glacier No. 1 was ~8.9 m during the period 1981–2009. Hence, the surface terrain of glacier may vary slightly over time.

4.2. Spatial distribution of surface albedo on Urumqi Glacier No. 1

In order to investigate the spatial variation of the satellite-derived albedo over an entire glacier in one mass-balance year, we selected six scenes of Landsat images from September 2013 to August 2014. According to the retrieval method described in Section 3, the albedo images are shown in Figure 5.

Fig. 5. Spatial distribution of Landsat-derived surface albedo for Urumqi Glacier No. 1.

On 13 September 2013, the surface albedo calculated from the image was high, with most values (69%) between 0.5 and 0.9. Spatial variation of the albedo appeared slightly as a fluctuation along the horizontal profile, such as at an elevation of ~4050 m a.s.l. for the west branch and 3900 m a.s.l. for the east branch. This was likely caused by a discrepancy in surface ablation intensity due to heat conduction from ambient bare rocks. On 24 March 2014, the mean surface albedo was 0.57–0.64, with ~63% of values between 0.5 and 0.8, indicating that the glacier surface was mainly old snow, based on the results reported as 0.63 (0.51–0.79) by Bai and Tetsuo (1989) in western China. However, spatial changes in surface albedo were not observed during this period. Some low albedo values resulting from terrain shielding were observed at high altitudes; hence, they were removed from our study.

In the early ablation season (12 June 2014), the derived surface albedo of the glacier was between 0.7 and 0.9 (82%), and the mean was 0.76. The spatial distribution of surface albedo was still not obvious. Higher surface albedo meant that fresh snow covered the entire glacier surface. In the middle ablation season (15 August 2014), the surface albedo was reduced to 0.48 and increased obviously with altitude in the lower part of glacier. As a consequence of the strong optical contrast between bare ice and snow-covered areas, a sharp transition of the surface albedo was visible, representing the position of the snow line. The altitude of the snow line was ~3950 m a.s.l. in the eastern branch and 4050 m a.s.l. in the western branch. The increase of surface albedo at the snow line was ~0.2 and occurred over a short distance. The surface albedo was higher than 0.6 above the snow line, while it varied between 0.2 and 0.5 below the snow line. Above the snow line, the snow albedo increased weakly with elevation. From August 15 to 31 2014, the snow line retreated and eventually reached the upper part of the glacier. The lowest surface albedo was also observed ~0.34. The spatial variation in surface albedo remained remarkably clear, characterized by an increase with elevation. The albedo images of 31 August 2014 and 5 September 2013 demonstrated how small the accumulation area had become by the late ablation season.

These findings suggest that the spatial distribution of surface albedo was different with respect to the transition in dominant surface type from snow to bare ice. When the entire glacier was covered by snow, the spatial distribution of surface albedo was not visible despite the presence of both fresh and aged snow; this has also been observed by Wang and others (2015) in Dongkemadi Glacier, central Tibetan Plateau. Alternatively, when the glacier surface was dominated by both snow and bare ice, partially overlain with sediment and debris, a trend of surface albedo increasing with altitude was observed at the ablation zones, was more rapid near the snow line, and fluctuated in the accumulation zones. This spatial distribution of surface albedo was also demonstrated by Wang and others (2014a) for three glaciers in western China and by Koelemeijer and others (1993) in Hintereisferner, Austria. However, it was very interesting that the surface albedo was higher on the tongue of the glacier on 5 September 2013 and 31 August 2014, as also reported by Brock and others (2000) for Haut Glacier d'Arolla, Switzerland, Greuell and others (1997) for Pasterzenkees, Austria, and Klok and others (2003) for Morteratschgletscher, Switzerland. The results indicate that the ice albedo did not show a dependence on altitude in localized areas. This will be further discussed in detail in the following section.

4.3. Temporal variation of the surface albedo of the glacial ablation area

Table 3 shows mean surface albedo values derived from satellite images for each date. The mean was very high from late April to June (>0.7), with a general trend of decreasing albedo as the melt seasons proceeded, and the lowest values was derived in the late ablation season (~0.35). Sudden pronounced increases in albedo were observed during melt seasons due to summer snowfalls. For example, between August 7 and 15, the mean albedo increased by 0.10 owing to a summer snowfall on August 14. Moreover, the values of free ice albedo (defined as the area below 4050 m a.s.l. according to the result of the mass-balance model) also showed significant changes, with a wide range from 0.06 to 0.44. Temporal variation of surface albedo was dominated by the transition of snow to ice and vice versa.

Table 3. Mean albedo derived from 18 satellite images from 2013 to 2014

The percentages (n) indicate the fraction of glacier pixels with values from 0.06 to 0.98 according to Cuffey and Paterson (2010), where ‘σ’ is the standard deviation. ‘*’ indicate the fraction of glacier pixels with values from 0.30 to 0.98. Here, we considered that values <0.3 were affected by shading.

5. DISCUSSION

5.1. The suitability of the derived method to extract albedo from Landsat-8 data

In this study, although the surface albedo values were successfully derived from Landsat-8 satellite images, based on the method proposed by Klok and others (2003), and the results were within an admissible range, it is notable that the retrieved procedure was proposed for Landsat ETM+. Using same protocol to derive albedo from Landsat-8 satellite images was based mainly on the following considerations.

Landsat-8, launched on 11 February 2013 carrying an OLI instrument, is intended to be broadly compatible with the previous Landsat-7 ETM+ instrument, although they have slightly different relative wavelength ranges (Flood, 2014). The equivalent bands in the near infrared wavelengths, in particular, have been narrowed to avoid some atmospheric attenuation features; due to the addition of an extra band at the blue end of the spectrum on the OLI instrument, the band numberings have changed between the instruments. However, this study used only bands 3 and 5 of the OLI instrument and only in the process of narrow- to broadband conversion; the between-sensor difference in these two equivalent bands is subtle, especially in band 3. Thus, we argued that the difference in relative wavelength range could be neglected. The other possible difference is from the BRDF correction; the measurements were made in TM band 2 and 4 for ice by Greuell and de Ruyter de Wildt (1999) and for snow by Reijmer and others (2001), but De Ruyter De Wildt and others (2002) argued that they were also applicable to AVHRR bands 1 (550–680 nm) and 2 (725–1100 nm), respectively. Hence, we also applied the same protocol to OLI bands 3 and 5. Other processing methods were proposed based on the radiation transmission model and were suitable for different sensors.

Moreover, the study of ETM+ and OLI reflectance in the Australian landscape suggests that the systematic component of the between-sensor difference is much smaller in the surface reflectance than in top-of-atmosphere reflectance (Flood, 2014). In the case of top-of-atmosphere reflectance, the average difference can be as much as 6% of the reflectance, while the average difference in surface reflectance can as low as 2%. This is because some of the difference is due to atmospheric effects interacting with the change in spectral response of the instrument, and the atmospheric correction is able to account for these differences directly, by proper modeling of radiative transfer though the atmosphere (Flood, 2014). It should be emphasized that all of the study of Flood (2014) was carried out in land pixels and does not include any pixels with snow and ice; thus, the between-sensor differences may depend on the nature of the target. Detailed differences for snow and ice pixels should be researched further.

5.2. Influence of terrain factors on spatial variation of the glacier surface albedo

Spatial variation in the surface albedo on a single glacier is primarily controlled by the glacier surface properties, which are affected by topographic factors such as altitude, slope and aspect. The glacier was divided into several altitude zones at 50 m intervals using the topographic map. Figure 6 demonstrates that the satellite-derived glacier albedo for each zone along two center lines generally increases with elevation during the late ablation season (5 September 2013 and 31 August 2014). A steep gradient is visible in the center of the glacier, which also represents the position of the snow line. However, the ice albedo in the lower parts of the glacier decreased weakly with elevation. This was probably due to an increase in the amount of supraglacial debris. The snow albedo in the highest section of the glacier decreased obviously. The standard deviation is rather large, showing values of 0.12 to 0.22 especially in the east branch, indicating that the physical composition of the glacier surface is variable. In the west branch, parts of glacier surface may not be covered by snow; mixed pixels in the higher altitude zones reduce the mean albedo. Furthermore, the surface of the glacier is very crevassed at the summit of the east branch; a glacier meltwater pool was discovered here (Li, 2005), possibly resulting in a decrease in surface albedo.

Fig. 6. Mean albedos along the center line of Urumqi Glacier No. 1, as derived from Landsat images as a function of altitude. The error bars show the standard deviation of the distribution. The location of the center line is shown in Figure 1.

In addition, the influence of the solar incident angles, which are caused by variation in slope and aspect, should also be considered. Figure 7 shows the relationships between surface albedo and solar incident angles in the ablation zone based on the assumption that the surface condition is uniform. The surface albedo obviously increases with the solar incident angle. This could be explained as follows: the larger solar incident angle enhances the chance that photons are scattered by snow particles of the surface layer, especially for the near infrared waveband; thus, the proportion of energy absorbed by surface snow decreases, which results in an increased surface albedo.

Fig. 7. Albedo variation with solar incident angle at the ablation zone based on the Landsat image derived albedo values and co-instantaneous solar angles for each pixel.

5.3. Influence of air temperature and precipitation on temporal variation of the glacier surface albedo

Although variation in solar incident angles can influence the spatiotemporal distribution of surface albedo, the deposition, change, and removal of snow and exposure of bare ice, are primary drivers of the temporal changes in surface albedo. Spatial and temporal distribution of snow and bare ice cover are controlled by air temperature and precipitation. To investigate the influence of air temperature and precipitation on surface albedo, monthly variations of air temperature, precipitation and surface albedo in 2013/14 are analyzed (Fig. 8).

Fig. 8. Monthly trends of the temperature, precipitation and glacier albedo on Urumqi Glacier No. 1. Air temperature and precipitation are from the Daxigou Meteorological Station (43°06′N, 86°50′E), located 3 km southeast of Urumqi Glacier No. 1 at 3539 m a.s.l.

As the air temperature rises, a positive feedback mechanism begins, involving exposure of bare ice, increased meltwater and impurities, accelerated snow metamorphism and reduction of surface albedo. Furthermore, glacier albedo is highly sensitive to snowfall; once a snowfall event occurs, albedo values rapidly increase. However, the effect of air temperature and precipitation on albedo varies during different ablation stages. From March to June, surface albedo was mainly affected by precipitation; frequent snowfall events caused an increasing trend. Although air temperature rose simultaneously, it did not compromise the impact of frequent snowfall events. Between June and August, however, air temperature was the major controller of changes in surface albedo. Surface albedo reduced as air temperature rose because solid precipitation occurred predominantly at higher elevations.

5.4. Relationship between light-absorbing impurities and surface albedo

In addition to shifts in dominant surface type from snow to bare ice that cause variations in glacier albedo, one of the most important factors related to surface characteristics driving albedo changes are light-absorbing impurities (LAIs). Recently, various research groups have revealed that LAIs including brown carbon (BrC), black carbon (BC), mineral dust (MD) and living organisms can reduce glacier albedo, thus accelerating glacier melt. Although the contribution of LAIs to albedo reduction varies in different regions (Qu and others, 2014; Li and others, 2017; Niu and others, 2017; Schmale and others, 2017), there is now a general consensus that both BC and mineral dust are important contributors. Urumqi Glacier No. 1, located in an arid to semi-arid region of central Asia, is surrounded by vast desert areas; the bare ice surface is mostly covered by brown dust during ablation season (Dong and others, 2009). Takeuchi and Li (2008) first revealed that dust deposited on the glacier surface contained high levels of organic matter, including living cyanobacteria. Moreover, the spectral albedo of the glacier surface showed spectrum curves typical of those of snow and ice contaminated with dust, indicating that the albedo was significantly reduced by the surface dust on the glacier. In August 2013, Ming and others (2016) collected fresh and aged snow samples on the east branch of Urumqi Glacier No. 1, and their simulation showed that BC and dust were responsible for 25 and 7% of the albedo reduction in the accumulation zone, respectively.

In order to further investigate the effect of LAIs on albedo, the spectrally resolved reflectivity of snow and ice surfaces with different impurity coverages (Fig. 9) is shown in Figure 10. Surface albedo was characterized by a large range, whether on snow or ice surface. At the ice sites, the albedo values varied from a maximum of 0.31 at clean-ice site III to a minimum of 0.06 at debris-rich-ice site VI; other values were 0.25 for slightly dirty-ice site IV and 0.14 for dirty-ice site V. Moreover, the ice-surface spectrum (IV, V and VI) exhibit the typical spectral curves of a dirty glacier ice surface, namely, the spectral reflectance in the visible wavelength range (380–750 nm) are generally low compared with those of clean bare ice surface (III) but slightly higher in the range of 550–750 nm (Zeng and others, 1983). The snow-surface albedo values are 0.74 for clean-snow site I and 0.57 for dirty-snow site II. The spectral curve (II) shows that reflectivity increases as wavelength increases from 380 to 600 nm, holds steady between 0.6 and 0.65 from 600 to 800 nm, and then decreases from 800 to 950 nm, which agrees with the typical snow spectral curves due to dust contamination (Zeng and others, 1983). Besides, the variation of surface reflectivity was highest in the visible range (380–750 nm) and only minor in the near-infrared (750–950 nm) at both the ice and snow sites. These results suggest that impurities significantly reduced the surface albedo of Urumqi Glacier No. 1, and the influence of impurities on albedo may largely be due to reflectivity changes in the visible range since sediment and impurities strongly absorb in the visible spectrum.

Fig. 9. Close-up images of the surface of Urumqi Glacier No. 1 and measured albedo (α) values. I–VI represent the clean snow, dirty snow, clean ice, slightly dirty ice, dirty ice and debris-covered ice, respectively. These photographs were taken on 30 July 2016 at sites I, J, B2, H1, G2 and E1 shown in Figure 1b.

Fig. 10. The spectral reflectance at the clean snow (I), dirty snow (II), clean ice (III), slightly dirty ice (IV), dirty ice (V) and debris-covered ice (VI) sites under clear-sky days. The conditions of the measurement surface are shown in Figure 9.

5.5. The effect of water to surface albedo

Water (derived from both melt and rainfall) plays a significant role in driving ice albedo changes; it washes out the finer sediments on the glacier (Oerlemans and others, 2009) and makes the ice surface smooth. The surface albedo was higher on the lower part of the glacier than on the middle part of the glacier on 5 September 2013 and 31 August 2014 (Fig. 3). This profile of spatial distribution contradicted the general assumption that the surface albedo increases with altitude. Moreover, from 15 to 31 August 2014, the area of exposed glacier ice increased; this indicated that the glacier experienced strong ablation. It was likely to cause the accumulation of impurities melted out from internal ice layers. As a consequence, the albedo decreased from >0.4 to <0.2 in the middle part of the glacier (~3900–3950 m a.s.l. in the east branch and 4000–4050 m a.s.l.in the west branch). However, in the lower part of the glacier (~<3900 m a.s.l. in the east branch and <4000 m a.s.l.in the west branch), the albedo still held steady between 0.2 and 0.4. A possible explanation for these higher albedo values on ice is related to removal of surface dust and debris by meltwater and rainfall, as a result of high air temperatures (>4°C) and steep slopes (~15–30°) in the lower part of glacier (Fig. 11). Brock (2004) proposed that the surface ‘weathering crust’ and cryoconite holes, within which small particles of rock and fine debris are normally incorporated, break down to leave a smooth ice surface from which debris may be easily washed by runoff. In addition, the occurrence of the rainfall washing out effect was found to be short lived (~1.8 days, ranging from 1 to 4 days) by Azzoni and others (2016). However, there was 0 mm of precipitation 5 days before 31 August 2014, according to records from Daxigou Meteorological Station. Therefore, the possibility of rainfall washing can be excluded; we can deduce that higher ice albedo in the lower part of the glacier is probably a result of meltwater washing, due to high air temperatures. This can be confirmed by the emergence of surface streams. Nevertheless, Azzoni and others (2016) also proposed that meltwater could decrease the surface albedo owing to its lower reflectivity (i.e. equal to 0.05–0.10; Hartmann, 1994). The effect of meltwater on ice albedo may be related to the amount of water involved; if it is not enough to yield surface run-off, the ice albedo may decrease. However, this process was not observed explicitly; further measurements of glacial albedo are required to quantitatively relate albedo changes and water.

Fig. 11. Slope of Urumqi Glacier No. 1 averaged over 50 m height intervals.

6. CONCLUSIONS

In this study, we derived surface albedo by applying, for the first time, the approach proposed by Klok and others (2003) to Landsat-8 images. This retrieval method took all possible sources of influence on the satellite signal into account when estimating albedo. The process includes geolocation, radiometric calibration, atmospheric correction, topographic correction, anisotropic correction, and narrow- to broadband conversion. Although some discrepancies were found between Landsat-derived albedo and albedo measured by ASD (within the range −0.024 to 0.049), the results were reliable for investigating the spatial and temporal variability of surface albedo of a mountain glacier. These deviations might have been caused by the different spatial resolutions between Landsat-derived albedo and albedo measured by ASD, employment of a standard atmospheric model instead of the actual atmospheric profiles when the Landsat images were acquired, and the time lag between the topographic map and satellite images.

We also found that the variability of the values of free ice albedo were significant, ranging from 0.06 to 0.44 due to the effects of topography including slope and aspect, and LAIs. These results emphasize that albedo values for ice and snow cannot be assumed to be constant. The distribution of albedo may not be completely dependent on altitude under conditions of intense ablation, and other factors, such as LAIs, may be active during melt season.

Although we successfully researched the spatiotemporal variability of surface albedo in one mass-balance year, we could not obtain images from November 2013 to February 2014 due to the influence of cloud and terrain shadow. Furthermore, due to the low satellite return period of 16 days, daily or hourly are still unclear and further research is needed to understand these variations of surface albedo over the entire glacier.