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

Meteorology

We reconstruct the meteorology at the South Col using observations from the automatic weather station (AWS) there (at 7945 m a.s.l)⁵ 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²⁹ 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⁻⁸W m⁻² K⁻⁴), 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⁵ 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³⁰. However, we recalibrate the regression using the Phortse and Basecamp AWSs because these new sites have weighing precipitation gauges protected by double alter shields⁵, 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²⁰ 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⁻³, the bottom density to 800 kg m⁻³, 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⁻²) 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⁻¹ or kg s⁻¹) 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). The resulting amplification in mean melt rate (left( {acute{A} = frac{{left( {{sum} {overline {Q_i} } } right)_{{{{mathrm{ice}}}}}}}{{left( {{sum} {overline {Q_i} } } right)_{{{{mathrm{snow}}}}}}}} right)), though, is by almost a factor of 500. To understand this result, note that the proportional increase in melt rate can be written:

where (overline {Q_{{{{mathrm{sw}}}}}}) and (overline {Q_{{{{mathrm{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 (acute{A}) is inversely proportional to the baseline melt rate in the snow simulation (left( {overline {Q_{{{{mathrm{sw}}}}}} – overline {Q_{{{{mathrm{sinks}}}}}} } right)). The very low melt rate in the snow scenario (3.3 mm w.e. a⁻¹), 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⁻³), U is the wind speed (m s⁻¹), 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₀ is the saturation vapor pressure at the melting point, L is the latent heat of sublimation (2.83 × 10⁶ J kg⁻¹), and R_v is the gas constant for moist air (461 J K⁻¹); 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 (left( {frac{{{mathrm{d}}x}}{{{mathrm{d}}t}}} right)) 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 (1950–2019), and ∆t to the number of seconds in 1 year (3.2 × 10⁷ s). We prescribe the turbulent exchange coefficient (C_e) using the output from the COSIPY snow simulation, dividing the simulated sublimation by (rho ;U(e_{{{mathrm{s}}}} – e_{{{mathrm{a}}}}Upsilon )frac{varepsilon }{{P_{{{mathrm{a}}}}}}Delta t) (see Eq. (12)).

Inserting these values into Eqs. (14) (16) yields:

We then estimate the time derivatives for Eq. (13) using the Theil-Sen slope estimation, yielding:

Summing these terms indicates a theoretical trend (left( {frac{{{mathrm{d}}S}}{{{mathrm{d}}t}}} right)) of 0.27 mm w.e. a⁻², in reasonably close agreement with the 0.22 mm w.e. a⁻² 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):

In which we have assumed melting conditions (T_s = 273.15, e₀ = 611.3 Pa; L_v is the latent heat of vaporization [2.5 × 10⁵ J kg⁻¹], and the longwave thermal radiation emitted by the snow surface is 312.5 W m⁻²); c_p is the specific heat content of the air (1004.7 J kg⁻¹ K⁻¹). It has been concluded³¹ 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₁ = −17 (−168) W m⁻², c₂ = 0.73 (2.12) W m⁻² %⁻¹ and c₃ = 0.57 (0.84) (dimensionless) for clear (cloudy) conditions.

The derivative of these fluxes with respect to air temperature is then:

Note that all terms in Eqs. (23)–(25) are positive, outlining the physical basis for why melt rates should increase with rising air temperature³².

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^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}) 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.

In summary, then, theory indicates that rising air temperatures should accelerate both sublimation and melt rates (Eqs. (15) and (24)). However, increases in wind speed and relative humidity will have opposite effects. Due to the persistence of freezing air temperatures during surface melt events, faster winds act to enhance sublimation (Eq. (14)) but reduce melting (Eqs. (28) and (29)), whereas increasing relative humidity amplifies melting (Eqs. (26) and (27)) but dampens sublimation (Eq. (16)).

Source: Resources - nature.com