Dramatic mass loss in extreme high-elevation areas of a western Himalayan glacier: observations and modeling

Rapid climate change at high elevations has accelerated glacier retreat in the Himalayas and Tibetan Plateau. However, due to the lack of long-term glaciological measurements, there are still uncertainties regarding when the mass loss began and what the magnitude of mass loss is at such high elevations. Based on in situ glaciological observations during the past 9 years and a temperature-index mass balance model, this study investigates recent mass loss of the Naimona’nyi Glacier in the western Himalayas and reconstructs a 41-year (1973/74–2013/14) equilibrium line altitude (ELA) and glacier-wide mass loss. The result indicates that even at 6000 m above sea level (a.s.l.), the annual mass loss reaches ~0.73 m water equivalent (w.e.) during the past 9 years. Concordant with the abrupt climate shift in the end of 1980s, the ELA has dramatically risen from ~5969 ± 73 m a.s.l. during 1973/74–1988/89 to ~6193 ± 75 m a.s.l. during 1989/90–2013/14, suggesting that future ice cores containing uninterrupted climate records could only be recovered at least above 6200 m a.s.l. in the Naimona’nyi region. The glacier-wide mass balance over the past 41 years is averaged to be approximately −0.40 ± 0.17 m w.e., exhibiting a significant increase in the decadal average from −0.01 ± 0.15 to −0.69 ± 0.21 m w.e.


Results and Discussion
Observed and modeled mass loss of the Naimona'nyi Glacier. Although stake measurements display significant spatial differences even at the same elevations due to the topographical effect (e.g., aspect, shading, snow drift), such data are helpful for deriving the mean mass loss as a function of elevation during the past decade (Fig. 2). The annual mass balance gradient between 5800 to 6100 m a.s.l. ranged from ~0.25 to ~0.38 m w.e. (100 m) −1 . From point measurements and mass balance gradients, it is calculated that the ELA during six out of nine years was above the elevation of the highest measuring stake (i.e., 6100 m a.s.l.). The two lowest ELAs were estimated to be ~6080 m a.s.l. and ~6060 m a The modeling performance of annual mass balance as a function of elevations is shown in Fig. 2. The grey lines envelop the first 100 optimized results among 1500 random simulations, with a root mean square difference (RMSD) less than 0.5 m w.e. a −1 between the measured and simulated values at the stakes (Supplementary Figs S1 and S2). Because many factors (e.g., topographic shading, snow drift) can contribute to the spatial heterogeneity of mass loss, the overall agreement between measurements and model simulations demonstrates both the effectiveness of mass balance reconstruction using the simply mass-balance model and its feasibility for evaluating glacier status during the past four decades. Thus, the 41-year ELA and glacier-wide mass balances were reconstructed (Fig. 3a,b). The modeled glacier-wide mass balances are compared with the geodetic mass balances in different periods 14,15   To further illustrate the temporal changes in surface mass loss at above 6000 m a.s.l., the accumulative mass balances from the two different altitudes (6060 and 6220 m a.s.l.) were generated (Fig. 3c,d). At the elevation closer to the ice core drilling site in 2006 11 , there was only a small mass accumulation during the period from 1973/74 to the late-1980s (Fig. 3c). However, after that, there was continuous negative mass balance resulting in a total mass loss of ~4.7 m w.e. by 2013/14. The distinctive characteristics are portrayed by the less magnitude of mass gain during 1973/74-1988/89 but an accelerating mass loss trend after that period. As a comparison, the ELA increases largely from ~5969 ± 73 m a.s.l. to 6193 ± 75 m a.s.l. with a dramatic enhancement of glacier-wide mass loss from − 0.06 ± 0.15 to − 0.62 ± 0.19 m w.e. The turning point corresponds to an abrupt climate change recorded at the Burang station (30°17′ N, 81°15′ E, 3900 m a.s.l., 1973-2014), including a dramatic temperature increase and significant precipitation decrease in the late-1980s ( Supplementary Fig. S3). The positive accumulations at 6220 ma.s.l. from 1973/74 to 2013/14 (Fig. 3d) likely indicates that the lower limit at which ice cores can provide successive, high-quality paleoclimatic records in the Naimona'nyi region. Similarly to this region, a few ice cores drilled in the southern and central Tibetan Plateau 22 also found no net accumulation even at 5800 m a.s.l. This suggests that the area of accumulation is becoming a ablation zone with the dramatic rise in ELA 22 . However, due to few in-situ glacio-meteorological observations, it is hard to precisely identify the ideal elevations for ice core drilling in these regions. Thus, more glaciological studies should be carried out to detect mass loss in extreme high elevations.
Due to the lack of long-term meteorological data before 1973 at the Burang station and in surrounding regions, the grid data of temperature and precipitation from the CRU TS 3.22 dataset 23 covering the Burang station were used to evaluate the climate background before 1973 ( Supplementary Fig. S4). As expected, compared to the period of the 1970s and 1980s, the Naimona'nyi region was under relatively warm, dry conditions between the 1940s and 1960s. Based on the variations in temperature and precipitation, it can be speculated that the mean ELA during the 1940s-1960s might have been higher than the mean value of ~5935 ± 76 m a.s.l. between 1973/74 and 1982/83, but should be much lower than that of ~6212 ± 77 m a.s.l. over the last decade. Thus, it is possible that mass loss above 6000 m a.s.l. on the Naimona'nyi Glacier has occurred since at least the 1940s. In this context, the fallout radionuclides from atmospheric thermonuclear bomb tests in the 1950s and 1960s were not contained initially in the ice body at the elevation band 6025-6200 m a.s.l. of this glacier. This finding does not provide a definitive conclusion, however, because the lack of radionuclides could be due to lack of accumulation in the 1950s and 1960s or subsequent mass loss of ice containing these particles.
The mass loss above 6000 m a.s.l. on the Naimona'nyi Glacier in the western Himalayas can be attributed to unique geographical and topographical conditions. Generally, annual precipitation has a negative gradient from east to west on the Tibetan Plateau 24 . In addition, the Himalayan range acts as an effective orographic precipitation barrier and creates a strong rain lee-effect for the regions north of the range, with drier air masses and less precipitation 25 . The mean annual precipitation was only 157 mm at the Burang meteorological station dur-ing1973-2014. Even at AWS2 (5950 m a.s.l), only 221 mm of precipitation was recorded by a T-200B precipitation gauge during October 2013 to September 2014. From the perspective of energy and mass balance, precipitation contributes to mass accumulation while solar radiation promotes glacier melt. In the Naimona'nyi region, on one hand, less precipitation positively limits the glacier accumulation; on the other hand, more solar radiation energy related to the low latitude and atmosheric moisture content could arrive and possibly be absorbed for snow and ice melting at the glacier surface by an accumulation-albedo-melt feedback mechanism 26 . Furthermore, the increased air tempeature in recent decades ( Supplementary Figs S3 and S4) has also contributed much energy for glacier mass loss in higher elevations. As a result, the enhanced energy supply with less precipitation for glaciers on the lee side of western Himalayas has caused some of the highest ELAs and net mass losses at high elevations worldwide.

Summary
In this study, our in situ measurements show that even at 6000 m a.s.l., the mean mass loss of the Naimona'nyi Glacier reached as much as ~0.73 m w.e. a −1 over the past 9 years. Using a calibrated mass balance model, the variations in mass loss and ELA over the past four decades were reconstructed and analyzed. We found that annual averaged mass loss of the Naimona'nyi Glacier had increased over this time period from ~0.01 ± 0.15 m w.e. (1973/74-1982/83) to ~0.69 ± 0.21 m w.e. (2003/04-2013/14). We also observed accelerating trends of mass loss and ELA rise beginning in the late-1980s, which are associated with a shift in regional climate condition. The ELA has risen approximately 280 m in the previous four decades and reaching an average elevation of ~6212 ± 77 m in the last decade. According to the regional climate condition, the net mass loss on the Naimona'nyi Glacier above 6000 m a.s.l. has been possibly occurring since as early as the 1940s.
Our study also sheds light on the possible minimum elevation at which future ice cores can be drilled to obtain successive, high-quality paleoclimatic records in the north western Himalayas. Concordant with a recent report on the mass loss at high elevations of southern and central Tibetan Plateau inferred from ice core records 22 , there is an urgent need to recover more ice cores from the Tibetan Plateau and the surrounding region; this area is a key region in climate research and such data must be collected before disappearance of the glacial records preserved therein.

Materials and Methods
Observations. Measurements of glacier mass balance were made using the glaciological method on the northern branch of the Naimona'nyi Glacier beginning in 2004 (Fig. 1b). The number of monitoring stakes over the entire glacier surface has progressively increased from four in 2004 to 31 in 2013. The heights of the individual stakes and measurements of snow pits were manually recorded in early October of every year to derive the annual mass balance. Due to logistical issues, no fieldwork was carried out in 2005, and the stakes drilled in 2004 were missing when the sites were revisited in October 2007. The annual mass balances for nine balance years, including the average value during the 2004-06 balance year, are available for this glacier over the past decade (Fig. 2). The temperature and precipitation used to force the mass balance model was recorded at the nearest meteorological station (Burang station), which is located approximately 20 km from the Naimona'nyi Glacier (Fig. 1).
Where M(i, t) is the daily melting (mm w.e.) and A(i, t) is the daily accumulation (mm w.e.) in i th elevations; DDF snow/ice is the melting factor for snow or ice; T(i, t) and P(i, t) are the daily mean air temperature and daily precipitation, respectively, which are extrapolated from the Burang meteorological station (T cs and P cs ) using the variable yearly temperature lapse rate (γ t ) and a constant precipitation gradient (γ p ); z cs and z i are the elevation of Burang station and the different elevational bands, respectively, and T M /T p are the threshold temperatures below or above which melt or solid snowfall is assumed to be zero. The area-averaged mass balance (B n ) is calculated by considering the area weight at each elevation bin. Based on the topographic map compiled in 1976 and two Landsat satellite images taken in 1999 and 2007, the temporal change of glacier area was considered in the model by linear area interpolation 28 .
In the mass-balance model, the annual γ t was calculated using ERA Interim reanalysis temperature data (June-September, 1979-2014) from 3 × 3 grids around the Naimona'nyi Glacier 29 . The γ t for the period from 1974 to 1978 is assumed to be same as the value in 1979. The temperature bias was corrected by comparison with temperature records at AWS1 in 5500 m a.s.l. (Supplementary Fig. S5). In the model, five parameters were first calibrated and evaluated to determine the introduction of uncertainies [30][31][32] (Supplementary Table S1). A Monte Carlo simulation consisting of 1500 realizations was performed to estimate the model uncertainties 33,34 . In this experiment, a pool of 100 out of the 1500 parameter combinations yielded the RMSD less than 0.50 m w.e. between the modeled and measured mass balance on the surface of the Naimona'nyi Glacier from 2004 to 2014 (Supplementary Fig. S1 and Fig. 2). Simulations with these 100 parameter combinations were selected for this study. The 100 individidual modeling outcomes represent the possible uncertainty ranges, and the average value was adopted as the accepted value for the Naimona'nyi Glacier.