in

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

[adace-ad id="91168"]

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 correction29 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:

$$tau = {{{mathrm{SW/}}}}{Psi}$$

(1)

$${{{mathrm{and}}}};alpha = {{{mathrm{LW/}}}}sigma T_{{{mathrm{a}}}}^4$$

(2)

Where Ψ is the theoretical top of atmosphere solar radiation, σ is the Stefan Boltzmann constant (5.67 × 10−8 W m−2 K−4), and Ta 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:

$${{{mathrm{SW}}}} = tau {Psi}$$

(3)

$${{{mathrm{LW}}}} = varepsilon T^4$$

(4)

Where T is the estimate from the bias-corrected ERA5 data.

Table 1 Predictor variables used in the machine learning downscaling.
Full size table

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 AWS5 and then compute the ratio of the total observed precipitation (Po) and ERA5 precipitation (PE) during the overlapping period (April 2019-November 2020). We then multiply all reanalysis precipitation by this scalar to produce a corrected precipitation series (PE‘) for Phortse 1950–2019:

$$P_{{{mathrm{E}}}}^prime = frac{{P_{{{mathrm{o}}}}}}{{P_{{{mathrm{E}}}}}}P_{{{mathrm{E}}}}$$

(5)

To extrapolate to the South Col, we assume that precipitation decays exponentially with increasing elevation30. However, we recalibrate the regression using the Phortse and Basecamp AWSs because these new sites have weighing precipitation gauges protected by double alter shields5, 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 model20 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.

Table 2 Parameterizations and parameter values used in the COSIPY model runs.
Full size table

Melt

The surface melt rate depends on the surface energy balance (SEB):

$$Q_{{{mathrm{h}}}} + Q_{{{mathrm{l}}}} + Q_{{{{mathrm{lw}}}}} + Q_{{{{mathrm{sw}}}}} + Q_{{{mathrm{g}}}} + Q_{{{mathrm{r}}}} – Q_{{{mathrm{m}}}} = 0$$

(6)

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 (Qm) is also defined as positive, meaning the melt rate (M; mm w.e. s−1 or kg s−1) can be calculated:

$$M = frac{1}{{L_{{{mathrm{f}}}}}}mathop {sum}Hleft( {Q_{{{mathrm{m}}}},T_{{{mathrm{s}}}}} right)Q_i$$

(7)

in which H(Q,Ts) is a Heaviside function that returns a value of one unless both the sum of the first six terms in Eq. (6) (Qm = ∑Qi, with i indexing terms Qh to Qr) 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:

$$M = Pfrac{{{Delta}t}}{{L_{{{mathrm{f}}}}}}{sum} {overline {Q_i} }$$

(8)

where P is the fraction of ∆t during which melting conditions occurred, and the overbar for energy component Qi 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:

$$acute{A} = frac{{koverline {Q_{{{{mathrm{sw}}}}}} – joverline {Q_{{{{mathrm{sinks}}}}}} }}{{overline {Q_{{{{mathrm{sw}}}}}} – overline {Q_{{{{mathrm{sinks}}}}}} }}$$

(9)

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−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:

$$S = rho ;U;C_{{{mathrm{e}}}}(e_{{{mathrm{s}}}} – e_{{{mathrm{a}}}}{Upsilon})frac{varepsilon }{{P_{{{mathrm{a}}}}}}{Delta}t$$

(10)

where ρ is the air density (kg m−3), U is the wind speed (m s−1), Ce is the turbulent exchange coefficient for moisture (dimensionless), ε is the ratio of gas constants for water vapor and air (0.622), Pa is air pressure (Pa), and Υ is relative humidity (fraction). The saturation vapor pressure for the surface (es) and the near-surface atmosphere (ea) are functions of the surface (Ts) and air temperature (Ta), respectively. If we assume that Ts = Ta (which is a reasonable simplification at the South Col where air temperature does not rise above 0 °C), and use the Clausius Clapeyron equation:

$$e_{{{mathrm{s}}}} = e_0e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}$$

(11)

in which e0 is the saturation vapor pressure at the melting point, L is the latent heat of sublimation (2.83 × 106 J kg−1), and Rv is the gas constant for moist air (461 J K−1); Eq. (10) becomes:

$$S = rho ;U;C_{{{mathrm{e}}}}(1 – {Upsilon})e_0e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}frac{varepsilon }{{P_{{{mathrm{a}}}}}}{Delta}t$$

(12)

which can be differentiated with respect to U, Ta, 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:

$$frac{{{mathrm{d}}S}}{{{mathrm{d}}t}} = frac{{partial S}}{{partial U}}frac{{{mathrm{d}}U}}{{{mathrm{d}}t}} + frac{{partial S}}{{partial Y}}frac{{{mathrm{d}}Y}}{{{mathrm{d}}t}} + frac{{partial S}}{{partial T_{{{mathrm{a}}}}}}frac{{{mathrm{d}}T_{{{mathrm{a}}}}}}{{{mathrm{d}}t}}$$

(13)

With:

$$frac{{partial S}}{{partial U}} = rho ;C_{{{mathrm{e}}}};e_0left( {1 – {Upsilon}} right)e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}frac{varepsilon }{{P_{{{mathrm{a}}}}}}{Delta}t$$

(14)

$$frac{{partial S}}{{partial T_{{{mathrm{a}}}}}} = L;C_{{{mathrm{e}}}};e_0;rho ;Uleft( {1 – {Upsilon}} right)frac{{e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}}}{{R_{{{mathrm{v}}}}P_{{{mathrm{a}}}}T_{{{mathrm{a}}}}^2}}varepsilon {Delta}t$$

(15)

$$frac{{partial S}}{{partial {Upsilon}}} = – rho ;C_{{{mathrm{e}}}};e_0;Ue^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}frac{varepsilon }{{P_{{{mathrm{a}}}}}}{Delta}t$$

(16)

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 × 107 s). We prescribe the turbulent exchange coefficient (Ce) 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:

$$frac{{partial S}}{{dU}} = 6.0;{{{mathrm{mm}}}};{{{mathrm{w}}}}{{{mathrm{.e}}}}{{{mathrm{.}}}};{{{mathrm{a}}}}^{ – 1}{{{mathrm{m}}}}^{ – 1};{{{mathrm{s}}}}^1$$

$$frac{{partial S}}{{d{Upsilon}}} = – 1.8;{{{mathrm{mm}}}};{{{mathrm{w}}}}{{{mathrm{.e}}}}{{{mathrm{.}}}};{{{mathrm{a}}}}^{ – 1}% ^{ – 1}$$

$$frac{{partial S}}{{dT_{{{mathrm{a}}}}}} = 6.7,{{{mathrm{mm}}}},{{{mathrm{w}}}}{{{mathrm{.e}}}}{{{mathrm{.}}}},{{{mathrm{a}}}}^{ – 1},^circ {{{mathrm{C}}}}^{ – 1}$$

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

$$frac{{{{d}}{U}}}{{{{d}}t}} = – 0.08;{mathrm{m}};{mathrm{s}}^{-1} ;{{{mathrm{decade}}}}^{ – 1};{{{mathrm{and}}}};frac{{partial S}}{{partial {U}}}frac{{{{d}}{U}}}{{{{d}}t}}; = ;0.05;{{{mathrm{mm}}}};{{{mathrm{w}}}}{{{mathrm{.e}}}}.;{{{mathrm{a}}}}^{ – 2}$$

$$frac{{{{d}}{Upsilon}}}{{{{d}}t}} = – 0.6% ;{{{mathrm{decade}}}}^{ – 1};{{{mathrm{and}}}};frac{{partial S}}{{partial {Upsilon}}}frac{{{{d}}{Upsilon}}}{{{{d}}t}}; = ;0.11;{{{mathrm{mm}}}};{{{mathrm{w}}}}{{{mathrm{.e}}}}.;{{{mathrm{a}}}}^{ – 2}$$

$$frac{{{{d}}T_a}}{{{{d}}t}} = 0.16;^circ {{{mathrm{C}}}};{{{mathrm{decade}}}}^{ – 1}{{{mathrm{and}}}};frac{{partial S}}{{partial T_{{{mathrm{a}}}}}}frac{{{{d}}T_{{{mathrm{a}}}}}}{{{{d}}t}}; = ;0.11;{{{mathrm{mm}}}};{{{mathrm{w}}}}{{{mathrm{.e}}}}{{{mathrm{.}}}};{{{mathrm{a}}}}^{ – 2}$$

Summing these terms indicates a theoretical trend (left( {frac{{{mathrm{d}}S}}{{{mathrm{d}}t}}} right)) 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:

$$S = alpha + beta _{{{mathrm{U}}}}U + beta _{{{mathrm{{Upsilon}}}}}{Upsilon} + beta _{T_{{{mathrm{a}}}}}T_{{{mathrm{a}}}}$$

(17)

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:

$$frac{{{mathrm{d}}M}}{{{mathrm{d}}x}} = Pfrac{{{Delta}t}}{{L_{{{mathrm{f}}}}}}{sum} {frac{{partial underline {Q_i} }}{{partial x}}}$$

(18)

Wind speed, air temperature, and relative humidity appear in the expressions for the sensible, latent, and longwave heat fluxes (Qh, Ql, and Qlw, respectively):

$$Q_h = rho ;U;c_{{{mathrm{p}}}};C_{{{mathrm{h}}}}(T_{{{mathrm{a}}}} – 273.15)$$

(19)

$$Q_l = rho ;U;C_{{{mathrm{e}}}};L_{{{mathrm{v}}}}({Upsilon}e_{{{mathrm{a}}}} – 611.3)frac{varepsilon }{{P_{{{mathrm{a}}}}}}$$

(20)

$$Q_{lw} = alpha sigma T_{{{mathrm{a}}}}^4 – 312.5$$

(21)

In which we have assumed melting conditions (Ts = 273.15, e0 = 611.3 Pa; Lv is the latent heat of vaporization [2.5 × 105 J kg−1], and the longwave thermal radiation emitted by the snow surface is 312.5 W m−2); cp is the specific heat content of the air (1004.7 J kg−1 K−1). It has been concluded31 that the incident longwave flux Qlw↓ in the Himalaya may be estimated from Υ and Ta:

$$Q_{{{{mathrm{lw}}}}} downarrow = c_1 + c_2{Upsilon} + c_3sigma T_{{{mathrm{a}}}}^4 – 312$$

(22)

where the cx terms are empirically determined coefficients, whose value depends on cloudiness. Optimizing this expression for the South Col AWS, we found c1 = −17 (−168) W m−2, c2 = 0.73 (2.12) W m−2 %−1 and c3 = 0.57 (0.84) (dimensionless) for clear (cloudy) conditions.

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

$$frac{{partial Q_{{{mathrm{h}}}}}}{{partial T_{{{mathrm{a}}}}}} = rho ;U;c_{{{mathrm{p}}}};C_{{{mathrm{h}}}}$$

(23)

$$frac{{partial Q_{{{mathrm{l}}}}}}{{partial T_{{{mathrm{a}}}}}} = rho ;U;C_{{{mathrm{e}}}};L_{{{mathrm{v}}}};L;{Upsilon};varepsilon ;e_0frac{{e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}}}{{R_{{{mathrm{v}}}}P_{{{mathrm{a}}}}T_{{{mathrm{a}}}}^2}}$$

(24)

$$frac{{partial Q_{{{{mathrm{lw}}}}}}}{{partial T_{{{mathrm{a}}}}}} = 4sigma c_3T_{{{mathrm{a}}}}^3$$

(25)

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

The derivative of these fluxes with respect to relative humidity is:

$$frac{{partial Q_{{{mathrm{l}}}}}}{{partial {Upsilon}}} = rho ;U;C_{{{mathrm{e}}}};L_{{{mathrm{v}}}};e_0e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)}frac{varepsilon }{{P_{{{mathrm{a}}}}}}$$

(26)

$$frac{{partial Q_{{{{mathrm{lw}}}}}}}{{partial {Upsilon}}} = c_2$$

(27)

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:

$$frac{{partial Q_{{{mathrm{h}}}}}}{{partial U}} = rho ;c_{{{mathrm{p}}}};C_{{{mathrm{h}}}}(T_{{{mathrm{a}}}} – 273.15)$$

(28)

$$frac{{partial Q_{{{mathrm{l}}}}}}{{partial U}} = rho ;C_{{{mathrm{e}}}};L_{{{mathrm{v}}}}frac{varepsilon }{{P_{{{mathrm{a}}}}}}(Ye_0e^{frac{L}{{R_{{{mathrm{v}}}}}}left( { – frac{1}{{T_{{{mathrm{a}}}}}} + frac{1}{{273.15}}} right)} – 611.3)$$

(29)

Because Ta 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

Reducing methane emissions at landfills

Students dive into research with the MIT Climate and Sustainability Consortium