The future of the region's glaciers and of the water resources they sustain is the focus of an intense debate and considerable research15,16,22. Our simulations reveal that changes in snowfall were the largest in glacier accumulation areas due to a high solid precipitation ratio and, therefore, higher sensitivity to total precipitation changes. The sensitivity of mass balance to precipitation changes was enhanced by the supply of mass through the headwalls of Kyzylsu Glacier (Figs. 5 and 6), and we found that simulations that do not include gravitational redistribution would result in a lower interannual variability of mass balance (Supplementary Fig. S23). Projected long-term (2081-2100) changes in temperature suggest a warming of 3 ∘C under SSP2-4.5 and 6 ∘C under SSP5-8.5 relative to 1995-201439, which will affect the phase of precipitation at high elevations46 and could prevent further periods of snowfall increase driven by increases in total precipitation. Our simulations show no change in snowfall fraction between 1999-2018 and 2018-2023, but substantial changes in snowfall fraction can be expected from March to October by the end of the 21st century, months for which the air temperature is not largely below 0 ∘C (Supplementary Figs. S36 and S39), and their consequences on glacier mass balance and runoff seasonality should be investigated. In general, a substantial effort is still needed to understand the current state of the Pamir's glaciers and their controlling factors, and this study provides an initial step into understanding the recent decline in glacier health and its causes. Longer term reconstructions of snow and glacier mass balances are needed to further understand the long-term sensitivity of these glaciers to temperature and precipitation changes, and project their future with confidence.
The Kyzylsu catchment (168 km) lies on the northwestern slopes of the Pamir mountain range in Tajikistan, with its outlet located at 2100 m a.s.l. It serves as a tributary to the Vakhsh River, downstream of which the Roghun Dam-one of the tallest in the world-is currently under construction, a project expected to substantially increase Tajikistan's hydropower energy production. The climate is continental and semi-arid, influenced by mid-latitude westerlies and local moisture recycling. Most of the precipitation falls in winter and spring, and the glaciers are considered as winter-accumulation type. The catchment ranges in elevation from 2110 m a.s.l. to 5806 m a.s.l.. The predominant land covers comprise pastures (50%), extensively utilized for grazing by local livestock, glaciers (20%), several lakes, and rocky surfaces (29%). The glaciers are predominantly debris-covered, have experienced surges in the last 30 years, and have their surfaces level with their lateral moraines in response to the limited mass losses in recent decades. We established a monitoring network (Fig. 1) in June 2021 consisting of several automatic weather stations (AWSs), one of which includes an all-type precipitation weighing gauge (referred to as pluviometer station), time-lapse cameras, lake and stream level gauges, glacier ablation measurements, and air temperature loggers. Continuous snow depth time series were retrieved at four locations, either from daily camera photos or ultrasonic depth gauges. Supplementary Table S1 summarizes the instruments' locations, variables recorded, and recording periods. The data collected from June 2021 to September 2023 were used to downscale and bias-correct ERA5-Land reanalysis and evaluate the model's performance.
In-situ measurements of near-surface meteorology were conducted to better characterize the climate of this high-elevation site and produce the time series needed for the statistical downscaling and bias correction of the ERA5-Land reanalysis. The air temperature was measured at various locations in the catchment (Fig. 1a) from 2021 to 2023 at hourly intervals, which we used to derive hourly air temperature lapse rates, later used in the downscaling of ERA5-Land reanalysis. The derived lapse rates show clear seasonal and diurnal cycles, consistent from one year to the other (Supplementary Fig. S2). One of the AWSs (pluviometer station) is equipped with an OTT all-weather precipitation gauge and recorded precipitation continuously from September 2021 to September 2023. A windshield was only added to the pluviometer in September 2023, such that the precipitation data used in this study was recorded without a windshield. To correct for undercatch, we used the empirical relationship established by Masuda et al., which is given below:
Where U is wind speed (m/s) at the height of the rain gauge and m is an empirical correction coefficient which depends on the precipitation type and the type of rain gauge. The empirical coefficients of this relationship were derived in Japan, with different climatic characteristics than our Central Asian study site, so we estimated the correction factor for solid precipitation using cumulative snow water equivalent derived from lake pressure measurements conducted from September 2021 to June 2022. Maidakul Lake is located 500 m away from the pluviometer station and freezes over during the winter. We followed the approach described in Pritchard et al., to convert increases in water pressure into snow water equivalent of snowfall (in mm w.e.). We then calibrated the pluviometer snowfall correction coefficient (m) to match the total accumulated snowfall amount estimated between November 2021 and February 2022, a period during which the lake surface was frozen, and the air temperature was well below 0 C (Supplementary Fig. S3). We obtained a correction coefficient for solid precipitation of 0.30, which is close to the empirical value given (0.34) in the original study. We further validated the choice of the precipitation undercatch correction coefficient by conducting snow depth simulations using the different options of precipitation undercatch correction (m = 0, m = 0.30, m = 0.34) (Supplementary Fig. S3). Shortwave radiation (incoming and outgoing), wind speed, air pressure, and relative humidity were recorded at the pluviometer station and at the on-glacier AWS, while the latter also recorded longwave radiation (incoming and outgoing).
Snow depths were continuously monitored at the two AWSs using ultrasonic depth gauges and at two additional locations using time-lapse cameras with graduated vertical stakes set in their field of view. For the camera-derived snow depth, the snow height was manually read from daily pictures and an uncertainty of +/-5 cm was attributed to each snow depth value. The measured snow depths were used to evaluate the model performance and to assess the effect of the precipitation undercatch correction (Supplementary Figs. S3, S27).
We derived snow cover maps from MODIS, Landsat 5/7/8/9, and Sentinel-2 to investigate the spatiotemporal changes of the catchment snow cover and evaluate our model performance against observed fraction snow cover and snow line altitude. We used the MODIS daily product MOD10A1 version 6.1 and identified snow-covered areas using a Normalized Difference Snow Index (NDSI) threshold of 0.40, but also with thresholds of 0.25 and 0.45 to quantify the uncertainty. We discarded all MODIS scenes having a cloud cover fraction greater than 0.10 to reduce the risks of cloud/snow misclassification. While the temporal resolution of MODIS is very high (daily), its spatial resolution (500 m) is rather coarse relative to the scale of our study domain, especially in the upper areas of the domain which have high spatial variability of slope, aspect, and surface conditions. We therefore also used all available cloud-free Landsat-5/7/8/9 and Sentinel-2 scenes to derive snow cover maps of higher spatial resolution (30 m). Outside of clean-ice glacier areas, we converted multi-spectral surface reflectance into a binary snow presence indication. This method, based on NDSI and reflectance in the red band, does not perform well for discriminating bare ice from snow, such that over clean-ice areas we used surface albedo derived from the surface reflectance to distinguish snow and bare ice. The critical albedo threshold was determined for each scene using the Otsu algorithm, bounded by 0.25 and 0.55 and assuming a bi-modal distribution of albedo values on glacier areas. Examples of the snow cover mapping are provided in the Supplementary Material (Supplementary Figs. S11, 12).
Digital Elevation Models (DEMs) derived from high-resolution (<5 m) optical stereo images (e.g. Pléiades) can be used to derive high-resolution snow depth maps over entire catchments. We acquired two Pléiades stereo images over the Kyzylsu catchment and derived DEMs (2 m spatial resolution) using the Ames Stereo Pipeline. The first scene was acquired on September 24th 2022 during snow-free conditions, and the second one was acquired on May 23rd 2023 when a large proportion of the catchment was snow-covered. We performed co-registration and corrected for tilt and along-track undulations (Supplementary Fig. S14). Co-registration and DEM differencing were performed using the open-source Python package xDEM. We used the surface elevation change to identify the location of avalanche deposits and qualitatively evaluate the simulated gravitational redistribution of snow.
The Tethys-Chloris model is an ecohydrological model simulating the interplay between energy, water, and vegetation dynamics on the land surface. It has recently been updated to include energy-balance calculations over snow, clean-ice, and debris-covered glaciers and simulate the water balance of high-elevation catchments. The model explicitly simulates lateral water transfer, surface runoff, and subsurface flow, capturing the complex dynamics of soil moisture and its interactions with vegetation. It uses a single prognostic surface temperature to calculate energy fluxes and determines aerodynamic resistances for turbulent heat and vapor fluxes using the Monin-Obukhov Similarity Theory. The partitioning of precipitation into rain, snow and sleet is based on the wet bulb temperature, the snow albedo dynamics are simulated using the parameterization of Ding et al. and the model also includes schemes for snow aging, snow settling and surface sublimation. It uses a 2-layer snowpack model, featuring a 6-mm thick surface skin layer that enables energy exchange with the atmosphere and heat transfer within the snowpack. While heat transfer is calculated across multiple layers, the model simplifies computations by using single prognostic variables for snow density and water content. Liquid water can infiltrate and be stored in the snowpack, within a volume given as a percentage of the total volume, and can refreeze depending on the latent heat flux. Snow is converted to ice at a constant rate once the snowpack reaches a threshold density. The model does not include a firn layer, but this should not affect our main results (cf. Supplementary Section 4.5 in Supplementary Material). The model uses the SnowSlide parametrization to route avalanches as a function of slope and maximum snow holding depth. Glacier flow was not accounted for in this study, as glacier areas and glacier surface elevation did not change substantially through that period (Supplementary Fig. S19). The surface elevation was therefore kept constant for the downscaling of the meteorological forcing, and the ice thickness, based on the ice thickness consensus estimate, was increased artificially to prevent the disappearance of the ice column in the ablation zone. We use the globally available AW3D Standard DEM resampled to 100m as a reference surface elevation to downscale the meteorological forcing, calculate runoff routine, delineate the catchment outline, and avalanche flow paths. The model is run at an hourly time-step from 1 October 1999 to 30 September 2023. In this study, we define the hydrological year as starting on October 1st and ending on September 30th and use this as a basis for our analyses.
ERA5-Land variables are statistically downscaled from 9-km to 100-m resolution on an hourly timescale. This process involves combining spatial interpolation methods with vertical gradients for each hourly timestep. To downscale air temperature, we use monthly-hourly mean lapse rates derived from the station data described above. Due to the absence of spatially distributed measurements, no precipitation gradients were prescribed, keeping only the ERA5-Land native horizontal variability which shows larger precipitation amounts for higher elevation pixels (Supplementary Fig. S9). Wind speed was downscaled by leveraging the DEM of the catchment to identify areas that are either exposed to or sheltered from synoptic wind gradients. Incoming shortwave and longwave radiation are downscaled using a vertical gradient approach found in the litterature. Bias correction is performed monthly with the empirical quantile mapping (EQM) method, using all overlapping periods of station and reanalysis data to capture seasonal variability effectively. For precipitation, the bias correction is performed on daily sums, using the observed precipitation corrected for undercatch. Comparisons of downscaled meteorological variables against in-situ observations are shown in the Supplementary Material (Supplementary Figs. S4-7).
Glacier outlines and debris-covered areas were based on the Randolph Glacier Inventory version 6.0 and manually adjusted where needed using the September 2022 Pléiades ortho-image. Distributed debris thicknesses were derived from a relationship between remotely sensed land-surface temperature and in-situ measurements of debris thickness conducted on Kyzylsu Glacier. Land cover was set up based on the PROBA-V land cover and plant species (macro-types) were selected based on field observations and literature. The SOILGrids soil product was used for soil thickness and soil composition. A visualization of the model spatial set-up is provided in the Supplementary Material (Supplementary Figs. S15, 16). Initial snow cover and snow albedo were based on the Landsat-7 scene acquired on 16 October 1999. Snow cover and snow surface albedo in shaded areas were set to the non-shaded area average of the corresponding elevation band. In the absence of snow depth information, we set the initial snow depth to 30 centimeters (Supplementary Fig. S18).
The model parameters were not calibrated through automatic calibration of multiple parameters against spatially integrated variables (e.g. discharge) but rather set from remote sensing data, in-situ observation, or literature. In our simulations, bare ice albedo is kept constant in time, but its spatial heterogeneity is accounted for by setting its value to the mean glacier albedo observed per 100-m elevation band for the 10 least snow-covered Landsat and Sentinel-2 scenes (Supplementary Fig. S17). While our model does not have a firn layer, using an altitudinally varying bare-ice albedo allows us to indirectly account for the fact that firn usually has a higher albedo than bare ice. The thermal conductivity and surface roughness of the supraglacial debris layer were optimized through energy balance simulations at the on-glacier AWS, where meteorological data and concurrent measurements of ablation were available. Phenological parameters are the same as those used in Fugger et al.. We tested three sets of parameters for the avalanche routine (a = 0.10 and C = 145, a = 0.12 and C = 145, a = 0.14 and C = 145) and selected the couple leading to the best agreement with remote-sensing observations of avalanche deposit locations (a = 0.12 and C = 145). We evaluated the model against measured surface albedo (Supplementary Fig. S24), glacier ablation stakes (Supplementary Fig. S25), glacier surface elevation change (Supplementary Fig. S26), observed snow depth (Supplementary Fig. S27), proglacial stream water level (Supplementary Fig. S28), remotely sensed snow cover (Supplementary Fig. S29-32) and glacier-wide geodetic mass balance for the period 2000-2019 (Fig. 5). We evaluate the gravitational snow redistribution against avalanche deposit outlines detected semi-automatically using all available Sentinel-1 scenes for the 2017-2023 period and against avalanche deposits visible in the spatially distributed snow depth derived from the high-resolution Pléiades DEM differencing (Supplementary Figs. S21, 22). The performance of the model is summarized in Supplementary Table S2, and was satisfactory enough to not require further adjustment of the meteorological forcing or model parameters.
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.