1. Introduction
Mountain glaciers constitute about 3% of the Earth’s glacierized area (Reference Anthwal, Joshi, Sharma and AnthwalAnthwal and others, 2006), and while the range of the estimated total area is quite small, the estimated volumes differ considerably, from 51 × 103 to 133 × 103 km3, representing a sea-level rise equivalent of 0.15–0.37 m. Reference Combes, Prentice, Hansen and RosentraterCombes and others (2004) reported that since the early 1960s, mountain glaciers worldwide have experienced an estimated net loss of >4000 km3 w.e. In the 1990s the loss rate was more than twice as fast as during previous decades. Glaciers on the Tibetan Plateau have been shown to be highly sensitive to climatic change (Reference Liu, Shen, Sun and LiLiu and others, 2002, Reference Liu, Zhang, Inoue and Horton2007), and several studies reveal that glaciers in this region have retreated remarkably in the past two decades (Reference Yao, Wang, Liu, Pu and LuYao and others, 2004; Reference PuPu and others, 2008). Although it is still ambiguous which climatic parameters play key roles in the glacier retreat in this region, Reference Ren, Qin, Kang, Hou, Pu and JingRen and others (2004, Reference Ren, Jing, Pu and Qin2006) reported that the current glacier retreat on the Tibetan Plateau is due to the combined effects of reduced precipitation and increased temperature.
Unfortunately, because of harsh climatic conditions, high altitude, steeply sloping terrain and difficult access, regular systematic surveys of glacier variations are almost impossible at most glaciers on the Tibetan Plateau. Much of the information we have is therefore based on glacier tongue observations in situ or remote-sensing techniques. Very little information is available on mass-balance variations. Thus we are forced to rely on estimates, but even these are difficult to make due to a shortfall of long-term precipitation and temperature data series in the region.
In this study we develop a model that calculates the energy balance on glaciers using minimal input data that are routinely measured over a few decades. Xibu glacier, in the Nyainqêntanglha mountain range, Tibet Autonomous Region (TAR), is selected to test the model because there are some previously observed data (area, total length of the glacier, etc.) and it is a relatively short distance (<100 km) from the meteorological station at Lhasa, which is one of the few that have recorded reliable observations since the 1950s.
The overall objective of the study is to present a simple approach to modelling the mass balance of Xibu glacier, in order to gain a better understanding of how this glacier may respond to short-term climate variability and longer-term change. The specific objectives are to investigate the sensitivity of the glacier to changes in temperature and precipitation, to explore which of the seasonal temperature trends plays the primary role in explaining the variability and trends in the modelled glacier mass balance and to examine how important changes in large-scale weather types may have been for the glacier mass balance over the past 50 years.
2. The Energy-Balance Model
From the conservation law of energy the incoming minus outgoing energy will equal the change in energy stored. If we assume there is no horizontal transfer of energy and neglect the small portion of heat from rain, we can write the surface energy equation for a point on the glacier surface at any instant,
where M is the net energy flux at the surface, Q S is the internal heat flux within the snowpack, α is the surface albedo, Q IS is the incoming solar radiation reaching the surface, Q IL is the incoming longwave radiation, Q OL is the emitted outgoing longwave radiation, Q H is the sensible heat flux and Q E is the latent heat flux. In addition to the surface fluxes, the model estimates the internal snow temperature and links the snow surface temperature to both the surface energy balance and the snowpack temperature. Full details of the model are provided in the Appendix.
3. Data and Methods
3.1. Study area
The Nyainqêntanglha mountain range stretches 740 km from west to east in the central part of the southern Tibetan Plateau. It joins Gangdise mountain in the west, and extends southeast to the Baxoila Ling peak of the Hengduan Shan. There are 7080 glaciers in total with an area of 10 701 km2 (Reference ShenShen, 2005); 28 valley glaciers are longer than 10 km. The Nyainqêntanglha serves as the watershed for two major water systems, the Yarlung Tsangpo and Nu rivers. The Yarlung Tsangpo river runs 2057 km through Tibet, and passes into India, where it is known as the Brahmaputra River. The Nu river runs through the plateau into Yunnan province, China, and then into eastern Myanmar (Burma). Figure 1a shows the Tibetan Plateau, the Nyainqêntanglha range and Xibu glacier.
Fig. 1. (a) Topographical relief of the Tibetan Plateau; the white box indicates the location of the Nyainqêntanglha range, with the arrow pointing to the terminus of Xibu glacier. (b) Xibu glacier; colour code indicates the altitude (m a.s.l.). (c) Map of Lhasa region, also showing Damxiong county (grey area), Lhado village and part of the Nyainqêntanglha. The circle and asterisk indicate Xibu glacier and the stations and village, respectively.
Xibu glacier lies in the central part of the Nyainqêntanglha. It is 12.7 km long, has an average width of 1.44 km and a surface area of 18.3 km2. Elevation ranges from 5070 m at the terminus to 7160 m at the head (Reference LiLi, 1986). Figure 1b shows the topography of Xibu glacier. The glacier is located between Lhado village and Damxiong county, northwest of Lhasa. Figure 1c is a map of the Lhasa region. The distance between the glacier and Lhado village is 23 km, between the glacier and Damxiong station is 49 km, and between the glacier and Lhasa station is 86 km.
3.2. Data sources
Daily mean temperature, precipitation, relative humidity, wind-speed and cloud-cover records from Damxiong station (49 km from the glacier) were obtained from the TAR Meteorological Bureau for the time period 2001–06. However, this dataset is short, and daily mean temperature, relative humidity, precipitation, wind speed and cloud cover (1955–2006) from Lhasa (86 km from the glacier) were also obtained from the TAR Meteorological Bureau. Table 1 shows geographical details for the selected stations and the glacier.
Table 1. Geographical details for Xibu glacier and meteorological stations
We estimate climatic information for the glacier over a longer time period by relating the Damxiong station data to the longer time series of Lhasa (section 3.4) and then relating the Lhasa data to large-scale weather patterns using k-mean clustering (section 3.3). The reason for the latter step is that we want to separate any long-term temperature trends from variability due to large-scale circulation (see Reference CaidongCaidong, 2008, for the effect of circulation variability on long-term temperature trends over the Tibetan Plateau). This procedure also helps fill gaps in the Lhasa dataset.
For information on atmospheric circulation above the plateau, we use the US National Centers for Environmental Prediction (NCEP)/US National Center for Atmospheric Research (NCAR) reanalysis data (Reference KalnayKalnay and others, 1996) for the period 1955–2006. The spatial resolution is 2.5° × 2.5°.
3.3. Linking Lhasa station data with large-scale circulation types
We separate the Lhasa temperature variability due to changes in large-scale weather types and long-term trends that may be assigned to other causes. This is done by applying a standard k-mean cluster method to the large-scale circulation. The k-mean method is an iterative clustering procedure that consists of partitioning the data into k-clusters to minimize the sum of within-cluster variance (Reference MacQueen, Le Cam and NeymanMacQueen 1967; Reference Diday, Simon and FuDiday and Simon, 1976; Reference CaidongCaidong, 2008).
The clusters represent different circulation types with variations in their frequencies from year to year. Each weather type has a deviation in temperature, precipitation, humidity and wind from the seasonal mean. To avoid influences of seasonal changes in the meteorological variables, we calculate the anomalies related to each weather type as deviations from the mean 1955–2006 state,
, in a 30 day window centred on the chosen day. The mean anomaly for the kth cluster,
, can then be written as
where Nk is the number of days in cluster k for a given year and Xn is the temperature/precipitation/humidity/wind speed at day n belonging to cluster number k.
As with all clustering methods, k-mean clustering is dependent on the choice of similarity measure (the distance measure) in assessing the similarity between daily atmospheric circulation patterns. In addition, the region over which to do the clustering, the meteorological parameter to use and the height level must all be chosen subjectively. In order to find a measure of how good the clustering is in describing local meteorological variables, we used temperature and precipitation criteria that will provide information on whether the clusters discriminate between warm and cold, as well as wet and dry, conditions. For temperature we use the information criterion IT , where
is the mean temperature anomaly for the kth cluster and Nk
is the number of days in cluster k. If IT
is zero, the cluster has no skill in differentiating between warm and cold days. For precipitation we use the information criterion IP
, where
n thr,k is the number of rainy days above a chosen threshold thr (1 mm) in the kth cluster, Nk is the total number of days within the kth cluster and p thr is the probability of total rainfall above the specified threshold at all times. If IP is zero, the clusters have no skill in differentiating between wet and dry days (the number of rainy days for each weather type is the same as the average number of rainy days).
The last point to consider is the area and meteorological parameter to choose for the clustering. As these choices are not obvious, three different meteorological parameters that are possibly related to the large-scale temperature and precipitation were tested: geopotential height, wind speed (the u and v components together) and vertical velocity. The clustering was tried at three different levels (700, 600 and 500 hPa) and over two different areas: a large area (0–50° N, 50–135° E), and a smaller area (20–45° N, 70–110°E) spanning just the Tibetan Plateau and surrounding regions. In addition, two distance measures, squared Euclidean and linear correlations, were tested (Reference CaidongCaidong, 2008).
The information criteria, IT and IP , were calculated for 2–20 clusters with the chosen regions, meteorological variables and pressure levels. This gave a total of 6684 possible cluster configurations. The conclusion, comparing IT and IP for the different configurations, indicated that clustering using the wind (u and v components together) gave the best results. The squared Euclidean distance was generally better than using linear correlations for the similarity measure. Using mid-troposphere data was generally better than using near-surface level data, and the result was usually better for the small than for the large area. Thus, we selected 500 hPa wind speed for the region 20–45° N, 70–110° E for the clustering. However, the information criteria did not clearly indicate the number of weather types to choose, so in addition to the information criteria we used the correlation between the monthly temperature and precipitation anomalies, based on weather types and observed Lhasa monthly temperature and precipitation anomalies, as a measure of goodness. Table 2 shows the correlation using different numbers of clusters. The correlations suggest that the number of clusters that may give fairly good descriptions of the monthly anomalies is 15 in the wet season and 12 in the dry season.
Table 2. Correlation coefficients between monthly temperature and precipitation anomalies based on weather types (Equation (2)) and the observed Lhasa monthly temperature and precipitation anomalies using different numbers of clusters
A poor correlation is found between observed precipitation anomalies and precipitation anomalies estimated based on the weather types during the wet season. This reflects the local nature of the Tibetan Plateau rainfall. Orographic effects and local convection have a strong influence on the distribution of precipitation. This leads to weak relationships between the amounts of local precipitation and large-scale circulation, especially during the wet season. However, we find that the small precipitation amounts in the dry season are weakly related to large-scale weather types (correlation of 0.47). Higher correlations are found between observed temperature anomalies and anomalies based on weather types (0.62 and 0.79 for the wet and dry seasons respectively).
The suggested weather types (15 and 12 in the wet and dry seasons respectively) and a predictor indicating longterm trends that are not related to variations in weather types (this term may be seen as a proxy for systematic climate changes due to an increased greenhouse effect or other systematic changes) are used to model the monthly Lhasa temperature:
The correlation coefficients between observed and modelled temperatures based on the weather types are 0.79 and 0.62 for the dry and wet seasons respectively (Fig. 2). However, changes in the weather types do not explain the long-term warming (0.44°C (10 a)−1 (dry season) and 0.23°C (10 a)−1 (wet season)) clearly visible in Figure 2 for the 1955–2005 period. Both warming trends are statistically significant at the 0.05 level. The regression indicates that the recent warming over the plateau is not related to any systematic changes in weather types, while the year-to-year fluctuations are.
Fig. 2. Lhasa observed (a) dry- and (b) wet-season temperature anomalies (solid curves) and their long term trends during 1955–2005 (dashed-dotted lines), along with a regression temperature curve (dashed) using weather types as predictors and a predictor indicating long-term trends not related to the variations in weather types (Equation (5)). r is the correlation coefficient between observed and regressed temperatures.
3.4. Linking Damxiong station data with the longer time series of Lhasa station data based on weather types
The Lhasa meteorological station is 86 km from the glacier. In order to link the temperature, humidity and precipitation over the glacier to the longer time series of Lhasa, we use the data from Damxiong station which is 49 km from the glacier; however, these data cover only a 5 year period (2001–06). We relate these data to the Lhasa data based on weather type through a simple linear regression
The correlation coefficient between the daily observed temperature and the regression equation is 0.99, indicating that temperatures vary in the same way at Lhasa as at Damxiong. The same is done for relative humidity (correlation of 0.88). However, there are poor correlations between the daily precipitation, wind speed and cloud cover at Lhasa and Damxiong, so a regression line for this is not useful and we do a simple correction between the rainfall amounts at the two stations:
where
and
are the average annual precipitation, wind or cloud cover for the 2001–06 period when data exist for Damxiong, and X
Lhasa,WT is the daily precipitation, wind or cloud cover at Lhasa.
3.5. Linking the modelled Damxiong station data with local glacier conditions
After linking the long-term data at Lhasa to the closer station of Damxiong, the final step is to link the estimated Damxiong parameters to conditions on the glacier. To do this we assume a linear lapse rate, Γ, of 0.64°C (100 m)−1 and a precipitation/elevation (increase) gradient, ΔP, of 5%(100 m)−1 based on the observations of Reference LiLi (1986) from the region. Thus, temperature at the glacier is given by:
where h is the height difference between the glacier point and Damxiong station, and T Glacier and T Dam,mod are the temperature at the glacier and the modelled temperature at Damxiong (Equation (6)) respectively.
Precipitation at the glacier is estimated from
where P Glacier and P Dam,mod are the precipitation at the glacier and the modelled precipitation at Damxiong (Equation (7)) respectively.
For relative humidity, we assume a constant mixing ratio with height up to cloud base. The daily mean relative humidity based on weather type, RH, at the glacier can then be calculated as
where e S is the saturation vapour pressure and p is the pressure given by the hydrostatic equation.
Since we have no information on wind speed or cloud cover at the glacier, we assume the wind and cloud cover are the same at the glacier as the modelled data at Damxiong.
4. Mean Climatic Conditions at the Glacier
The glacier is divided into 68 elevation zones with a 100 m vertical resolution from 5000 to 5500 and 6100 to 7200 m a.s.l. and a 10 m vertical resolution near the equilibrium-line altitude (ELA). The mass balance is calculated for each zone. To calculate the long-term mean mass balance of the glacier, we force our model repeatedly with a mean seasonal cycle consisting of the daily 30 day running mean meteorological parameters averaged over the 1966–2005 period (Equations (8–10)).
Using the mean climatic conditions, the calculated ELA is 577 ± 5 m a.s.l. Daily means of the meteorological parameters over the ablation and accumulation zones are plotted in Figure 3. Annual mean air temperature in the ablation zone is −5.8°C, exceeding the annual mean air temperature in the accumulation zone by 7.1°C. The modelled snow surface temperature indicates that the glacier starts to melt on 8 May and ceases melting on 23 September (this is an average for the ablation zone).
Fig. 3. (a) Seasonal cycle of daily mean air temperature, T a, daily mean snow surface temperature, T snow, and daily maximum temperature, T max. The dotted lines link the snow surface temperature (at 273.15 K), the surface melt start day (day number 128) and the melt end day (day number 266). (b) Daily mean precipitation. (c) Daily mean relative humidity over the accumulation (grey curve) and ablation (black curve) zones. (d) Daily mean wind speed. Daily mean cloud cover in (c) and daily mean wind speed in (d) are the same values for all 68 elevation zones on the glacier.
Daily precipitation in both the accumulation and ablation zones is shown in Figure 3b. The average annual precipitation is 965 and 709 mm w.e. over the accumulation and ablation zones respectively. There is a peak during August, with low values throughout the dry season. Based on Equation (A19) in the Appendix, ∼63% of annual precipitation falls as snow in the ablation zone and ∼97% falls as snow in the accumulation zone.
Figure 3c shows that the relative humidity is nearly 100% in the accumulation zone and ranges between 77% and 98% in the ablation zone. The near 100% value in the accumulation zone might be an overestimation due to the constant mixing ratio assumption of Equation (10).
Cloud cover and wind speed are assumed constant with elevation because of the lack of data, i.e. they are the same in the accumulation and ablation zones. Figure 3c shows the seasonal cycle in cloud cover, which is similar to the seasonal cycle in precipitation, with a peak during August (83%) and a lower value during the dry season. The average cloud cover is 55%.
Figure 3d shows day-to-day fluctuation of wind speed from 1.6 to 3.4 m s−1. The average wind speed over the period is 2.4 m s−1. For the higher parts of the glacier this may be an underestimate.
The radiative energy fluxes (defined as positive towards the surface) in the ablation zone (average over the zone) are shown in Figure 4a. The largest flux over the seasonal cycle is the outgoing longwave radiation, with a daily average of −289 W m−2. It is nearly constant over the period 8 May–23 September since the surface was at the melting temperature. The second largest energy flux is incoming shortwave radiation, ranging from 169 to 271 W m−2 (daily average 215 W m−2). The average reflected shortwave radiation is −92 W m−2, implying an average surface albedo of 0.43. Incoming longwave radiation shows the largest day-to-day fluctuations (151–269 W m−2). The turbulent fluxes are much lower (Fig. 4b), being slightly negative (0 to −25 W m−2) during winter (due to lower snow temperatures than air temperatures) and positive in summer (generally <10 W m−2). These low values are due to both low wind speed and stable stratification in winter.
Fig. 4. Average radiation fluxes and turbulent heat fluxes over the ablation zone for the period 1966–2005. (a) Daily incoming shortwave radiation, Q IS, incoming longwave radiation, Q L, net shortwave radiation, Q NS, net longwave radiation, Q NL, reflected shortwave radiation, Q RS, and outgoing longwave radiation, Q OL. (b) Daily sensible heat fluxes, Q H, and latent heat fluxes, Q E.
The average energy flux in the accumulation zone is given in Figure 5. The daily average net shortwave radiation is 51% less in the accumulation zone than in the ablation zone, due to more snow and therefore higher albedo. The outgoing longwave radiation is 13% less in the accumulation zone due to colder conditions, and the net longwave loss is also 13% lower in the accumulation zone than in the ablation zone.
Fig. 5. Average radiation fluxes and turbulent heat fluxes over the accumulation zone for the period 1966–2005. (a) Daily incoming shortwave radiation, Q IS, incoming longwave radiation, Q IL, net shortwave radiation, Q NS, net longwave radiation, Q NL, outgoing shortwave radiation, Q OS, and outgoing longwave radiation, Q OL. (b) Daily sensible heat fluxes, Q H, and latent heat fluxes, Q E.
The annual average net mass balance, b n,z, is calculated by running the model with the seasonal cycle of the meteorological parameters (Fig. 3) repeated for 40 years. In addition, using the Shuttle Radar Topography Mission data (Reference FarrFarr and others, 2007) and the total glacier area data of Reference LiLi (1986), we have divided the total area into different topographic zones and derived the Xibu glacier area-weighted specific net mass balance, B n,z (m3 w.e), calculated by integrating the net mass-balance values over the glacier surface area within different elevation zones. The results (Fig. 6a) show a steep gradient in the annual average net mass balance in the ablation zone, with more melting in the lower part of the glacier.
Fig. 6. (a) Average annual net mass balance versus elevation. (b) Area-weighted specific net mass balance versus elevation. (c) Glacial surface area ratio for every 100 m height interval.
The calculations indicate that mass balance in the lower part of the ablation zone is about −3 to −3.5 m w.e. a−1.In the accumulation zone, mass balance is less variable with elevation, and ranges between 0.6 and 1 m w.e. a−1.The ELA is defined where the b n,z curve intersects zero. For the mean condition, the calculated ELA is 5770 ± 5 m. The values of b n,z show a large jump from ablation zone to accumulation zone. We interpret this as possibly reflecting two properties: (1) a larger portion of the precipitation falls as snow above the ELA, increasing the accumulation, and (2) the amount of snow strongly influences the surface albedo, thus reducing the melt rate. On Xibu glacier, only 5% of precipitation falls in winter (October–April). On glaciers in higher latitudes, a higher percentage of precipitation falls during winter. There is therefore relatively more snow in the ablation zone of glaciers with larger wintertime precipitation than Xibu glacier. Thus, Xibu glacier has a more constant snowline elevation late in the ablation season than glaciers in higher latitudes, resulting in a distinct jump between the ablation and accumulation region.
The area-weighted net mass balance, B n,z, demonstrates that the total net ablation is slightly larger than the total net accumulation using our estimates of mean climatic conditions. Figure 6c shows the surface area ratio of Xibu glacier sharply increasing just above the terminus, indicating that the Xibu glacier net mass balance is very sensitive to small changes in ELA.
5. Sensitivity to Changes in Climatic Parameters
To test the sensitivity of Xibu glacier to climatic changes, the ELA, annual average net mass balance, b n,z, area-weighted specific net balance in different elevation zones, B n,z , specific net mass balance in the ablation zone, B n,ab, and in the accumulation zone, B n,ac, were calculated with prescribed changes to air temperature and precipitation. The mass-balance model was run with temperature changes of ±0.5°C and ±1 °C and precipitation changes of ±20% and ±40%. The resulting changes are summarized in Table 3 and Figure 7.
Fig. 7. Annual net mass balance versus elevation and its dependence on (a) change in air temperature and (b) change in precipitation.
Table 3. The net, annual and cumulated mass balance, B n, B n,ab and B n,ac, respectively, of Xibu glacier, ELA variation with precipitation and air-temperature changes over the 40 year period 1966–2005, and their deviations, ΔB n, ΔB n,ab, ΔB n,ac and ΔELA, with respect to the mean condition (ΔP= 0% and ΔT= 0°C)
When the rate of precipitation is assumed to be constant and air temperature is increased, B n,ab becomes more negative, with a reduction of ∼15%(10 a)−1 if a 1°C warming is imposed for 40 years. The accumulation in the ablation zone is also reduced, but is somewhat less sensitive to the temperature change (∼10%(10 a)−1 reduction). The ELA is raised by 140 m with a 1°C warming, and the enlargement of the ablation area is 20% (Fig. 8).
Fig. 8. Changes in ELA and cumulative area with ablation for changes in air temperature and precipitation.
Unlike temperature change, which has more effect on the ablation area, precipitation change affects the specific net mass balance of the ablation and accumulation areas almost equally. A 40% change in precipitation gives changes in ELA of −100 and +155 m for positive and negative precipitation anomalies respectively. A 40% reduction in precipitation enlarges the ablation area by 22%. The ELA is more sensitive to a precipitation decrease than to a precipitation increase (Table 3), because of the non-linear effect of precipitation due to its impact on both accumulation and melting through the albedo changes. The results imply that the change in ELA due to a ±1°C change in temperature is equivalent to a ±35% change in precipitation.
6. Influence of Circulation Variability and Long-Term Temperature Trends on Xibu Glacier Mass Balance
The previous estimates indicate that air temperature is the dominant influence on mass-balance changes on Xibu glacier. In this section, we discuss how the estimated mass balances react to observed seasonal temperature trends and observed large-scale circulation variability.
In section 3.3 we showed how a simple regression line (Equation (5)), using weather types and a long-term trend as a proxy for trends not related to circulation changes, represented fairly well the temperature series for Lhasa (the closest station with long time series of meteorological parameters). We now use the same regression line to show how the mass balance of Xibu glacier may have reacted to temperature change during the last 50 years. We run the model based on the links between weather type and meteorological parameters given in sections 3.3 and 3.4 for the period 1955–2005 at 68 elevation zones. Using the first 10 years, 1955–65, as a spin-up period, we perform four simulations: (1) We include only the variability due to weather types (WT), i.e. keeping only the first term in Equation (5), which indicates the impact of changes in frequency of different weather types on the glacier mass balance over the past 40 years. We then perform the same calculation but add (2) the observed wet- (May–September), (3) the observed dry- (October–April) and (4) the observed wet- and dry-season temperature trends not accounted for by the weather types, i.e. the second term in Equation (5). The dry-season temperature trend is 0.44°C(10 a)−1 and the wet-season trend is 0.23°C(10 a)−1 (Fig. 2). These simulations are designated WTBTT for the simulation with both the dry- and wet-season temperature trend, WTDTT for the simulation with the dry-season temperature trend but no wet-season trend and WTWTT for the simulation with the wet-season trend but no dry-season trend.
Figure 9 shows the modelled cumulated wet-season, dry-season and annual specific net mass-balance variations over Xibu glacier for the period 1966–2005. The simulation based on weather types with no temperature trends shows fairly strong year-to-year variations in the cumulated wet-season mass balance, with a maximum accumulation of 2.5 × 106 m3 w.e. in 2003/04 and a minimum of 6.3 × 105 m3 w.e. in 1974/75. Thus the 2003/04 wet-season circulation made the glacier accumulate four times more mass than in the 1974/75 season. There is a small positive trend in the WT simulation, indicating that if there were no temperature trends, the circulation variability over the plateau should have made the glacier gain mass over the simulation period if no additional forcing changes took place during the period. The correlation between the wet-season variability and the annual net mass-balance variability is 0.98, indicating that the wet season is the main driver in the mass-balance variability of the glacier.
Fig. 9. Simulation of the mass-balance sensitivity of Xibu glacier to weather types. Wet-season mass balance (WET), dry-season mass balance (DRY) and net mass balance (NET) based on the four simulations, WT, WTDTT, WTWTT and WTBTT during 1966–2005. The grey and black colours indicate mass balance based on weather type and weather type with seasonal temperature trends respectively. The correlation coefficient between net mass balance and wet-season balance, r nw, and between net mass balance and dry-season balance, r nd, are indicated.
The WTDTT simulation nearly overlaps the curves based on the WT simulation, indicating that the dry-season (wintertime) temperature trend does not greatly affect the annual glacier net mass balance, even if the trend has been high (0.44°C(10 a)−1). In contrast, the WTWTT simulation shows a strong deviation from the WT simulation. The specific net annual mass balance decreased significantly, with a mass-balance trend of about −4.3 × 105 m3 w.e. (10 a)−1 over the period. The decrease is statistically significant at the 0.05 level. The simulation using both the wet- and dry-season trends is very similar to that with only the wet-season trend.
This result is in line with Reference Wang, Jiao, Li, Jing and HanWang and others (2004) and Reference PuPu and others (2008) who report decreased mass balance over Da Dongkemadi and Xiao Dongkemadi glaciers, Tibetan Plateau, over the past two decades.
Our simulations indicate that the wet season is more important than the dry season in explaining variations in the net balance of Xibu glacier. Changes in wet-season temperature affect both surface melting and the partitioning between snow and rain, which again affects surface melting through albedo changes.
Comparing the mean (1966–2005) wet-season (May–September) energy fluxes in the ablation zone between the WT and WTWTT simulations shows a large change in net shortwave radiation gain, a trend of 13 W m−2 (10 a)−1. This is due to changes in albedo. Net longwave radiation loss is reduced by 9 W m−2 (10 a)−1, resulting in an increased net radiative gain of 22 W m−2 (10 a)−1. The turbulent fluxes add another 10 W m−2 (10 a)−1. Thus the 0.23°C (10 a)−1 wet-season warming results in a total wet-season net energy flux increase of 32 W m−2 (10 a)−1 in the ablation region.
Figure 10 shows b n,z and B n,z for the period 1966–2005 as a function of elevation. ELA is defined where the b n,z curve intersects zero. Based on the WT simulation, the mean ELA over the period 1966–2005 is 5492 ± 5 m. This is approximately 100 m lower than the mean for the WTWTT simulation (5590 ± 5 m). Thus a constant wet-season trend of 0.23°C (10 a)−1 results in an ELA rise of 49 m (10 a)−1.
Fig. 10. Annual net mass balance (a) and area-weighted specific net mass balance (b) versus elevation, and calculated ELA, based on the four simulations, WT, WTDTT, WTWTT and WTBTT during 1966–2005.
Based on the WT simulation, the mean 1966–2005 cumulated net glacier mass, B n, for Xibu glacier is 5.1 × 106 m3 w.e. (1121 mm w.e. of b n). B n over the ablation area is −4.6 × 106 m3 w.e. (–11 482 mm w.e. of b n). B n over the accumulation area is 9.7 × 106 m3 w.e. (12 603 mm w.e of b n). If the wet-season air temperature is increased by the seasonal trend (0.23°C(10 a)−1), the cumulated net glacier mass loss, B n, is increased by 40% over the ablation area and decreased by 21% over the accumulation area (Fig. 10b).
In the model the snow surface starts to melt when the daily mean snow surface temperature reaches melting point. The mean 1966–2005 start and end dates for the melting season, the simulations with weather types and no trends (WT) and weather types and wet-season trend (WTWTT) are listed against elevation in Table 4.
Table 4. Mean length of ablation season, melting start and end dates for different elevations for the period 1966–2005, based on simulations with weather types with no trend, WT, and weather types with wet-season trends, WTWTT
Based on the WT simulation, the lengths of the ablation season at the glacier terminus (5000 m a.s.l.) and at the ELA (5492 m a.s.l.) are 132 and 45 days respectively. The ablation season is just 3 days at an elevation of 5625 m. By adding the wet-season trend (simulation WTWTT) the ablation season was lengthened by 2 days (10 a)−1 at the glacier terminus and by 2.75 days (10 a)−1 at the ELA.
7. Discussion
Lack of data from the glacier is the primary problem for this study. The altitude of the glacier makes it almost impossible to carry out traditional fieldwork. Mass-balance studies at high, remote glaciers can be made by remote sensing or by calculating the energy balance based on data from lower altitudes and from other glaciers. These data can then be used to create synthetic time series that have some relevance for the glacier in question. We use this method in this study and we take a further step, linking the data from lower altitudes with weather types to distinguish between circulation-induced variability and long-term trends that might have other causes (e.g. global warming). In the following we discuss the accuracy of the different input parameters relevant to our simulation of the mass balance of Xibu glacier.
7.1. Reliability of input parameters and main assumptions
The correlation of daily precipitation between Lhasa and Damxiong stations was only 0.47 and the link with weather types was even weaker. This demonstrates that the weather systems delivering the precipitation are small and, to a large extent, modified by the topography and local convection. On a monthly basis, however, the precipitation in Tibet is largely affected by large-scale circulation systems (the monsoon system).
The change in precipitation with height in an area with complex topography such as Tibet is complex and probably very dependent on wind direction and season. As most of the precipitation falls in the summer season (90%), correct estimation of precipitation is important for the accumulation area, which is the only part of the glacier where the precipitation falls as snow during summer.
Our study of the relationship between observed air-temperature anomalies at Lhasa and estimates based on weather type gave surprisingly good results (correlation of 0.62 and 0.79 for the wet and dry seasons respectively). An air-temperature lapse rate of 0.64°C (100 m)−1 is often assumed and is based on a world average. As our lapse-rate correction takes into account both altitude and location differences it is not clear that this is a reasonable value. In winter, strong inversions take place in this region and it may be that a lapse rate of 0.64°C (100 m)−1 is too high. However, we have demonstrated that in winter the rate of melting can be neglected, making the winter lapse rate less important.
We assume that the cloud cover is the same for the whole area and equal to the modelled Damxiong station data. Cloud reduces incoming shortwave radiation and enhances incoming longwave radiation. In the accumulation area, with high surface albedo, the net radiation will not vary greatly with cloud cover. In the ablation area the average surface albedo was estimated to be 0.43, and here our cloud-cover assumption is more important. In the ablation area the albedo is dependent on the thickness and structure of the snow and on the amount of accumulated dust and debris. Thus the albedo is variable. Here we face the same problem as in most other studies of energy budgets on glaciers. The variability in glacier albedo is usually not well understood, and simple parameterizations must be used to identify the main features. We believe our albedo estimates are probably as good (or bad) as in similar studies.
Our estimates of air humidity are based on data from Damxiong station and the assumption that the mixing ratio (i.e. mass of water vapour/mass of dry air) does not vary with height. This may be a fairly good estimate if no condensation takes place in the air volume, but the relative humidity in the accumulation region seems to be overestimated.
We have demonstrated that the mass of ice increases above the equilibrium line and decreases in the ablation area. This will lead to increasing glacier surface slope. We have assumed that the slope of the glacier remains constant so that the vertical component of the ice flow at the equilibrium line is equal to the accumulation rate in the accumulation area. Also, the area of the glacier remains unchanged, despite melting and accumulation, over the 1966–2005 period (fixed geometry approach). This is not correct, but no information about the bottom topography or ice thickness of Xibu glacier exists, making more accurate simulations of the glacier flow and area variations difficult.
7.2. Sensitivity tests: comparison with other studies
We carried out sensitivity simulations in which we perturbed air temperature by ±1°C and precipitation by ±40%. An air-temperature change of ±1°C resulted in an ELA change of 140 ± 125 m. This is strongly dependent on our assumption of a 0.65°C (100 m)−1 lapse rate. If no other changes occur, this gives a 153 m change in ELA. Our estimates also indicate that a ±1°C change is compensated by a precipitation change of approximately ±35%. This percentage precipitation change corresponds to an absolute change over the glacier region of ∼290 mm. This is in good agreement with the studies in Xinjian province by Reference Yao and ShiYao and Shi (1988) who showed from mass-balance modelling that a 30–40% increase in annual precipitation is needed to compensate the mass loss due to a warming of 1°C.
Comparison with other studies indicates that the sensitivity of the ELA of Xibu glacier to air temperature is slightly higher than for glaciers in the European Alps. Reference OerlemansOerlemans (2001) found a typical value of 80 m °C−1 for a set of Alpine glaciers. For 35 glaciers in Europe, Reference Zemp, Hoelzle and HaeberliZemp and others (2007) reported that a modelled summer air-temperature change of 11°C could cause the ELA to change by 137 ± 125 m. This is in line with our results.
7.3. Results based on weather types, 1966–2005
Simulations using variations in large-scale weather types (closely linked to air-temperature fluctuations) and a longterm trend (not linked to changes in the weather) revealed that wet-season (May–September) air temperature is the dominant influence on mass-balance variability on Xibu glacier. This is in line with palaeo-records (e.g. from tree rings), which show that the major glacier advances in the northeast Tibetan Plateau correspond with periods of low summer temperature (Reference Gou, Chen, Yang, Jacoby, Peng and ZhangGou and others, 2006). Reference Wang, Jiao, Li, Jing and HanWang and others (2004) noted that the annual mass balance of Ürümqi glacier No. 1, Tien Shan, China, was correlated with summer mean air temperature, suggesting summer temperatures are a major cause of mass-balance variation. The observed variation of glacier mass balance, ELA and glacier terminus position of Xiao Dongkemadi glacier (northeast Tibetan Plateau) has also been attributed to changes in summer temperature (Reference PuPu and others, 2008). Our simulations indicate a rapid decrease in the net mass balance, and we would expect shrinkage of the glacier. Unfortunately, however, we are not aware of any in situ data to support this reduction. Reference CaidongCaidong (2004) compared changes in the terminus of 16 glaciers in the Nyainqêntanglha range, from topographic maps (produced in 1977), with images dated 2001, from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER). There are large uncertainties both in the quality of the topographic maps and in the dating of the source information. The analysis did, however, show an overall retreat of all the glacier tongues, with an average retreat of >2000 m. This is in line with our estimates indicating a mass loss of Xibu glacier. In addition, Reference Yao, Wang, Liu, Pu and LuYao and others (2004) indicate a 35 m a−1 retreat of glacier tongues (based on two glaciers) in the area.
Recently, new mass-balance observations from Zhadang glacier, located on the northeast slope of the Nyainqêntanglha, have been reported (Reference KangKang and others, 2009). These observations span the period 2006–08 and will provide new and valuable information. We have also established a new observation site in the region, which measures meteorological parameters, and also flow from one of the many glacier-fed rivers in the region. This will hopefully aid mapping of the local meteorological conditions and the diurnal and seasonal cycle of the glacier river runoff, providing better input to the mass-balance modelling.
8. Summary and Conclusions
We have introduced a method for estimating the mass balance of a glacier based on sparse data, by linking the data to the frequency of large-scale weather types. An objective classification procedure using the k-mean algorithm was used to link the observed daily temperature and precipitation records to large-scale weather types over the Tibetan Plateau (20–45° N, 70–110° E). A clustering scheme was based on daily wind at 500 hPa. As a result, 12 weather types during the dry season and 15 weather types during the wet season were selected to account for the large-scale flow over the plateau. It was shown that the temperature variability was fairly closely related to the weather types, but the precipitation variability was not, especially in the wet season during which local convective activity seems to be the main source of variability.
An energy-balance model was applied to model the mass balance of Xibu glacier. Due to its high altitude, systematic surveys of glacier mass balance are lacking and it is impossible to obtain direct validation of the model’s performance. However, the results of sensitivity tests are in good agreement with other studies, suggesting that our main conclusions are fairly robust.
A sensitivity study based on running the mass-balance model with the average seasonal cycle, based on observed meteorological parameters from surrounding stations, shows that an air-temperature deviation of ±1°C results in an ELA deviation of 140 ± 125 m and an ablation zone area change of 20 ± 14%. An annual precipitation deviation of ±40% led to an ELA deviation of −100 ± 155 m and an ablation zone area change of −13 ± 22%.
A sensitivity study based on weather type from the period 1966–2005 indicates that variability in the frequency of weather type during the wet season is the main contributor to year-to-year variability in the annual mass balance. In addition, a long-term trend in air temperature, not related to change in weather type, gave a systematic reduction in the mass of the glacier. Inserting the observed wet-season temperature trend ((0.23°C (10 a)−1 observed at Lhasa, 86 km from the glacier) into the simulation led to an ELA rise of 49 m (10 a)−1. The ablation season was lengthened by 2 days (10 a)−1 at the glacier terminus (5000 m a.s.l.) and by 2.75 days (10 a)−1 at the ELA.
Acknowledgements
This work was supported by University Network Norway–Tibet and the quota scholarship at the University of Bergen. We thank the TAR Meteorological Bureau for providing station data. We thank Y. Gjessing for useful comments that helped improve the paper.
Appendix
The Mass-Balance Model
Incoming shortwave radiation
The solar radiation striking a horizontal surface, Q IS, is given by
where S is the solar irradiance at the top of the atmosphere, T K is the net sky transmissivity and the zenith angle, ω, is defined as
where ϕ is the latitude, τ is the hour angle and δ is the solar declination. τ is calculated as
where dec is decimal hour, and δ is approximated by (Reference BourgesBourges, 1985)
with D = 2π/365.25(dy/79.346), dy being the number of days after 1 January.
To account for the scattering, absorption and reflection of shortwave radiation by clouds, we use the simple transmissivity scaling of Reference Burridge and GaddBurridge and Gadd (1974):
where k represents the clearness index. Reference Olseth and SkartveitOlseth and Skartveit (1993) suggested the empirical formula
where c is cloud cover in octas.
Snow depth and surface albedo
Measured precipitation rate, P, is partitioned into rain, P r, and snow, P s, (mmw.e.) based on air temperature, T a, using (Reference Tonnessen, Williams and TranterTonnessen and others, 1995)
where T r =3°C is a threshold air temperature above which all precipitation is rain and T b=−1°C is a threshold air temperature below which all precipitation is snow. Note that the T r and T b thresholds depend on the frequency of the temperature data. We use daily data, and therefore values that deviate significantly from zero to take into account the daily amplitude in temperature.
As no wind-blown drifting snow or accumulation due to avalanches is accounted for, the equation for the snow thickness, H (mmw.e.), is
where ρ w is water density, L f is specific latent heat of fusion, P s is the rate of precipitation as snow and H 1 is the snow thickness of the day before. Q m is heat used to melt the snow and ice:
The surface albedo equals the ice albedo, α ice (0.43 Reference PatersonPaterson, 1994), if snow is absent on the glacier, and i dependent on snow thickness if snow is present (Liu an others, 2007), i.e.
where ξ is a step function defined as
Here snow albedo, α s, is assumed to be independent of the age of the snow. To compensate for this simplification, the snow is assumed to have a standard age, and the value for ξ is reduced by 29% (Reference OerlemansOerlemans, 1992).
Longwave radiation
The downward longwave radiation is based on Reference SwinbankSwinbank (1963) with the corrections of Reference Iziomon, Mayer and MatzarakisIziomon and others (2003). To take into account clouds and altitude dependence,
where T a is the air temperature (K), σ is the Stefan–Boltzmann constant and e is the vapour pressure at screen level (hPa). X s = 0.35, Y s = 11.5 K hPa−1 and Z s is taken as 0.005 for mountain sites (Reference Iziomon, Mayer and MatzarakisIziomon and others, 2003).
The longwave radiation emitted by the snow surface is computed under the assumption that snow emits as a grey body:
where T S is the snow surface temperature and δ S is the surface emissivity, assumed to be 0.99.
Turbulent fluxes
Sensible (Q H) and latent (Q E) heat fluxes between the snow surface and air above are modelled using bulk formulas:
Using the relationship between specific humidity and vapour pressure and the ideal gas law, one obtains
where ρ a is the air density, Cp is the specific heat capacity, q a is the air specific humidity, q S is the specific humidity at the snow surface, e S is the vapour pressure at the snow surface, assumed saturated at T S, e is the vapour pressure at screen level and R d is the gas constant for dry air. K h and K e are turbulent transfer conductivities for sensible and latent heat respectively, with the adjustments of Reference Price and DunnePrice and Dunne (1976) for non-neutral conditions:
Ri is the Richardson number and z is the nominal measurement height for air temperature and wind speed. According to Reference PatersonPaterson (1994), z 0 ≈ 10−3 m for snow and ice. u is the wind speed and k v is von Kármná’s constant.
Snowpack temperature
To take into account the snowpack temperature, we use the simple formulation of Reference Fontaine, Cruickshank, Arnold and HotchkissFontaine and others (2002), but have replaced the air temperature with surface temperature in their relationship:
where T SN1 is the snowpack temperature from the day before and the coefficient (λ = 0.5) indicates the lag of the snow-pack temperature on the surface temperature response.
Snow surface temperature
Reference Tarboton, Chowdhury and JacksonTarboton and others (1995) computed the snow surface temperature, T S, from a simple heat conduction, Q, that is near steady state over an effective depth, Z e (m). The effective depth is equal to the depth of penetration of a diurnal temperature fluctuation (Reference Rosenberg, Blad and VermaRosenberg and others, 1983) considering a uniform temperature gradient between T S and snowpack temperature, T SN, approximated by
where κ. is the snow thermal diffusivity (m2 s−1), ρ s is the density of snow and C S is the specific heat of ice. The ratio κ/Z e is denoted by K S and termed snow surface conductance (m s−1). Assuming that equilibrium at the surface and the surface energy balance are given as a function of T S,
where Q le, Q h and Q e are dependent on T S. According to Reference Tarboton, Chowdhury and JacksonTarboton and others (1995), Equation (A19) can be Simplified to equation (A20) and solved for T S:
where the surface temperature, T S, is solved iteratively with an initial T* = T a and in each iteration T* is replaced by the latest T S until the difference using a new T* is very small.
