Firn line altitude is an indicator of equilibrium line altitude (ELA) and plays an important role in predicting climate behavior and global change. Landsat 30-m resolution data from clear sky at the end of 21 summer melt seasons (July and August, 1990–2011) were used to obtain the firn lines and extract the glacier boundary at Qiyi Glacier in the Qilian Mountains of China. High-resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer 30-m resolution Global Digital Elevation Model (ASTER GDEM) was employed to obtain the firn line altitudes. The firn line altitudes on the Qiyi Glacier increased from 4540 m a.s.l. to 5000 m a.s.l. during the study period, a rise of 460 m. Changes in firn line altitudes had good positive correlations with air temperature data and were in good agreement with measured and simulated ELA changes. The calculated firn zone area decreased (increased) with the increase (decrease) of firn line altitude. Linear regression was used to fit the variational trend, indicating that the firn zone area decreased by ∼1.19 km2, with a change rate 0.057 km2 a-1, during the study period. The firn zone area ratio decreased from 76.2% to 18.1% during the same time.
Introduction
The snowline is the hub of glacier development (Shi et al., 1988), which divides glacier surfaces into ice and snow/firn surfaces and delineates the minimum snow-covered areas (Seidel et al., 1997). In the summer months, ice crystals in glaciers deform quickly due to meltwater percolation through them. As a result, some of the snow becomes firn (granular, partially consolidated snow that has passed through one summer melt season but is not yet glacial ice) at the end of a summer on a glacier. The lower limit of accumulated firn is called the firn line. A firn line maps the firn area and separates the firn from glacier ice on a glacier surface at the end of summer melt season, and delineates the minimum firn-covered area. The equilibrium line altitude (ELA) is the altitude at which the annual mass balance equals zero in a hydrological year (Cogley et al., 2011). When there is no superimposed ice on a glacier, the snowline altitude will be coincident with the ELA at the end of the melt season. In equilibrium and steady state, the firn line altitude can conform with the ELA in the absence of superimposed ice. The firn line altitude will be higher than the ELA in a year of negative mass balance, and it will be lower than the ELA in a year of positive mass balance. Air temperature and precipitation play important roles in snowline altitude, firn line altitude, and ELA change, all of which rise with temperature increase and descend with precipitation increase.
The firn line altitude is a crucial index for evaluating future glacial changes (Kaur et al., 2009). The distribution of firn line altitudes also depends on the local temperature, precipitation, topography, and other factors that increase the complexity of the spatial distribution of the firn line (Kerr and Sugden, 1994). Firn line depression can theoretically relate to climatic perturbations (Kuhn, 1989; Ohmura et al., 1992). For large-scale areas and long-term periods, an altitude shift of the firn line is a response to certain climate change behavior and can indicate the hydrologic balance (Droz and Wunderle, 2002). Firn line altitudes are also an appropriate proxy for ELA and, consequently, for mass balance and climate reconstructions (McFadden et al., 2011).
Variation in firn line altitude is crucial in a local area because glaciers are important water sources for agriculture and industrial production, especially during dry seasons (Kaser et al., 2003; Mark, 2005; Bradley et al., 2006). Mass balance observations of glaciers could provide valuable information about climate conditions and possible changes (Braithwaite and Zhang, 1999). The firn line altitude is not influenced by yearly variations in ELA, but variation in the average ELA will eventually lead to firn line altitude change (König et al., 2000).
Glacier monitoring has been implemented on only a small number of glaciers because field observations are highly labor-intensive (Heiskanen et al., 2003). Firn lines were already obtained from microwave and optical remote sensing technology, the latter of which mainly includes statistical and classification methodology.
The methodology of identifying snowlines using remote sensing images in mass balance assessment was first presented by Ostrem (1975). Hall et al. (1989) and Williams et al. (1991) found that snowlines were easily identified in satellite imagery during the melt season. Seidel et al. (1997) defined “snowline” as a belt of approximately 50% snow coverage, and they used Landsat TM (Thematic Mapper) and SPOT ( Système Probatoire d'Observation de la Terre) data to evaluate its statistical objectivity. De Angelis et al. (2007) calculated the TM4/TM5, TM2/TM5, and TM4/TM7 band ratios to classify glaciers, which is more accurate than single band satellite imagery because the method reduces the effects of the non-Lambertian snow reflectance and shadows. Still another method to determine the snowline is through the NDSI (Normalized Differ ence Snow Index) and classifying the reflectance map (gained from preprocessing the image) by supervised classification, unsupervised classification, or decision-tree methods (Heiskanen et al., 2003). Dyurgerov (1996) estimated transient accumulation area ratios using observed snowlines and calculated transient mass balance values.
In this study, we used Geographic Information System (GIS) technology and remote sensing satellite images acquired at the end of the summer melt seasons (July and August) at Qiyi Glacier from 1990 to 2011. We preprocessed them, including by radiometric calibration and atmosphere correction, and then we produced reflectance maps. Narrowband-to-broadband albedo conversion was carried out to gain the “real albedo” from all the solar radiation. All of this, combined with preregistered ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model), enabled us to determine the firn line altitudes and firn zone areas in the Qiyi Glacier, elucidate the variations thereof, and discuss the influence of temperature and precipitation on firn line altitude. Finally, based on the calculated firn zone area and glacier area, we investigated changes in the relative firn zone area ratio.
Study Area
The Qiyi Glacier (39°14.22′N, 97°45.34′E, numbered CN5Y437C18) is located in the Qilian Mountains in the north edge of the Qinghai-Tibet Plateau (Fig. 1). The meltwater of the Qiyi Glacier flows into a branch of the Beida River. This glacier was named in memory of the Chinese glaciologists who first climbed it on 1 July 1958 (Snow and Ice Research Team of Chinese Academy of Sciences, 1958), and the scientific study of Chinese glaciers began with this glacier in the same year. By its morphology, the Qiyi Glacier is a cirque-valley glacier, and according to its physical properties, it is categorized as a subcontinental glacier. From ground stereophotogrammetry in 1975, the Qiyi Glacier was 2.817 km2 in area and 3.8 km long, and its summit altitude and terminus altitude were 5158.8 m a.s.l. and 4304 m a.s.l., respectively.
Methods and Data
IMAGE SELECTION
Images were preferentially selected at the end of the summer melt seasons (July and August) to maximize the probability of extracting firn lines at the times of minimum firn covers. We had to take into account the effect of cloud, so we excluded the Landsat images with mean cloud coverage greater than 20%. Landsat TM/ETM+ images with 30-m spatial resolution were used for analysis, which include 10 sensor recordings in 1990, 1992, 1994, 1995, 2001, 2003, 2006, 2007, 2009, and 2011 (Table 1). Landsat TM imagery was used to derive the firn lines in 1990, 1992, 1994, 1995, 2001, 2006, and 2011, and Landsat ETM+ imagery was used to derive the firn lines in 2003, 2007, and 2009.
TABLE 1
Landsat TM and ETM+ data for firn line altitude measurement, 1990–2011, for path 135, row 33.a
The firn line position was obtained by using ASTER GDEM, which are the most accurate DEM for firn line altitude extraction. A 30-m resolution ASTER GDEM was acquired on 17 October 2011. The vertical and horizontal accuracies of the ASTER GDEM were 20 m and 30 m, respectively.
RADIOMETRIC CALIBRATION
Radiometric calibration, which includes prelaunch and post-launch tasks using laboratory and onboard methods, enables the use of full Landsat data in a quantitative sense (Saunier and Rodriguez, 2006; Thorne et al., 1997; Barsi et al., 2007). Radiometric calibration converts 8-bit satellite image digital numbers (DNs) to the top of atmosphere (TOA) reflectance or spectral radiance at the satellite's aperture. For atmospheric corrections, we calculated the radiance at the satellite's aperture for Landsat TM using the following equation:
where i is the Landsat TM specific band, L is the radiance at the satellite's aperture (w m-2 sr-1 µm-1), DN is the satellite image digital number, B is the bias of the sensor, and G is the gain of the sensor (w-1 m2 sr µm).
For Landsat ETM+, we derived the radiance at the satellite's aperture using Equation 2:
where Lmin and Lmax are the minimum and maximum spectral radiance band limits; DNmax is the maximum satellite image digital number (255); and DNmin is the minimum satellite image digital number.
ATMOSPHERIC CORRECTION
FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) is an atmospheric correction software that is based on MODTRAN4+ and was developed by the U.S. Air Force Research Laboratory at Hanscom Air Force Base and Spectral Sciences, Inc. in Burlington, Massachusetts, to support hyperspectral sensors (e.g., Hyperion, AVIRIS, and HyMap) and multispectral sensors (e.g., Landsat, SPOT, and ASTER) (Adler-Golden et al., 1999). FLAASH quantifies atmospheric properties, including the water vapor column, aerosol, and cloud optical depths, as well as surface properties such as albedo, altitude, and temperature (Kruse, 2004). FLAASH first calculates the spectral radiance, L, for any given sensor pixel in the solar wavelength range. The equation is as follows:
where ρ is the pixel reflectance, ρe is an average reflectance for the pixel and its adjacent pixels, S is the spherical albedo of the atmosphere, La is the radiance backscattered by the atmosphere, and A and B are coefficients that rely on atmospheric and geometric conditions but not on the surface conditions.
After we performed the water retrieval, we solved the pixel surface reflectance for all sensor channels by Equation 3. The average reflectance, ρe, was computed from the averaged radiance image Le using the approximate equation
CONVERTING NARROWBAND ALBEDO INTO BROADBAND ALBEDO (NTB CONVERSION)
From the Landsat narrowband calculation, broadband albedo depends mainly on the relation between TOA reflectance and the measured ground broadband albedo. We used the following linear equation:
where CH1, CH2,…, CHn is the specific band albedo of the satellite, and a0, b1, b2,…, bn is the polynomial regression coefficient that is determined by the surface and atmosphere conditions. This method has common problems of poor universal application because of large temporal and spatial differences. Greuell et al. (2002) used aircraft and near-surface measurements (including Landsat TM Bands 2 and 4), data sets from Iceland and Greenland, and data sets with ground-based albedo measurements from Switzerland and Antarctica to develop a method for converting narrowband albedo to broadband albedo for glacier ice and snow. The residual standard deviation of this method is 0.011 and it can be used without surface classification. The equation is as follows:
where αTM2 and αTM4 represents the albedo of the Landsat TM/ETM+ Band 2 and Band 4, respectively.
FIRN LINE ALTITUDE, FIRN ZONE AREA, AND FIRN ZONE AREA RATIO EXTRACTION
ASTER GDEM data registration is conducted by Envi 4.7 software and 30-m resolution Landsat TM/ETM+ data. First, the ground control point (GCP) in each Landsat image was selected, and then the same point was selected in the ASTER GDEM. The GCPs had to be feature points, such as a ridge, a river bend, or bifurcation and a road junction. The GCPs were distributed around but excluded glacial margins because of their annual variability (McFadden et al., 2011). Once we obtained at least 30 GCPs that were uniformly distributed in the whole image, and by keeping the root mean square to less than 1, then we could fit them with a quadratic polynomial method.
Contours were extracted at intervals of 10 m from the ASTER GDEM. In addition, each extracted contour layer was overlaid with the broadband albedo image, and then the glacier surface was classified according to the difference between the firn and bare ice albedo by Arc GIS software (Fig. 2). There was a transition zone between the bare ice zone and firn zone, and the boundary between the firn zone and the transition zone was firn line. The threshold criteria for separating the bare ice, firn, and transition zones on the glacier by albedo are shown in Table 2. The position of each firn line was determined by the contour layer. The firn line did not necessarily coincide with the contour, which tended to span several contours. In this case, the mean contour nearest to the firn line was decided as the firn line altitude.
This method predicted the variations of ice and firn areas on the Qiyi Glacier from 1990 to 2011 (Fig. 2). The firn area was minimal while the ice area was maximal in 2001, and there was severe melt with fewer firn zones on the upper parts of the glacier in that year. In 1994 and 2007, the differences between the firn zones were less than 50 m. In contrast, in 2001 the firn zones notably decreased and the firn line altitude increased by ∼200 m.
TABLE 2
Albedo thresholds of different boundaries on the Qiyi Glacier from 1990 to 2011.
The firn line was not necessarily always distributed along a contour, but it was located along the firn zone. We used the estimated firn line altitude to separate the glacier surface into firn zones and others based on Figure 2. The glacier boundary was mapped by the Landsat data and the glacier area was calculated, and the firn zone area was estimated using the firn zone area and the glacier area. The firn zone area ratio was obtained by the ratio of the firn zone area and the glacier area.
Results and Discussion
TEMPORALAND SPATIAL CHANGES OF THE FIRN LINE ALTITUDE
Overall, the firn line altitude of the Qiyi Glacier increased slowly from 1990 to 2011 (Fig. 3) and ranged from 4540 m a.s.l. in 1992 to 5000 m a.s.l. in 2001, a variation of 460 m.
Meteorological data from the Tuole weather station (38°48′ N, 98°25′ E, 3368 m a.s.l.) were used to verify the retrieved firn line altitude because it is the nearest weather station to the Qiyi Glacier. Variations of summer (June to August) air temperature and annual precipitation at the Tuole station from 1990 to 2011 are shown in Figure 4, and they both increased during this period. The maximum summer air temperature and annual precipitation were 11.4 °C in 2011 and 404.4 mm in 1998, respectively.
Variation of firn line altitude correlated very well with the summer air temperature. The firn line altitude increased with summer air temperature increase (correlation coefficient r = 0.7), while r between the firn line altitude and annual precipitation was -0.4. This means that the summer air temperature was the key factor in firn line altitude change. The firn line altitude was at maximum in 2001, and according to the summer air temperature and annual precipitation data at the Tuole station, during that year the summer air temperature was higher than in any other years except 2006 and 2010. The annual precipitation was lower in 2001 than in 2010, but was higher in 2005 relative to 2000.
The largest decline in the firn line altitude was observed during 2001 to 2003 (down by 330 m); the second largest decline was observed during 1990 to 1992 (down by 220 m). The firn line altitude decreased by 90 m from 2007 to 2009. The reason was that the summer air temperature decreased obviously from 2001 to 2003 (1.2 °C), 1990 to 1992 (0.6 °C), and 2007 to 2009 (0.8 °C), and the annual precipitation increased during these periods.
FIGURE 3.
Variation in retrieved firn line altitude in Qiyi Glacier from 1990 to 2011. The dashed lines are linear trends of the variations in firn line altitude.
FIGURE 4.
Variations in summer air temperature and annual precipitation at the Tuole meteorological station from 1990 to 2011.
Conversely, the firn line altitude increased in some years, including 1992–1994, 1994–1995, 1995–2001, 2003–2006, and 2009–2011. The most significant increase of firn line altitude occurred from 1992 to 1994 (230 m) compared to other years when it increased only 110 m to ∼130 m. This is because the air temperature was at minimum in 1992 and increased the most (1.2 °C) from 1992 to 1994, while the annual precipitation reduced obviously (79 mm) relative to other years.
Field observations of the Qiyi Glacier mass balance began in 1958 (Wang et al., 1985) and were discontinued in 1978 (Liu et al., 1992). Since 2001, the measurements have resumed (Pu et al., 2005). Typically, mass balance is determined by digging snow pits in the accumulation zone and using measuring stakes in the ablation zone on Qiyi Glacier. In a study by Pu et al. (2005), 14 cross-sectional grids for mass balance measurement were laid on the glacier surface. Below the ELA (∼4700 m), the spacing of the stakes was ∼100 m in the cross section and ∼50 m in the longitudinal profile, while above the ELA, two stakes and one longitudinal profile were laid. The stakes generally covered the whole glacier, and fallen stakes had to be replaced. General observations, such as the stake length, snow depth, snow density, the thickness of super-imposed ice, and the depth of the dirty layer, were recorded every five days after the stakes were laid.
However, the study of Pu et al. (2005) produced much less data compared to our longer study period. Statistical models using data from meteorological stations, and correlating ELA and temperature and precipitation, are useful for reconstructing ELA. Wang et al. (2010) used 16 years of observed ELA data on the Qiyi Glacier to create a statistical model between ELA and warm season (September, July, and August) air temperature (Tw) and cold season air (January, February, March) precipitation (Pc). Their equation is:
The ELA was reconstructed of the Qiyi Glacier from 1961 to 2008 according to Equation 7, and simulated ELA was close to the measured ELA (Wang et al., 2010). Thus, we used the measured ELA and the simulated ELA to assess the relationship between firn line altitude and ELA (Fig. 5). The actual firn line altitude retrievals correlated fairly well with the ELA simulations (r = 0.6), and the variation in the estimated firn line altitude was in good agreement with measured and simulated ELA changes. There are large differences between retrieved FLA and measured ELA. The retrieved firn line altitudes were lower than the simulated ELAs in 1992, 1994, 2003, 2006, and 2007. The simulated ELA was 318 m higher than the retrieved firn line in 2006, and the simulated ELA was 140 m higher than the retrieved firn line in 2003. In 2001, little difference was observed in a range of one pixel. Retrieved firn line altitudes in 1990, 1995, and 2001 were greater than both the measured ELAs and the modeled ELAs. Because the satellite image is only a snapshot picture of the glacier surface, it is quite sensitive to the glacier surface snowing. Moreover, it should be noted that the firn line altitude does not coincide with the ELA in the presence of superimposed ice. Simulated ELAs are generally higher than retrieved firn line altitudes; in other words, the firn line altitude calculated from Landsat imagery underestimated the ELA. This is possibly because snow covered the glacier in June, July, or August, and superimposed ice was attached on the glacier. Therefore, the retrieved firn line altitude may not be an approximation of ELA for a continental glacier. However, the variations in retrieved firn line altitude agreed with the ELA changes. Consequently, the retrieved firn line altitude can be used to indicate ELA.
Linear fitting showed that summer air temperature and annual precipitation increased by ∼0.66 °C and 46.8 mm during 1990–2011, respectively. In general, increasing precipitation led to decreasing firn line altitude. This means that the increasing tendency in firn line altitude from 1990 to 2011 could be attributed to rising air temperature. Based on the sensitivity of the Qiyi Glacier ELA to warm season (September, July, and August) air temperature (172 m °C-1), the 0.66 °C rise in air temperature could have caused the firn line altitude to increase by 113 m at Qiyi Glacier, which was reasonably close to the estimated firn line altitude increase of ∼99 m according to the linear fitting method.
FIGURE 5.
(a) Comparison of retrieved firn line altitude (FLA) with measured equilibrium line altitude (ELA) and simulated ELA from the statistical model, (b) Variations in the Qiyi Glacier FLA and simulated ELA from 1990 to 2011.
These results reconfirmed that summer air temperature is the dominant factor leading to firn line altitude change. Observations from glaciers all over the world and averaged over the long term indicate that ELA increased ∼200 m from 1961 to 1998 (Dyurgerov, 2001). This is due to global warming (Wang et al., 2010).
VARIATION IN FIRN ZONE AREA AND FIRN ZONE AREA RATIO
Based on the firn line altitude, we calculated the firn zone area and the ratio between the firn zone area and the glacier area. Here, this ratio is called the firn zone area ratio for convenience. As shown in Figure 6, the firn zone area decreased with the increase of the firn line altitude, and vice versa. The firn zone area markedly declined from 1990 to 2011. We used a linear regression to fit its variation trend, and found that the firn zone area decreased by ∼1.19 km2 at a change rate of 0.057 km2 a-1. The firn zone area of the glacier ranged from 3.08 km2 (1992) to 0.58 km2 (2001), a decline of 2.5 km2.
The change trend of the firn zone area ratio was basically identical to the variation of the firn zone area. The firn zone area ratio decreased from 76.2% to 18.1% during the study period. Temperature was the major factor controlling the firn zone area and the firn zone area ratio.
Conclusions
Firn line altitude is an indicator of ELA, and it plays an important role in predicting climate behavior and global climate change. In our study, we retrieved the firn lines and determined the firn line altitudes on Qiyi Glacier from 1990 to 2011 using Landsat TM/ETM+ and ASTER GDEM imagery. We calculated the firn zone areas and the firn zone area ratios, and analyzed variations of the firn line altitude, the firn zone area, and the firn zone area ratio. We have drawn the following conclusions:
The firn line can be identified by assessing the differences in albedo between various surfaces on a glacier.
The firn line altitude of the Qiyi Glacier increased during 1990 to 2011, generally moving from 4540 m a.s.l. to 5000 m a.s.l. The firn line altitude correlated well with measured and simulated ELA.
The firn zone area and firn zone area ratio decreased (increased) with the increase (decrease) of firn line altitude during 1990 to 2011 (the firn zone area and the firn zone area ratio had a negative correlation with the firn line altitude). The firn zone area of the Qiyi Glacier declined from 3.08 km2 to 0.58 km2 during the study period. Linear fitting demonstrated that the firn zone area decrease rate was 0.057 km2 a-1. The firn zone area ratio decreased from 76.2% to 18.1%- during the study period.
Air temperature rise over the past 21 years was the main cause of firn line altitude increase and firn zone area and firn zone area ratio decrease.
Acknowledgments
This work was supported by the “Strategic Priority Research Program (B)” of the Chinese Academy of Sciences (Grant No. XDB03030204), the National Natural Science Foundation of China (Grant Nos. 41190084 and 41371095), and the State Key Laboratory of Cryospheric Science (NO. SKLCS-OP-2014-03). We thank the Earth Resources Observation and Science (EROS) Data Center—National Aeronautics and Space Administration (NASA)/ U.S. Geological Survey for the use of available satellite imagery. We would also like to thank the National Oceanic and Atmospheric Administration (NOAA) National Climatic Data Center for providing the meteorological data.
