Hydrological regimes and evaporative flux partitioning at the climatic ends of high mountain Asia

To simulate the coupled dynamics of water, snow and glaciers, energy, vegetation physiology and carbon in the three catchments, the Tethys-Chloris model (T&C; Fatichi et al 2012a, 2012b, Manoli et al 2018, Mastrotheodoros et al 2019, Botter et al 2021, Paschalis et al 2022) (figure S3), was applied to the period October 2010 to September 2022. The period October 2010 to September 2012 was considered a spin-up period and was not included into the analysis of model outputs. For each grid cell, the energy balance is solved for the surface temperature of any type of land cover, including vegetated and bare land, rock, surface water, snow, glaciers and supraglacial debris, and includes the computation of net radiation, sensible and latent heat, heat consumed by photosynthesis and heat exchange with the subsurface. The model also resolves infiltration and exfiltration, deep percolation, lateral surface and subsurface water flows, soil moisture dynamics in the unsaturated zone including root water uptake at different depths at the grid level. Glaciers are represented as a single-layered ice-pack including a retention capacity for water. Glacier outlines and debris covered areas were manually corrected, based on the Randolph glacier inventory (RGI-6) and the debris thicknesses product was prepared combining remote sensing data and in-situ observations (supplementary section 4). The snow-, ice and ice-under-debris mass balance models are described in detail in Fugger et al (2022) and Fyffe et al (2021). A more detailed description of T&C is provided in Fatichi et al (2012a, 2012b, 2020), and a setup in the upper part of Langtang and additional model details are described in Buri et al (2023). Advances implemented for this study include a 2-layer snowpack model, where the surface energy balance is solved at a skin layer, and a parameterization for glacier geometry changes as detailed in supplementary section 5. The sensitivity of the energy and mass fluxes involved in the computation of the glacier and snow mass balances to model parameters and meteorological forcing were presented in Fyffe et al (2021) and Fugger et al (2022).

The model parameters were not calibrated in the traditional sense, e.g. through automatic calibration of multiple parameters against spatially integrated variables like runoff. Instead they were estimated from the literature and using expert knowledge (e.g. Fatichi et al 2016). Parameters requiring optimization and manual adjustments from the T&C standard setup were related to the thermal properties of the supraglacial debris, the glacier geometry updates, bare-ice albedo, the vegetation phenology and to high-elevation precipitation at one of the sites (supplementary section 5).

The model was forced using downscaled hourly air temperature (Ta), relative humidity (RH), total precipitation (P), incoming shortwave radiation (S_in), incoming longwave radiation (L_in), wind speed (W_s) and air pressure (Pa) data from the ERA5-Land 9 km reanalysis (Muñoz-Sabater et al 2021). Statistical downscaling to 100 m spatial resolution was performed following the approach of Machguth et al (2009), which combines simple parameterisations with spatial interpolation (supplementary section 2, figures S1 and S2). Bias correction of forcing variables was performed using empirical quantile mapping against local station data (supplementary section 3). Independent stations were used, where possible, to assess the efficacy of the bias-correction in space and time (figures S4–S21). Despite minor discrepancies, the forcing data are deemed acceptable for the representation of spatial and temporal variability in the catchment-wide meteorology, as indicated by low biases, RMSE values and high correlation coefficients against observations (figures S4–S21, table S3).

Model performance was evaluated using a combination of in-situ and remote sensing data. Snowlines were derived from the MODIS Terra snow cover V6.1 product (Fugger 2018, Hall and Riggs 2021) and compared to the modelled snowlines. Snow cover was additionally evaluated scene-by-scene against Sentinel-2/Landsat-8 and Dice coefficients (Dice 1945) were calculated. The modelled surface mass balance of glaciers was evaluated against two references: the geodetic mass balance results of Hugonnet et al (2021); and altitudinal values derived following the approach of Miles et al (2021) that uses the thinning data of Hugonnet et al (2021) and the surface velocities of Millan et al (2022) at Langtang. At Kyzylsu and 24 K the velocity and geodetic mass balance data quality prevented the second approach. Modelled snow depth (all sites) and ice melt (Kyzylsu, Langtang) were evaluated using automatic measurements of surface height. Catchment mean leaf area index (LAI) was compared against MODIS (MCD15A3H v.6, Myneni et al 2015) and VIIRRS (VNP15A2H v.1; Myneni and Knyazkhin 2018) LAI products. Modelled discharge was evaluated against data from gauging stations in proglacial streams of Kyzylsu, Langtang and 24 K glaciers.

Modelled snowlines correspond well with remote sensing observations are typically within the uncertainty ranges of the observed snowlines, with biases (R² values) of −52.2 m (0.74), −16.6 m (0.43), −194.4 m (0.41) for Kyzylsu, Langtang and 24 K, respectively (figure S22).

Modelled snow cover agreed well with Landsat-8/Sentinel-2 scenes in a scene-by-scene comparison, and average Dice coefficients for each site were 0.66, 0.66, 0.76 (figure S23), with higher performance during periods of more extensive snow cover and lower performance during patchy snow cover.

Compared to observations, the timing of local snow accumulation and glacier melt was very well reproduced, while some under- and overestimations of the snowpack height up to 45 cm occurred at 24 K and Kyzylsu (figure S24). These local discrepancies are likely related to the absence of a wind-driven snow redistribution scheme in the model.

The model agrees very well with the remote-sensing-based glacier-wide and altitudinal mass balance values at Langtang (0.65 m w.e.yr⁻¹ vs 0.62 ± 0.13 w.e.yr⁻¹), to within observational uncertainty at 24 K (−0.42 m w.e. yr⁻¹ vs. 0.25 ± 2.07 m w.e. yr⁻¹), and with a slight negative overestimation at the Kyzylsu glaciers (−0.46 m w.e. yr⁻¹ vs. 0.18 ± 0.22 m w.e. yr⁻¹) (figures S25–S27).

LAI values show overall good agreement with observations in terms of magnitude (figure S28), while there are mismatches in terms of timing at Kyzylsu (RMSE: 0.52) and Langtang (0.71) and peak LAI magnitudes at 24 K (0.36).

Both the magnitude and timing of streamflow variations are well reproduced at the pro-glacial stream gauges, with NSE values of 0.5 at Kyzylsu, 0.55 at 24 K, and a somewhat weaker agreement at Langtang with 0.27 (figure S29).

The annual water throughput of the catchments, which is defined as the total water added to or removed from the water balance within a year in either liquid, solid or gaseous form, increases from west to east (figures 2 and 5). The annual water throughput is highest in 24 K, followed by Langtang and Kyzylsu: Langtang and 24 K respectively receive around 1.8 times (1742 mm yr⁻¹) and 3.7 times (3523 mm yr⁻¹) more precipitation than Kyzylsu (961 mm yr⁻¹) (figures 2 and 5). The partitioning of inputs into rain and snowfall is controlled by the precipitation seasonality: with a winter-accumulation regime, 70.6% of annual precipitation falls as snow at Kyzylsu, while those portions are smaller in Langtang (42.7%) and 24 K (39.7%) (figure 1(f)), where the bulk of the precipitation falls in the summer half-year.

Zoom In Zoom Out Reset image size

Across the common elevation range of the three catchments (2500–5700 m asl), snowmelt contributes 67% of annual runoff in Kyzylsu, 35% in Langtang, and 41% in 24 K (figure 3(a)). Once most of the snow has melted out in Kyzylsu between the middle and end of August, ice melt becomes the dominant runoff contribution with 55% in September (annually 24%). In Langtang, which has a similar proportion of glacier area, the ice melt contribution peaks in post-monsoon with 21.7% in October, while later during summer, it is dwarfed by monsoonal rainfall so that the ice melt contribution is only 14% annually. At 24 K, which hosts around half the proportion of glacier area compared to the other sites, but gets substantially higher liquid precipitation inputs, ice melt contributes only 8% to the runoff in August (3% annually) (figure 3(a)). Note that snowmelt above glacier surfaces is counted towards the catchment-wide snowmelt rather than to the ice melt.

Zoom In Zoom Out Reset image size

Rainfall runoff contributes the largest share (56%) to total runoff in 24 K, followed by Langtang (51%), and Kyzylsu (9%) (figure 3(a)). During the period March to August rainfall is transformed into runoff in Kyzylsu, while from September to October the majority of rainfall is evaporated (figure 3(a)).

In all three catchments, catchment-wide, modelled glacier mass balances were negative over the study period with −0.46, −0.65 and −0.42 m w.e. yr⁻¹ in Kyzylsu, Langtang and 24 K respectively (figures S25-S27). The imbalance melt, defined as the glacier melt not replenished by new ice formation, accounted for 13, 9 and 2% of the mass inputs to the water balance in Kyzylsu, Langtang and 24 K, respectively (figure 2). In all catchments, sub-debris ice melt is larger than bare-ice melt, with 66.8, 62.5, 57.9% for Kyzylsu, Langtang and 24 K, respectively (figure 2). The total E (combined ET and S) accounts for 28, 19 and 13% of the catchment water output (figure 2(a)) of Kyzylsu, Langtang and 24 K, respectively. ET alone corresponds to 76% of the liquid precipitation falling in Kyzylsu, 28% in Langtang, and 19% in 24 K. The relative importance of glacier melt and evaporative fluxes therefore decreases with an increasing water throughput from west to east. A spatial distribution of annual melt, evapotranspiration and sublimation is shown in figure S30.

Across the common elevation range, the evaporative fluxes follow the seasonality of net radiation (Rn) that is available for the turbulent heat transfer in the case of ET and additionally the snow covered area and temperature regime in the case of S (figure 3(b)). Soil water limitations do not affect the ET rates in either of the catchments during average years, since the latent energy (LE) and ET follow the Rn regime throughout the season, and soil water limitations do not occur, as indicated by the soil moisture seasonality (figure 3(b)). A high availability of intercepted water on leafs and rock surfaces however boosts evaporation rates in the wet climate of 24 K (figure 3(b)). Evaporation from supraglacial debris is a small flux for the catchment (not shown), but helps insulate the glacier ice by removing energy in the form of LE from entering the glacier (Fugger et al 2022).

In all three catchments, sublimation removes mass from the ground and canopy intercepted snowpack year-round (figures 3(b) and (c)). Relative to the snowfall amounts, sublimation is largest in Kyzylsu, followed by Langtang and 24 K (15, 13, 6%, respectively) (figure 2(b)). In absolute terms and for the common elevation range, sublimation is also largest in Kyzylsu (91 mm yr⁻¹), but it is second largest in 24 K (71 mm yr⁻¹) and smallest in Langtang (57 mm yr⁻¹). In Langtang, however, where the snowline is generally higher, most sublimation happens above the elevation range of the other basins, amounting to 84 mm yr⁻¹ over the whole elevation range of Langtang (figure S32).

During the warmest months of the year, sublimation rates are low (e.g. Kyzylsu 0.2 mm d⁻¹, Langtang <0.01 mm d⁻¹, Parlung24K 0.03 mm d⁻¹), since snow retreats to the highest elevations, and temperatures are often above the melting point during the periods of high radiation inputs, favouring snow and ice melt rather than sublimation (figure 3(c)). Only in Kyzylsu there is notable snowpack sublimation also in the lower reaches of the catchment, owing to the lower temperatures and heavier snowfalls there during winter (figure 3(c)).

Under warmer- and drier-than-average conditions, such as during the hydrological years 2021/22 (Kyzylsu) and 2016/17 (Langtang) and 2017/18 (24 K), additional ice melt of +35 mm, +33 mm and +10 mm partially compensates for the supply deficits. The supply deficit is here defined as the combination of precipitation deficits and ET losses, which amount to −225 mm, −153 mm and −368 mm compared to the multi-year average (2010/11-2021/22), at Kyzylsu, Langtang and 24 K, respectively (figure S31). Snowmelt under warm conditions has an even greater replenishing effect on the water yield, since old snow, e.g. from glacier accumulation zones, is being depleted: the snowmelt deficit (−32 mm, −22 mm, −170 mm) is considerably smaller than the snowfall deficit (−141 mm, −61 mm, −236 mm) (figure S31), implying a net snowmelt increase of 109 mm, 39 mm and 66 mm. While ice melt alone offsets only around 16%, 22% and 3% of the supply deficits, the ice and snow melt enhancements together offset 64%, 47%, and 21% of the precipitation deficits. As a result of compensations due to enhanced glacier and snow melt, the driest and warmest hydrological years (such as 2021/22 at Kyzylsu, 2016/17 at Langtang and 2017/18 at 24 K) do not correspond to the greatest runoff deficits (figure 4(a)).

Zoom In Zoom Out Reset image size

Evapotranspiration anomalies can reach values similar to those of enhanced ice melt under warmer conditions (e.g. Kyzylsu: 2015, 2021; Langtang: 2021; 24K: 2017; figure 4(b)), implying that evapotranspiration counterbalances the additional ice melt contribution to streamflow. Note however, that the level of this balancing is strongly related to the proportion of glacier cover versus the remaining land cover in a catchment.

Sublimation is highest in both relative and absolute terms and across the common elevation range in Kyzylsu, supporting our first hypothesis (I) that sublimation is most important in semi-arid, high elevation catchments (figures 2 and 3). Clear-sky conditions, with low winter temperatures and high snowfall rates create favourable conditions for sublimation. In contrast, cloudy, warm conditions during the accumulation period drastically reduce sublimation but favour melt (Mandal et al 2022, Fugger et al 2022). Snow sublimation has long been recognized as an important mass flux in snowy alpine environments, including in the Alps (Strasser et al 2008), glacierized catchments in the Canadian Rockies (Pradhanga et al 2022, Aubry-Wake et al 2022), the Himalaya (Buri et al 2023), and at the glacier scale in the Andes (Ayala et al 2017). The share of snowfall removed by sublimation ranges between 5% and 90% in previous literature, depending on the study site's environment, elevation range and topography. Here, we estimate this share for the three catchments as 15% (Kyzylsu), 13% (Langtang) and 6% (24 K). Modelling scenarios not including sublimation would overestimate the annual snowmelt runoff by this share, affecting especially summer season runoff and peak flows (Stigter et al 2018).

Our second hypothesis (II), that ET is higher in monsoonal catchments, can also be confirmed, since Kyzylsu shows less overall ET than Langtang and 24 K. This is mostly a result of the energy regimes and snow cover dynamics, rather than a result of soil water limitations (figure 3). However, intense rainfalls over the greater part of the year provide additional water to evaporate from leaf-, rock- and other interception-storages, boosting ET at 24 K.

In none of the catchments and years did ice melt alone overcome the supply deficit during warm-dry years, although this may happen at the seasonal or event timescale (Van Tiel et al 2021). While ice melt enhancements reach similar magnitudes as the evaporative loss enhancements, the third hypothesis, that ice melt can compensate for the combination of both precipitation deficits and enhanced evaporative losses during warm-dry years, is not supported by our results. Ice melt combined with melt from old snow can however offset as much as 64% of the supply deficit in Kyzylsu. These results reveal important implications for the system behaviour under climate change and should contribute to the discussion around future water yield under a vanishing cryosphere (e.g. Immerzeel et al 2013, Ragettli et al 2016, Huss and Hock 2018, McCarthy et al 2022).

The interannual anomalies of sublimation and evapotranspiration can reach similar magnitudes as the ice melt anomalies (figure 4(b)), depending on the catchment configuration. The interannual variability in evaporative fluxes is generally smaller than the interannual variability of rain- and snowfall inputs at all sites, in agreement with previous studies (Roberts 1983, Fatichi and Ivanov 2014). Under a warming climate, therefore, shifts in precipitation volume and phase (Jouberton et al 2022, Shaw et al 2022), as well as the decline in glacier cover will probably be more important to streamflow variability (e.g. Ragettli et al 2016, Huss and Hock 2018) than changes in the evaporative fluxes. The streamflow at the semi-arid site is likely the most sensitive to climate change, due to the important role of the cryosphere in runoff generation for two reasons: (i) glacier melt sustains summer runoff and glaciers are shrinking, even in the domains affected by the Pamir–Karakoram anomaly (Hugonnet et al 2021). Where 'imbalance' ice melt at present compensates for runoff deficits under warm-dry conditions, this effect will initially increase, but eventually vanish, once ice masses have shrunken further (e.g. McCarthy et al 2022); (ii) Our analysis shows that the melting of old snow can buffer streamflow deficits during warm-dry years to an even greater extent than the melting glacier ice, with the highest degree of compensation at Kyzylsu. These old snow reserves would however vanish quickly over a few consecutive warm-dry years. Additionally, seasonal snowpacks are thinning and global warming is bringing the bulk snowmelt events earlier in the season (Unger-Shayesteh et al 2013). Modelling experiments based on future scenarios are needed, in order to fully understand the High Mountain Asian headwaters' sensitivity to climate change.

Uncertainties in model outputs were not explicitly quantified due to the high computational demand of the model (table S4) limiting our ability to perform a large number of simulations. Instead, the model results were compared against a large range of ground and remotely sensed data, confirming its ability to replicate the necessary physical processes. While sublimation from the ground- and canopy-intercepted snowpack is simulated, blowing snow sublimation is not, and so sublimation on wind-exposed mountain ridges is potentially underestimated in this study (Strasser et al 2008). Similarly, due to the absence of a blowing snow transport scheme, the heterogeneity of snow depths between ridges and valleys is underrepresented (Aubry-Wake et al 2022). The model does not have a separate firn layer representation. This might affect the mass balance of the accumulation zones, since the thermal and surface properties of firn differ from those of snow and ice (Pradhanga and Pomeroy 2022). Supraglacial ice cliffs and ponds were not modelled, but these features are expected to enhance melt-rates of debris-covered glaciers (Miles et al 2018, 2022, Buri et al 2021, Kneib et al 2022). Finally, measurements of the surface latent heat fluxes are needed to better evaluate these fluxes in high elevation catchments and should be a goal of future fieldwork.