Mt. Everest’s highest glacier is a sentinel for accelerating ice loss

Mountain glacier systems are decreasing in volume worldwide yet relatively little is known about their upper reaches (>5000 m). Here we show, based on the world’s highest ice core and highest automatic weather stations, the significant and increasing role that melting and sublimation have on the mass loss of even Mt. Everest’s highest glacier (South Col Glacier, 8020 m). Estimated contemporary thinning rates approaching ~2 m a−1 water equivalent (w.e.) indicate several decades of accumulation may be lost on an annual basis now that glacier ice has been exposed. These results identify extreme sensitivity to glacier surface type for high altitude Himalayan ice masses and forewarn of rapidly emerging impacts as Mt. Everest’s highest glacier appears destined for rapid retreat.


INTRODUCTION
The worldwide retreat of mountain glaciers in recent decades is well documented 1 and the resulting loss of water storage capacity for agriculture, hydropower, and both human and ecosystem consumption have significant impacts on the 250 million people living near mountain glaciers 2 . In addition, more than 1.6 billion people are recipients of water from mountain regions and 50% of Earth's biodiversity centers are in mountain regions 2 . However, the upper reaches (>5000 masl) of mountain landscapes have received relatively little scientific attention leaving gaps in knowledge concerning the key drivers influencing atmospheric circulation, changes in snow and ice extent over time, and climate model verification that together will decrease uncertainty in the climate change projections needed to plan a sustainable future 3 . Within the Hindu Kush Himalaya (HKH) lie Earth's highest mountains including the foremost Mt. Everest (Sagarmatha, Qomolangma, 8848.86 m). Despite its iconic status and having been climbed >7000 times, Mt. Everest remains in its highest reaches poorly understood in terms of weather, climate, and glacier health. To fill some of these knowledge gaps National Geographic and Rolex's Perpetual Planet Everest Expedition mounted in April/May 2019 the most comprehensive scientific investigation of the Nepalese side of Mt. Everest thus far undertaken, including studies in biology, geology, glaciology, meteorology, and mapping 4 . This expedition resulted in (Fig. 1): the world's highest ice core -recovered from Mt. Everest's highest glacier (27.977211, 86.929861; 8020 m, South Col Glacier (SCG), mean annual air temperature 1991-2020 = −22.6°C); and the world's two highest automatic weather stations (AWSs) located along the southern climbing route of Mt. Everest at 7945 m on South Col and 8430 m on the Balcony 5 . In the following, we investigate the timing and cause of the significant SCG mass loss documented during our investigation and the implications for high mountain glacier systems.

Estimated SCG thickness decrease and timing
Moraines surrounding SCG attest to the past larger extent of the glacier (Supplementary Fig. 1) assumed to mark its Little Ice Age position (1300s to late 1800s 6 ). Except for seasonal snow and a perennial snow/firn apron on the flanks of Mt. Everest, which form the upper reaches of this southerly facing glacier, SCG's surface is primarily exposed ice. Logistic constraints (oxygen and weather) restricted the time available to drill and package the ice core to two hours resulting in a surface to 10 m depth retrieval of ice (average density~0.89 gm/cm 3 ) from what is estimated to be its current 30-50 m thickness. Aerosol-based micro-radiocarbon 7 dating of the upper 10-69 cm of the SCG ice core reveals an age of 1966 ± 179 years ago. Identification of annual layers in the SCG ice core, using ultra-high resolution (153 μ) sampling 8,9 of magnesium (winter/spring maxima), as previously demonstrated for the East Rongbuk Glacier ice core 10 (ERG; 6518 m,~5 km north of SCG), reveals a net annual layer thickness of~27 mm w.e. a −1 for the SCG ice core (Supplementary Fig. 2). If we assume this is representative of the last 2000 years, multiplying by the near-surface age this yields an estimated net thinning of SCG of~55 m w.e. ERG ice core annual layer counted dating reveals an age of~500 years ago at~80 m w.e. depth 11 and a model estimated depth of~2000 years ago at 108 m suggesting~50% less ice thickness for the same age at SCG. In addition, estimated modern precipitation for SCG (Methods section) is~50% of the~480 mm w.e. a −1 derived for ERG modern precipitation 11 . While a definitive ice thickness loss cannot be determined we suggest that by comparison with the ERG ice core our SCG ice thinning estimate is reasonable.
To determine the timing of the ice loss on SCG we refer to written records, photogrammetry, ice core and meteorological records, and climate reanalysis for the region. A synthesis of written and photographic records of retreat/advance states of 112 Himalayan glaciers covering the period 1812-1965 demonstrate that some glacier recession started as early as 1860, but with considerable temporal and spatial variability 12 . More recent syntheses suggest that glacier mass loss has been the overall trend since the 1950s with the greatest rate of glacier area loss since~2000 13 . For the south side of Mt. Everest, combined photogrammetric and satellite imagery reveals thinning greater than 100 m up to 5700 m since 1962, with a near doubling in rate since 2009 14 . Climate reanalysis data show June/September freezing level heights rising~7 m per year since 2005 15 and an AWS operated May/July 2005 near the ERG ice core site reveals the increased ablation effect of cloud cover 16 . Ice cores collected in 1980 show the effects of melting at 6100 and 6400 m on the Khumbu Glacier 17 . Air content used as a proxy for summer temperature measured in the ERG ice core indicates that the last 100 years have been the warmest in the last 2000 years 18 .
While SCG is at a significantly higher elevation than sites previously noted we suggest that its thinning (~55 m w.e.) likely follows the overall ice mass loss trajectory in the region by starting in the mid-late 1800s, increasing as of the 1950s, and most marked since~2000. The newly extended fifth major global reanalysis of the European Centre for Medium-Range Weather Forecasts (ERA5) that extends back to 1950 19 suggests large-scale warming of annual surface temperature over Asia starting in the late 1990s ( Supplementary Fig. 3) with the most significant warming over the Tibetan Plateau and northern slope of the Himalayas for the period 2001-2020 minus 1979-2000 approaching~2°C in winter and a smaller increase in summer (Fig. 2). ERA5 reanalysis focused over Mt. Everest reveals a sharp demarcation of warming below 150 hPa as of the late 1950s with the most intense warming starting in the late 1990s ( Fig. 3 and Supplementary Fig. 4).

Modeling SCG snow and ice loss
To explore whether recent meteorological trends can be reconciled with SCG ice loss we perform experiments using the COSIPY model (COupled Snowpack and Ice surface energy and MAss balance model) 20 applied to downscaled ERA5 data ( Fig. 4; Methods section). First, we assess the extent to which changes in climate (Supplementary Fig. 5; air temperature, relative humidity, wind speed, and precipitation, short-and longwave radiation) explain the glacier mass loss. We do so by simulating the energy and mass fluxes for a prescribed snowpack (1950-2019; Methods section), assumed to be in surface mass balance (SMB) for the first decade of the experiment (1950)(1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958)(1959). During this period of equilibrium, we find sublimation to be by far the main process of mass loss (mean loss = 77 mm w.e. a −1 ), over 35 times the (minimal) surface melt (2 mm w.e. a −1 ), in agreement with energy balance simulations from relatively high altitude, low-latitude glaciers elsewhere 21,22 . We find no significant trend in snowfall over the full simulation period (1950-2019), but upward in sublimation of 0.22 [95% confidence interval: 0.10-0.38] mm w.e. a −2 . Whilst surface melt remains negligible throughout the snow simulation, it does exhibit a significant upward trend of 0.04 [0.01-0.07] mm w.e. a −2 , and it combines with the uptick in sublimation to drive a significant decline in the SMB of −0.33 [−0.51 to −0.14] mm w.e. a −2 ( Supplementary Fig. 6). Both physical and empirical lines of evidence (Methods: Sublimation & Melt Sensitivity to Climate Forcing) identify the increasing air temperature as playing a key role in explaining the sublimation trend. As the surface warms, saturation vapor pressure grows at the Clausius Clapeyron rate, generating an ever-greater moisture gradient from the surface into the atmosphere for fixed relative humidity. At the SCG declining relative humidity amplifies this gradient further, with strengthening winds also helping to increase sublimation by enhancing turbulent heat exchange. This experiment demonstrates that even if SCG snowpack was in SMB in the middle of the 20th century, changing climate likely drove considerable thinning since then, with a cumulative SMB of −1530 mm w.e. by the end of 2019.
A second COSIPY experiment involved running with no initial snowpack, enabling estimation of the potential rate of mass loss once ice is exposed at SCG. The results indicate a mean 1950-2019 SMB of −1929 mm w.e. a −1 , with a significant downward trend of −3.7 [−5.9 to −1.5] mm w.e. a −2 . Mean annual ablation (the sum of sublimation and melt) is over 20 times greater in the ice simulation (1964 mm w.e. a −1 ) compared to the snow simulation (96 mm w.e. a −1 ), due to sublimation increasing by a factor of 4.9 (from 93 to 456 mm w.e. a −1 ) and the initiation of widespread melting as a tipping point is crossed, rising from negligible amounts when the SCG was snow covered (3.3 mm w.e. a −1 ) to significant meltwater generation for the ice surface (1508 mm w.e. a −1 ). Once the snowpack is replaced by ice that reflects less than half the insolation, the SCG surface is more often raised to 0°C. After being warmed to the melting point, additional energy gains from solar heating are less easily lost to the turbulent and longwave fluxes (which would otherwise amplify as the glacier surface warms relative to the atmosphere). Energy must therefore be dissipated through melting once the glacier surface reaches 0°C, much like the spillway that diverts runoff once a reservoir has filled ( Supplementary Fig. 10). The high receipts of insolation at the South Col 5 help explain this acute sensitivity of the surface energy balance to surface reflectivity, consistent with other lowlatitude glaciers at relatively high elevation [23][24][25] .
The COSIPY experiments therefore outline a plausible mechanism for dramatic mass loss at SCG. First, once glacier ice is regularly exposed, the results suggest that the~55 m w.e. thinning could occur in~25 years, over 80 times faster than the~2000 years it took to form the ice now exposed at the surface of SCG. Second, the COSIPY snow simulation indicates that the transition from a permanent snow/firn surface to majority ice-cover could have been triggered by climate change since 1950, with sublimation enhanced by rising air temperatures playing the critical role. Although the total mass of snow lost at the SCG is unknown and would be sensitive to the onset of negative mass balances (which could pre-date 1950) 12 , and the physical properties of the snowpack (e.g., its temperature and capacity for refreezing of meltwater) 20 , the modeled thinning provides a plausible mechanism to expose glacier ice. We also note that whilst the magnitude of the snow loss is sensitive to parameter uncertainty in our simulations, the interpretation of thinning is not ( Supplementary  Information). Similarly, whilst the ability of ERA5 to capture historic climate variability for the high Himalaya is unknown, warming of the upper troposphere is a robust and widespread feature of anthropogenic climate change 26 and physical considerations of the surface energy balance identify that enhanced sublimation should follow from this heating (Methods section). We therefore highlight that the increasing ablation of snowpack identified here  19 . Maps generated using Climate Reanalyzer 33 . is, in the absence of other changes in meteorology, the expected response to greenhouse gas forcing. The surface mass balance in the ice simulations is also impacted by parameter uncertainty, but our interpretation of rapid loss (at least 1200 mm w.e. a −1 ) is robust across all scenarios considered.

DISCUSSION
Beyond illustrating that climate change may have driven dramatic mass loss at SCG in recent decades, our simulations also highlight mechanisms that may be of much broader significance for glacier retreat across the Himalaya. First, we find that the region's extreme insolation means that ablation can accelerate by over a factor of 20 as snow cover gives way to glacier ice. This is particularly critical for glaciers like SCG that have relatively small snow accumulation rates. Second, climate trends have led to mainly sublimation-driven thinning of high-altitude snowpack since 1950, making such transitions more likely. This mechanism is significant because our results indicate that whilst warming air temperatures did the most to enhance sublimation, declining relative humidity and strengthening winds also played a role, which should have acted to suppress melting (Methods section). This contrasting behavior underlines the theoretical shortcoming of widely used empirical glacier ablation models 27 and highlights the need to differentiate between these processes when projecting future glacier mass balance. The latter is essential to anticipate the type of highly non-linear mass loss identified here. In identifying sublimation as a major control on the mass balance, our results add to the consensus regarding the high sensitivity of low-latitude ice masses to moisture variability 28 . We also note that the acceleration of upper-atmosphere winds we report may have enhanced snow ablation at the SCG by increasing physical deflation, a process not considered in our simulations, yet also capable of triggering rapid glacier loss if it leads to more frequent ice exposure.
Warming will also change the experience associated with Mt. Everest ascents as loss of high elevation snow and ice cover continues to thin exposing bedrock; warm thicker air increases oxygen availability 29 ; ice block movement in the Khumbu icefall and avalanches becomes even more dynamic; and glacier melt destabilizes the Khumbu base camp that is home to~1000 climbers and logistics teams during the climbing season. Climate predictions for the Himalaya all suggest continued warming and continued glacier mass loss 13 . We find the transition from snow covered to exposed ice can change the state of the glacier from one of equilibrium to one of extremely rapid mass loss. At an estimated thinning rate approaching 2000 mm a −1 even glaciers such as SCG that are above 8000 m may disappear by mid-century. Our study points to the critical balance afforded by snow-covered surfaces and the potential for loss throughout high mountain glacier systems as snow cover is depleted by changes in sublimation and surface melt driven by climate trends. Everest's highest glacier has served as a sentinel for this delicate balance and has demonstrated that even the roof of the Earth is impacted by anthropogenic source warming.

METHODS Meteorology
We reconstruct the meteorology at the South Col using observations from the automatic weather station (AWS) there (at 7945 m a.s.l) 5 to downscale ERA5 reanalysis via a parsimonious blend of bias correction and machine learning. Initial screening indicates strong correlations between hourly ERA5 pressure level data bilinearly interpolated to South Col and air temperatures (r = 0.98), wind speed (r = 0.94), and relative humidity (r = 0.80) observed by the AWS. We therefore apply a simple empirical quantile mapping correction 29 to remove systematic bias for these variables.
Incident shortwave (SW) and longwave (LW) radiation are not available on ERA5 pressure levels, so we reconstruct them by downscaling the transmissivity (τ) and emissivity (α) of the atmosphere, defined: Where Ψ is the theoretical top of atmosphere solar radiation, σ is the Stefan Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), and T a is the 2 m air temperature (Kelvin). Observed values for τ and α are evaluated using AWS measurements of incident radiation and air temperature (using calculations of solar geometry to compute Ψ). We then train a Random Forest model (with 100 trees and a minimum leaf size of three) using Python's Scikit Learn (version 0.20.1), modeling τ and α as a function of the ERA5 predictors in Methods Table 1. SW and LW could then be computed from: Where T is the estimate from the bias-corrected ERA5 data. We calibrate the bias correction and RF models using between 5012 (wind speed) and 12,810 (air temperature) overlapping hours of AWS observations and ERA5 data (May 2019 to December 2020). We evaluate the performance using a fivefold cross-validation, with results indicating very strong agreement between the observed and downscaled meteorology: hourly Pearson correlations range from 0.83 (relative humidity) to 0.98 (air temperature), translating respectively to root mean square errors between~31 and 8% of the observed means (Supplementary Fig. 7). We also detect no sign of a seasonal dependence in the performance of bias correction and RF models (Supplementary Fig. 8). The resulting downscaled ERA5 data provide a complete annual series of hourly values for 1950-2019.
We estimate precipitation at the South Col also using ERA5. First, we linearly interpolate the reanalysis data to the location of the Phortse AWS 5 and then compute the ratio of the total observed precipitation (P o ) and ERA5 precipitation (P E ) during the overlapping period (April 2019-November 2020). We then multiply all reanalysis precipitation by this scalar to produce a corrected precipitation series (P E ') for Phortse 1950-2019: To extrapolate to the South Col, we assume that precipitation decays exponentially with increasing elevation 30 . However, we recalibrate the regression using the Phortse and Basecamp AWSs because these new sites have weighing precipitation gauges protected by double alter shields 5 , and hence are less prone to under-catch error ( Supplementary  Fig. 9). Note that, as described below, the precipitation estimate is adjusted to an "effective" flux before being used to simulate glacier mass balance changes.

Mass balance
We use the precipitation, along with the other downscaled meteorological variables, to force the COSIPY model 20 at hourly resolution for 1950-2019. First, we compute the effective precipitation (which implicitly includes the net effects of avalanching and wind transport, as well as correcting for any systematic bias in the downscaling/extrapolation method described above) required for the glacier to be in equilibrium for the period 1950-1959, iterating until the surface mass balance is zero. This is achieved when the precipitation is decreased by 65%. We then run two simulations with COSIPY. The first (referred to as the "snow" simulation) assumes a starting snowpack that is arbitrarily deep (20 m) to ensure that it remains present throughout the entire (70-year) simulation. We set the initial surface density of the snowpack to 350 kg m −3 , the bottom density to 800 kg m −3 , and linearly interpolate between. The second simulation (hereafter the "ice" simulation) uses the same effective precipitation but assumes no initial snowpack. However, snow is free to accumulate in the model in response to meteorological forcing. The algorithms and parameter values used in our application of COSIPY are outlined in Methods Table 2.

Melt
The surface melt rate depends on the surface energy balance (SEB): where Q denotes energy flux (W m −2 ) and the subscripts h, l, lw, sw, g, and r refer to the sensible, latent, net longwave radiative, net shortwave radiative, ground, and precipitation heat fluxes, respectively. The fluxes are defined as positive when directed towards the surface. The energy consumed in melting (Q m ) is also defined as positive, meaning the melt rate (M; mm w.e. s −1 or kg s −1 ) can be calculated: in which H(Q,T s ) is a Heaviside function that returns a value of one unless both the sum of the first six terms in Eq. (6) (Q m = ∑Q i , with i indexing terms Q h to Q r ) is positive and the surface temperature is also at the melting point; otherwise, it returns zero. The melt total over a period of Δt seconds can then be expressed: where P is the fraction of Δt during which melting conditions occurred, and the overbar for energy component Q i indicates the mean value calculated during melting conditions. In terms of energy components Eq. (8) is the major driver of the amplification in melt totals between the snow and ice COSIPY simulations, increasing by a factor of 4.4; the energy sinks (sum of the remaining terms) amplify by a factor of 3.6 (Methods Fig. 3 increase in melt rate can be written: where Q sw and Q sinks are the mean energy gains and losses, respectively, during melt conditions in the snow simulation, and k and j are the proportional increases in these terms when transitioning to an ice surface (4.3 and 3.6, respectively). Critically, Eq. (9) reveals that A is inversely proportional to the baseline melt rate in the snow simulation Q sw À Q sinks À Á . The very low melt rate in the snow scenario (3.3 mm w.e. a −1 ), therefore acts to amplify the numerator of in Eq. (9).

Sublimation and melt sensitivity to climate forcing
Total sublimation (S, mm w.e. or kg) can be written: where ρ is the air density (kg m −3 ), U is the wind speed (m s −1 ), C e is the turbulent exchange coefficient for moisture (dimensionless), ε is the ratio of gas constants for water vapor and air (0.622), P a is air pressure (Pa), and Υ is relative humidity (fraction). The saturation vapor pressure for the surface (e s ) and the near-surface atmosphere (e a ) are functions of the surface (T s ) and air temperature (T a ), respectively. If we assume that T s = T a (which is a reasonable simplification at the South Col where air temperature does not rise above 0°C), and use the Clausius Clapeyron equation: in which e 0 is the saturation vapor pressure at the melting point, L is the latent heat of sublimation (2.83 × 10 6 J kg −1 ), and R v is the gas constant for moist air (461 J K −1 ); Eq. (10) becomes: which can be differentiated with respect to U, T a , and Υ to explore the sensitivity of sublimation to changes in these meteorological parameters.
In turn, the contribution of temporal trends dx dt À Á in these variables to the trend sublimation can be evaluated via the chain rule: With: To evaluate Eq. (13) we compute the derivatives (Eqs. (14)-(16)) using the mean meteorology at the South Col during the ERA5 reconstruction , and Δt to the number of seconds in 1 year (3.2 × 10 7 s). We prescribe the turbulent exchange coefficient (C e ) using the output from the COSIPY snow simulation, dividing the simulated sublimation by ρ Uðe s À e a ϒÞ ε Pa Δt (see Eq. (12)  Summing these terms indicates a theoretical trend dS dt À Á of 0.27 mm w.e. a −2 , in reasonably close agreement with the 0.22 mm w.e. a −2 derived from the COSIPY snow simulation and reported in the main text. This theoretical analysis also indicates that 82% of the trend can be attributed to increasing air temperature (41%) and declining relative humidity (41%), with strengthening winds explaining the remaining 18%.
An alternative (empirical) method to estimate the sensitivity of sublimation in the COSIPY snow simulation is to use multiple linear regression: Where the slope coefficients (β x ) are linear approximations of the derivatives, relating the changes in the annual mean of the meteorological variables to the total annual sublimation. Performing the regression (Supplementary Fig. 11) lends support to the interpretation from the theoretical analysis above, attributing 49, 26, and 25% of the sublimation increase to the trends in air temperature relative humidity, and wind speed, respectively. In the main text, we highlight that sublimation and melt rates may differ in their response to climate forcing. To support this assertion, we repeat the sensitivity assessment above, evaluating the derivatives of the melt rate with respect to air temperature, wind speed, and relative humidity.
We simplify the analysis by assuming that the proportion of time that the surface is melting (P) is constant (see Eq. (8)). Although physically incomplete, we note that there is no temporal trend in P for the COSIPY snow simulation (p > 0.05 according to Seil-Then slope estimation). With this simplification, the sensitivity of the melt rate to changes in meteorological component x can then be written: Wind speed, air temperature, and relative humidity appear in the expressions for the sensible, latent, and longwave heat fluxes (Q h , Q l , and Q lw , respectively): Q l ¼ ρ U C e L v ðϒe a À 611:3Þ ε P a (20) Q lw ¼ ασT 4 a À 312:5 In which we have assumed melting conditions (T s = 273.15, e 0 = 611.3 Pa; L v is the latent heat of vaporization [2.5 × 10 5 J kg −1 ], and the longwave thermal radiation emitted by the snow surface is 312.5 W m −2 ); c p is the specific heat content of the air (1004.7 J kg −1 K −1 ). It has been concluded 31 that the incident longwave flux Q lw ↓ in the Himalaya may be estimated from Υ and T a : where the c x terms are empirically determined coefficients, whose value depends on cloudiness. Optimizing this expression for the South Col AWS, we found c 1 = −17 (−168) W m −2 , c 2 = 0.73 (2.12) W m −2 % −1 and c 3 = 0.57 (0.84) (dimensionless) for clear (cloudy) conditions. The derivative of these fluxes with respect to air temperature is then: ∂Q l ∂T a ¼ ρ U C e L v L ϒ ε e 0 e L Rv À Note that all terms in Eqs. (23)- (25) are positive, outlining the physical basis for why melt rates should increase with rising air temperature 32 .
The derivative of these fluxes with respect to relative humidity is: Because all terms are positive in Eqs. (26) and (27), increases in relative humidity also drive increases in the melt rate.
The derivatives of the sensible and latent heat fluxes with respect to wind speed are then: Because T a is always less than 273.15 K during melt events at the South Col in the COSIPY snow simulation ( Supplementary Fig. 12), e 0 e L Rv À 1 Ta þ 1 273:15 ð Þ must be less than 611.3 Pa. Hence, the last terms in Eqs. (28) and (29) are negative, and so increases in wind speed act to reduce the melt rate.

DATA AVAILABILITY
The downscaled meteorological data and COSIPY simulations are available from the corresponding authors on reasonable request.