Glacier-influenced hydrological regimes in the Afghanistan Hindu Kush Himalaya under current and future climate

Glaciers cover around 30,000 km² or about 17% of the Hindu Kush Himalaya (HKH) region (Sarfaraz et al. 2004) and are one of the most sensitive systems to climate change (Immerzeel et al. 2020). They are also important water resources for drinking water supply, irrigation, hydropower generation, and downstream ecosystems. Notably, they can sustain streamflow during meteorological drought conditions (Marshall 2014; Li et al. 2015; Clason et al. 2023). Such increases are more severe in more arid regions of the HKH where glaciers may compensate for precipitation deficits (Kaser et al. 2010; Tiel et al. 2021).

Changes in glacier volume or area directly affect meltwater runoff in terms of both annual totals and seasonal and inter-annual variability. In the short to medium term, higher air temperatures generally enhance summer streamflow until glaciers become small enough that this ‘glacial subsidy’ ends and summer runoff reduces (Rees & Collins 2006; Scherler et al. 2011). The transition from glacial subsidy to glacial penury is known as ‘peak water’ (Huss & Hock 2018). Such changes will influence water availability in downstream rivers and affect longer-term groundwater recharge. Therefore, determining the impacts of climate change on the cryosphere is essential to understanding future water supply changes (Quincey et al. 2018; Kraaijenbrink et al. 2021). This requires catchment-scale glacio-hydrological models that take into account the uncertainties relating both to models and the climate predictions used to force them (Azam et al. 2021).

In this paper, we focus on Afghanistan in the northwestern HKH, which has a complex climate and has rarely been studied in terms of glacio-hydrological regimes. Winter snowfall is predicted to decline in this region, and summer temperatures will continue to rise in the immediate future (Savage et al. 2009; NEPA & UNEP 2016; Aich et al. 2017). Such changes will impact both the annual mean and the intra-annual distribution of streamflow. However, they will also impact glacier mass balance, the volume of water stored as ice, and hence the magnitude of glacial meltwater supply.

Determining the impact of these changes on runoff in Afghanistan has three broad challenges. First, climate change has spatially variable impacts on glaciers within Afghanistan and surrounding countries. The more strongly westerly-influenced Karakoram glaciers (in the western HKH) are typically stable or advancing slightly (Scherler et al. 2011; Bhambri et al. 2013; Minora et al. 2016; Azam et al. 2018; Farinotti et al. 2020). Those in Pamir-Alay have strongly negative mass balances and are retreating (e.g., Gardner et al. 2013; Kääb et al. 2015; Brun et al. 2017; Shean et al. 2020; Denzinger et al. 2021). This regionally differentiated change in glacier response to climate change will translate into future streamflow regimes. This regional variation remains a primary source of uncertainty in assessing the regional hydrological impacts of climate change (Lutz et al. 2014).

Third, available models applied over multi-decadal timescales rarely consider debris cover effects (Fyffe et al. 2019). Debris-covered glacier response to climate change is variable (Mihalcea et al. 2006; Benn et al. 2012; Reid et al. 2012; Banerjee 2017; Rounce et al. 2018). Debris cover is generally assumed to reduce melt rates as compared to bare ice, except where debris is thin or patchy, where it may locally increase ice melt (Östrem 1959; Mihalcea et al. 2006; Lejeune et al. 2013; Fyffe et al. 2014; Minora et al. 2015; Fyffe et al. 2019). The effect on downstream streamflow can vary significantly even in the same geomorphological, geological, and climatological settings (Miles et al. 2020; Shafeeque et al. 2022). Very few glacio-hydrological modeling studies account for the effect of debris (Fyffe et al. 2019;  2020); and these have tended to focus on smaller catchments with areas <100 km² (Reid et al. 2012; Ayala et al. 2016; Carenzo et al. 2016).

Given the likely spatial variability of catchment hydrological response to climate change in the HKH, the relatively underdeveloped representation of glaciers, and the need to include a debris cover representation, the aims of this study are (1) to implement a suitable glacio-hydrological precipitation-streamflow model suitable for data-poor glaciated mountain settings that includes a debris cover representation and (2) to use this model to estimate how glacier-influenced streamflow in Afghanistan could evolve under future climate change, taking into account glacier melt. This work was undertaken for three catchments in Afghanistan (Figure 1) chosen to represent different climatic regions because they had the most reliable hydro-meteorological and glaciological data. In the paper, we (a) use in situ discharge and climate data to set up a streamflow model, (b) couple it to bias-corrected historical and future climate forcing, (c) incorporate remotely-sensed-based glacier observations to account for current debris cover and debris cover trends, and (d) use a simple approach to represent glacier retreat effects.

Afghanistan had an estimated glacier extent of 2,690.7 ± 108.2 km² in 2020, with 75 ± 0.7% bare ice and 25 ± 3.0% debris-covered ice (Shokory & Lane 2024). It has an arid and semi-arid continental climate characterized by desert, steppe, and mountain temperature and precipitation regimes (Humlum et al. 1959; Shroder 2014). Mediterranean cyclonic systems and westerlies mainly cause winter precipitation (Glantz 2005; Sorrel et al. 2007). In the summer, the South Asian monsoon associated with the intertropical convergence zone influences the country (Shroder 2014; Shokory et al. 2017), with occasional snowfall during summer in the highest mountain peaks in the northeastern region.

Afghanistan has five major river basins (Shokory et al. 2023): (1) Panj-Amu, (2) North, (3) Harrirud Murghab, (4) Helmand, and (5) Kabul (Figure 1) with very different contributions to the total volume of streamflow leaving the country (Bromand 2017). The Panj-Amu and Kabul River Basins produce the highest water volumes (38 and 35% of the total national water, respectively) despite having smaller percentage basin areas (14 and 11% of the country). This is due to glaciers, and these basins account for 96 and 4% of the total glacier area of Afghanistan, respectively (Shokory & Lane 2024). The other three basins have lower streamflow volumes than expected given their catchment area, with 17% by volume but 52% by area for the Helmand River Basin, 5% by volume but 12% by area for the Harrirud Murghab River Basin, and 4.8% by volume for the Northern River Basin but 11% by area (Mahmoodi 2008; Bromand 2017; Shokory et al. 2023) (Figure 1).

The three catchments used for this study are located (1) in the Panj-Amu River Basin, in the central region (the Bamyan catchment); (2) in the northern region (the Taqchakhana catchment); and (3) in the eastern region (the Sust catchment) (Figure 1). The catchments were selected for both practical and scientific reasons. Practically, these three basins have relatively reliable streamflow data, necessary for model calibration and evaluation. Scientifically, they represent three different degrees of glacier cover (Table 1) reflecting climatic differences. The Sust catchment has the highest glacier extent (15.6% by area), the lowest daily mean temperature (−0.3 °C), and the lowest mean annual precipitation (320 mm year⁻¹) (Table 1). The Taqchakhana catchment has the highest daily mean temperature (+9.2 °C) and mean annual precipitation (1,397 mm year⁻¹) but a lower glacier extent (2.8% by area). The Bamyan catchment has the lowest glacier extent (0.7% by area), an intermediate mean daily temperature (+4.2 °C), and mean annual precipitation (414 mm year⁻¹) (Table 1).

This section presents the glacio-hydrological and debris cover evolution models, the data used for model calibration and evaluation, and the data and its bias correction used for climate change scenario simulation.

The model structure is a modified version of the Glacier and SnowMelt – SOil CONTribution (GSM-SOCONT) glacio-hydrological model, which is a semi-lumped, reservoir-based precipitation-streamflow model (Schaefli et al. 2005). It transforms precipitation into streamflow at the scale of hydrological units here defined as elevation bands of 50 m elevation, based on the 30 m resolution digital elevation model (DEM) (Shimada et al. 2009, 2024). These hydrological units are further subdivided into landscape units. Here, the selected units are ice and non-ice-covered catchment parts, the ice unit being further subdivided into debris-covered and non-debris-covered parts. The three units are defined as dynamic fractional land cover units. The model is semi-lumped, which means that the model parameter values are not distributed in space, but all elevation bands share the same parameter values for a given landscape unit, obtained via calibration on observed streamflow. The model uses a fuzzy threshold for snow accumulation (<1 °C) and a temperature-index-based snow and ice melt model. The transformation of melt and rainfall into streamflow is achieved here with three linear reservoirs (Figure 2), i.e. with reservoirs where the outflow is proportional to the degree of filling. Of these reservoirs, two are in parallel, the third one is in series, which allows an emulation of a quick runoff component and two slower runoff components. The model runs at a daily time step and has input precipitation, air temperature, and potential evapotranspiration (PET), which is obtained with the Hargreaves & Samani (1985) method. Air temperature and precipitation are interpolated to the mean elevation of each band based on the given temperature lapse rate and precipitation gradient (see Section 3.4.3). Precipitation-streamflow transformation is then computed at the scale of each elevation band for each of the land cover units (bare ice, debris-covered ice, and non-ice) and summed up to build streamflow at the catchment outlet. For a detailed description of the model, the reader is referred to Horton (2023). Here we focus on summarizing the essential components.

Compared to the underlying original GSM-SOCONT model, the Hydrobricks version used here can handle multiple land cover types and a dynamic evolution of the land cover fractions over time; the model also accounts for the temporal evolution of debris-covered ice areas. Furthermore, we use here three reservoirs instead of two to enhance the representation of groundwater-fed baseflow. As in almost all similar high-elevation modeling studies, the model does not account for precipitation interception by vegetation, which we assume does not play a major role in this climate due to minimal vegetation cover.

We considered three primary land cover types to characterize streamflow contributions: non-ice land cover, bare ice, and debris-covered ice. Non-ice land cover includes all land cover types except glaciers, which are represented by bare ice and debris-covered ice, with glacier outlines produced by Shokory & Lane (2024). These distinctions are needed for accurate hydrological modeling, as some models do not differentiate between bare ice and debris-covered ice, nor do they allow for temporal changes in ice cover (Nepal et al. 2014; Azam et al. 2021). A failure to represent debris cover and to allow glacier contributions to decline may substantially overestimate streamflow (Tiel et al. 2020).

To address this challenge, we adopted two strategies. First, we accounted for different degree-day factors (Hock 2003) for bare ice and debris-covered ice (e.g. Kayastha et al. 2000; Hagg et al. 2008; Immerzeel et al. 2015). The limits of the degree-day approaches for modeling glacier melt are well known (Nicholson & Benn 2006); but, given the poor data availability and requirements of more complex methods, the degree-day method was deemed to be parsimonious and also allowed us to include uncertainty analysis in the modeling.

Second, we developed simple empirical glacier evolution scenarios based on remotely sensed measurements of bare ice and debris-covered ice evolution for our catchments from 2000 to 2020 (Shokory & Lane 2024). Future glacier evolution was calculated using a simple linear change to the end of the 21st century. While physically based models depend heavily on detailed observational data from individual glaciers to initialize, calibrate, and evaluate model accuracy. Consequently, these methods are generally applied to single glaciers rather than large samples (Zekollari et al. 2022). Additionally, physical models require detailed bed topography, reconstructed through field techniques such as borehole drilling or radio-echo sounding, to calibrate and evaluate the results (Zekollari et al. 2022). Given the labor-intensive nature of these methods, they are limited to only a few point measurements and are rarely used to derive ice thickness over a large number of glaciers (Li et al. 2017; Zekollari et al. 2022). Given security issues, undertaking fieldwork in Afghanistan is almost impossible. Due to these data and resource constraints, simplified empirical methods remain widely used for modeling large numbers of glaciers globally (Zekollari et al. 2022).

To assess the effects of debris-covered ice on streamflow prediction, the model is set up and calibrated twice: (1) with the distinct representation of bare ice and debris-covered ice and (2) with the representation of bare ice and debris-covered ice summed as a single land cover type.

The glacio-hydrological model presented here has 12 parameters to calibrate (Table 2) against observed streamflow or other observed data (such as snow data or glacier mass balance data), which are, however, not available for our study. Some parameters were kept constant, based on other studies, or were estimated based on observed data (Schaefli et al. 2005; Horton 2023). The monthly temperature lapse rates and elevation threshold for precipitation were estimated from observed data (Table S1, Supplementary Material). Given the well known challenges of equifinality (Beven & Freer 2001), we did not select for an algorithm-based parameter optimization but generated 30,000 parameter sets in a priori possible parameter ranges taken from previous work (Schaefli et al. 2005). The performance of the model for each of the sets was evaluated with the Kling Gupta efficiency (KGE) value closer to 1, and the best 100 parameter sets were retained for further analysis. We use KGE rather than Nash–Sutcliffe efficiency (NSE) (Kling et al. 2012) because KGE provides a balanced evaluation of the model performance between low flows and high flows by incorporating correlation, variability, and bias, making it particularly suitable for our study that focuses on streamflow components estimation rather than flood peaks (see Section 5.5). In addition to KGE, we provide additional metrics like Percent Bias (PBIAS), Root Mean Square Error (RMSE), and Mean Absolute Error (MAE) to give a full picture of the model performance (Moriasi et al. 2007) (Table S2, Supplementary Material).

The calibration and evaluation periods varied among catchments depending on data availability (Table 1). For the Sust catchment, calibration and evaluation were conducted for 2013–2015 and 2016–2017, respectively. For Taqchakhana, these periods were 2015–2017 and 2018–2019, and for Bamyan, 2012–2015 and 2016–2019. We used a 6-month spin-up period, with the datasets starting from the first of January that consider the model effective run in the middle of the summer (which reduces the impact of not initializing the snowpack); such spin-up periods are commonly used to remove state-variable initialization problems.

We used the Shuffled Complex Evolution Algorithm (SCE-UA) for automatic parameter optimization (Duan et al. 1992), which is now considered one of the most robust and efficient algorithms for parameter calibration (Seong et al. 2015). The model runs on daily time steps and is iterated 30,000 times for each catchment to determine permissible parameter values. It uses the modified KGE criterion as the objective function (Kling et al. 2012). The KGE is based on the equal weighting of three subcomponents: linear correlation between observed and simulated discharge, the bias ratio, and the ratio of standard deviations of the simulated and observed streamflow (Kling et al. 2012). The SCE-UA was applied using the Statistical Parameter Optimization Tool for Python (Houska et al. 2015). To assess the performance of the model beyond a calibration period, we divided the available data (Table 1) into two periods, one for model calibration and one for evaluation, i.e. for the assessment of the quality of model predictions beyond the calibration period. The model with two optimized parameter sets, one with the debris cover melt treated separately and one with it combined, was then run from 2000 to 2100 using the bias-corrected climate scenarios for two Representative Concentration Pathways (RCPs, Table 3) which we define in the following.

The model requires four main input datasets for each catchment: (i) hypsometric distributions of elevation for each land cover (bare ice, debris-covered ice, and non-ice), (ii) daily precipitation, (iii) daily mean temperature, and (iv) PET estimated for at least one meteorological station. The availability of observed streamflow and climate data, such as precipitation and temperature, is limited to the periods outlined in Table 1. The glacier change percentages reflect variations between two distinct points in time (2000 and 2020) rather than continuous annual data (Table 1).

The catchment hypsometry was computed based on a DEM with a resolution of 30 m using Advanced Land Observing Satellite (Shimada et al. 2009) data. Glacier shape files and retreat rates between 2000 and 2020 were provided by Shokory & Lane (2024). Meteorological data (temperature and precipitation) and daily streamflow were obtained from two Afghanistan governmental entities, the Ministry of Energy and Water and the Ministry of Agriculture Irrigation and Livestock (Figure 1, Table S1, Supplementary Material). A total of 41 meteorological station data time series were obtained (at elevations <∼3,000 m.a.s.l). Each catchment contained one of these stations at the catchment outlet where streamflow was also measured. These were used in the modeling. These and the remaining meteorological stations were also used to calculate precipitation gradients and temperature lapse rates for the catchments during data preparation. The PET data were calculated based on the Hargreaves–Samani method using maximum and minimum daily temperatures (Hargreaves & Samani 1985) recorded at the meteorological stations.

Subsequently, to analyze the impacts of climate change on future water availability, we forced the model with the latest ensemble of climate models. We collected regional climate model (RCM) outputs driven by eight general circulation models (GCMs) from the Coordinated Regional Climate Downscaling Experiment (CORDEX-CMIP5) (Table 3) and applied bias. We chose CORDEX-CMIP5 over CMIP6 data due to its higher spatial resolution (approximately 22 km) and regional focus, necessary to capture localized climate details in our study area. While CMIP6 provides valuable global-scale data, it typically operates at a coarser resolution (100–250 km) (Sediqi et al. 2022; Rahimi et al. 2024), which limits its effectiveness for regional analysis (Gutowski et al. 2016; Giorgi 2019). Additionally, CORDEX-CMIP6 data for the South Asia domain (WAS-22) is still under development and unavailable for our study region (https://wcrp-cordex.github.io/simulation-status/CORDEX_CMIP6_status.html#WAS-22).

Thus, CORDEX-CMIP5 currently offers the most suitable dataset for our research objectives. Two RCPs were used; 2.6 (the concentration of carbon that delivers global warming at an average of 2.6 W/m² across the planet 2.6); and 8.5 (8.5 W/m²). We used RCMs from the South Asia domain (WAS-22) at a daily time step with a horizontal resolution of 0.22° (about 22 km) provided by the Earth System Grid Federation (ESGF https://esgf-node.llnl.gov/search/cmip6/). We used three climate variables from CORDEX; near-surface temperature (tas), precipitation (pr), and PET (evspsbl) at the daily time scale. Data were obtained for the location of the outlet of each catchment (Table 2). Historical periods (1950/1970–2005) were not considered in the analysis, and we had only data for recent years to perform bias correction. Two periods were considered for the climate change assessments: mid-future (2040–2060) and late future (2080–2100) to compare with a baseline period of 2000–2020.

In this study, observed station data is used to set up the mode; however, we avoid using reanalysis datasets. While these datasets provide consistent meteorological data globally, they are not free from uncertainties, particularly in regions with a complex topography and limited observational data, such as Afghanistan (Rahil et al. 2024). And we had limited data to perform bias adjustment if reanalysis data was used (Gudmundsson et al. 2012; Beck et al. 2019); therefore, we decided to rely on observed data for this study.

The quality of observed precipitation and streamflow data of the catchment stations was assessed using double-mass curves at the daily time scale to evaluate data consistency (Searcy & Hardison 1960).

For precipitation data, this was achieved by comparing the individual catchment stations with the average of many records from other stations of the same basin (Figure 1) to identify unexpected slope breaks in the double-mass curve. Then, the best reference station with no inconsistency was selected to adjust the catchment station datasets. For streamflow data, we also used a double-mass curve analysis for individual streamflow datasets of the catchment stations with the average of many station records (Figure S1, Supplementary Materials). This technique was applied only to the Bamyan and Taqchakhana stations, where neighboring stations were available, but not to the Sust station (Figure 1).

RCMs have biases compared to ground stations, partly inherited from the GCMs used to drive them (Rummukainen 2010; Hall 2014; Maraun 2016). Such biases were also observed for the three catchments used in this study for precipitation and temperature data (Figures S2, S4, Supplementary Material). To remove these biases and to adjust the distribution of RCMs to resemble the observed data, we applied the quantile mapping bias correction method, frequently used in water resources studies (Gudmundsson et al. 2012). Gudmundsson et al. (2012) presented three transformation approaches: (1) distribution-derived, (2) parametric, and (3) nonparametric, in which they obtained the highest skill in systematically reducing biases with nonparametric transformation, which we applied in our study. RCM data were first downloaded for the grid cell closest to the observed station. The station was always located at the lowest elevation point of the catchment (the outlet) and served as the reference point. We then performed bias correction of the RCM data using observed station data at this location. The first (lowest) elevation band was subsequently set as the reference for the RCM grid cell. Using this reference, we applied precipitation gradients and temperature lapse rates to interpolate climate-forcing data across all elevation bands within the catchment. The temperature lapse rate from station data is used for future scenarios; however, the precipitation gradient has been optimized. This method ensured that the spatial variability in climate forcing due to elevation was represented for this data-sparse setting.

Precipitation gradients and temperature lapse rates have a significant impact on the water balance in glacio-hydrological models, especially in mountain regions, and may vary based on elevation and season (Immerzeel et al. 2014). Afghanistan rivers are generally fed by snow and ice at elevations > 3,000 m a.s.l. At these elevations, data are scarce (Shokory et al. 2023). This makes the estimation of lapse rates challenging. In addition, complex topography, the dominance of solid precipitation, and strong winds at high elevations likely result in unreliable precipitation measurements (Freudiger et al. 2017). Moreover, there is a known considerable regional variation in the precipitation gradient in the mountainous region of Afghanistan (Azizi & Asaoka 2020).

To obtain reasonable catchment-scale precipitation estimates, elevation-dependent precipitation gradients have to be estimated based on available station data (Immerzeel et al. 2014). In our study, the precipitation gradient was calculated based on the observation from all 41 stations. Observed data showed that at 2,000 m.a.s.l. there was a breakpoint in the relationship between precipitation and elevation (Figure S5, Supplementary Material); however, above this elevation the precipitation points were fuzzy. To define this relationship above 2,000 m.a.s.l., we followed Immerzeel et al. (2014), who reported a similar breakpoint in the HKH. They noted that this breakpoint coincides with a transition from plains to mountainous terrain, where elevation increases but area coverage decreases, found to be close to the modal elevation of the hypsometric curves. We used this as the breakpoint, but further work is needed to assess the suitability of our approach and to understand how sensitive model predictions are to it.

In each catchment, there was only one precipitation recording station at the lowest elevation. Therefore, it was not possible to determine (lower elevation precipitation gradients, δ_precip,1) and (higher elevation precipitation gradients, δ_precip,2), and these had to be treated as parameters that required calibration. The calibration process profited from the observation that if precipitation was too high at high elevations, then the model would progressively over-accumulate snow at high elevations over a number of years, in what are known as ‘snow towers’ (Freudiger et al. 2017). To avoid building ‘snow towers’, initial model simulations focused only on snow accumulation to identify maximum possible precipitation gradients for both δ_precip,1 and δ_precip,2 ; as an outcome, in the calibration process; 0.1 mm/100 m in the Sust; and 0.2 mm/100 m for the other two catchments were identified. The minimum possible values were set to 0 in all cases.

The need for monthly variation in temperature lapse rates is well known (Aizen et al. 2000; Immerzeel et al. 2014; He et al. 2015). Thus, for the temperature lapse rate, we assessed the observed mean monthly temperature data for 23 out of the 41 stations with available temperature measurements (Figure 1) (Table S1, Supplementary Materials) for each month of the year to calculate the lapse rates.

The double-mass curves suggested a slight inconsistency in observed precipitation and streamflow data for Bamyan station (Figure S1b, d, Supplementary Materials). These were corrected (Figure S1a, c, Supplementary Materials) using the Keshem and Balaye Kelagai stations (Figure S1e, f, Supplementary Materials). In addition, we observed a significant inconsistency in the streamflow of the Taqchakhana catchment (Figure S1g, Supplementary Materials). This inconsistency was adjusted by comparison with the Khwajaghar station (Figure S1h, i, Supplementary Materials).

Substantial overestimations of precipitation for all three catchments, mainly for March, April, and May (MAM) and June, July, and August (JJA), were suggested as compared with RCM output (Figure S2, Supplementary Materials). The projected temperature data from the original RCMs underestimated mean temperature compared with observed data for the Sust and Taqchakhana catchments (Figure S4, Supplementary Materials) and overestimated them for the Bamyan catchment. Bias correction considerably improved RCM precipitation and temperature estimates (Figures S3, S4, Supplementary Materials). However, we found that precipitation data from RegCM4-7_v0 continuously overestimated summer precipitation after the bias correction. Therefore, we removed it from further consideration in the modeling and used only five GCM-derived projections. The temperature data of three catchments showed a more acceptable fit with observations after the bias correction (Figure S4, Supplementary Materials).

The breakpoint in the precipitation–elevation relationship for the Taqchakhana catchment was taken as 2,000 m.a.s.l. The threshold was higher for the Sust (2,800 m a.s.l.) and much higher for the Bamyan (4,800 m a.s.l.) (Figure S6, Supplementary Materials).

The optimized lower elevation precipitation gradients (δ_precip,1) were smaller than the higher elevation gradients (δ_precip,2) for the Taqchakhana. However, for the Sust, they were more similar to each other, and for the Bamyan catchment in the center of Afghanistan, optimal precipitation gradients were the same for both low and high elevations (Table 4). A comparison of the calibrated precipitation gradient with observational data was conducted to evaluate model performance. Analysis of precipitation data from 21 stations below 2,000 m a.s.l. (Figure S5, Supplementary Materials) showed a precipitation gradient of 0.04 mm/100 m with an R² value of 0.63. The model calibration for Taqchakhana resulted in a higher gradient of 0.07 mm/100 m (Table 4). It should be noted that the 21 stations cover a broader area than Taqchakhana, where local conditions may be different.

The monthly mean temperature data for 23 observed stations showed important differences in the lapse rates for each month of the year (Figure S7, Supplementary Materials). A lower lapse rate was observed in winter and early spring months: 0.28 °C/100 m in December, rising to 0.55 °C/100 m in June.

Optimal model parameter sets (Table 4) were obtained after 30,000 iterations for each catchment (Sust, Taqchakhana, and Bamyan) for the model version that included bare ice and debris-covered ice separately. After calibration, the KGE values for the daily streamflow series were more than 0.8 for all catchments (Table 4). For the Sust, the KGE improved for evaluation, but it declined for the Taqchakhana and the Bamyan catchments, although the values remained high (Table 4).

Table 4 compares separate bare ice and debris-covered ice representation (Table 4(a)) with the single glacier land cover representation (Table 4(b)). When the separate debris-covered ice representation is removed, the melt factor for ice declines strongly in the most glaciated catchment (the Sust) and hardly changes in the least glaciated catchment (the Bamyan). This is expected as a failure to represent debris-covered ice would be compensated by a lower catchment-wide ice melt factor. Most other parameterizations do not change much, the notable exception being some of the Taqchakhana parameters, likely reflecting the calibration process having difficulties in finding a single best parameter set.

The distribution of water production by elevation band (weighted by the land cover fraction and the elevation band area, Figure 5) suggests that the Sust catchment produces water mainly through glacier melt in the 4,000–6,000 m.a.s.l elevation range, with a peak at 5,000 m.a.s.l. The Taqchakhana catchment streamflow originates predominantly between 4,000 and 5,000 m.a.s.l., along with a substantial contribution from rain and snow runoff between 3,000 and 4,000 m.a.s.l elevations, mainly in the spring months (Figure 5(b) and 5(e)). The Bamyan catchment has a substantial contribution between 4,000 and 5,000 m.a.s.l and substantial baseflow at lower elevations (Figure 5).

This section presents future projected climate (precipitation and temperature) and streamflow for two RCP (low emission 2.6 and high emission 8.5) averaged over two timespans; 2040–2060 labeled as the 2060s; and 2080–2100 as the 2100s. The baseline is considered as 2000–2020, the 2020s. The simulations for each time line compared to the baseline are tested for significant differences with the Mann–Whitney test (Tallarida et al. 1987). The results are illustrated for two timescales: (1) at the annual scale with uncertainties calculated, including the standard deviation of the five ensembles (i.e. Table 2, except RegCM4-7_v0) and (2) at a seasonal scale that shows all data.

Table 5 also shows changes in precipitation. These are less clearly statistically significant and only clear to the 2100s under RCP 8.5. The exception is for the Bamyan catchment, where, under both RCPs (2.6 and 8.5), significant decreases in mean annual precipitation are predicted by the 2060s (−37.0, −62.4 mm year⁻¹ by 2060s) (Table 5, f), primarily in spring and autumn (Table 5). The Sust is predicted to have statistically significant precipitation increases, but only under RCP 8.5 (Figure 6(d), Table 5). The Taqchakhana is expected to have significant precipitation decreases in spring under RCP 8.5, at the annual scale at both the 2060 and 2100s horizons, and in summer at the horizon 2100s. The magnitude of changes for the Taqchakhana catchment is large but not always significant; this results from the fact that this catchment has the highest inter-annual variability of the baseline and future scenarios.

The prediction for the Bamyan catchment suggests decreases in all mean annual streamflow components under both RCP 2.6 and RCP 8.5 and for both simulated periods (Figure 8 and Table 6), except for projections of rain/snow runoff in spring under RCCP 2.6 to the 2100s. The Taqchakhana catchment is predicted to have increasing glacier runoff under RCP 8.5 to both the 2060s and 2100s; but under RCP 2.6, the 2100s are marked by decreases in glacier runoff. There is a marked change in seasonality associated with increased rain/snow runoff in winter and with declining rain/snow runoff in summer (RCP 2.6 to the 2100s excepted). Under RCP 2.6, there is some evidence of higher winter and spring baseflow and lower summer and autumn baseflow. Table 5 shows that the dominant effects of these changes are winter increases and the summer/autumn decreases in snow/rain runoff, which are then seen in major changes in total streamflow in the Taqchakhana in summer, even if the annual streamflow changes are only significant for RCP 8.5 to the 2100s. The Taqchakhana does not seem to have sufficient glacier cover at ‘productive’ elevations (where temperature, even if increasing, is sufficient to cause melt) to compensate fully for declining winter snow accumulation and subsequent summer runoff.

The most glacier-covered catchment, Sust, differs in many senses. Annual total streamflow is predicted to rise under both RCPs, but notably RCP 8.5, and to the 2060 and 2100s. Under RCP 2.6, the changes are much smaller than under RCP 8.5. The generalized increases in total streamflow can be found across all seasons, except under RCP 2.6 and for winter and summer in the 2100s. These increases in flows can also be seen to differing degrees in all streamflow components, and this reflects the implications of wetter conditions, notably for RCP 8.5 (Table 5), which will impact both baseflow and rain/snowmelt runoff, and warmer conditions (Table 5), which should increase glacier melt in this catchment with the highest ice cover.

We conducted a seasonal comparison of observed, model-simulated, and RCM-simulated streamflow for the periods during which observed data were available in each catchment. Although the observational data are not extensive enough for a comprehensive climatological assessment of historically simulated streamflow from RCMs, this analysis, detailed in the Supplementary Material (Figure S8), shows that the historical RCM simulations effectively replicate the climatological patterns in the Sust and Taqchakhana catchments. However, a slight overestimation of autumn streamflow was observed in the Bamyan catchment. This comparison provides a basis for evaluating historical model performance within the constraints of the available data.

In contrast, for the Taqchakhana catchment, we obtained considerable differences in glacier runoff with and without debris cover representations (Figure 10(c), and (d)). More significantly, the differences were obtained in the summer and autumn seasons, showing about 100% more reduction in glacier runoff without debris cover representations than with debris representation (Table 7). Glaciers in this catchment are 40% debris-covered (calculated from Table 1) and have a higher ice melt factor than the Sust (Table 4(a) and 4(b)).

Table 7

Model results without debris cover representations, changes in streamflow runoff components (mm) under two RCPs (2.6 and 8.5) in two timespans (2060 and 2100s) as compared with the 2000–2020 baseline

Note. Statistically significant changes (p < 0.05) are marked for a two-tailed student's t-test.

View Large

In the last 20 years (Figure 10(a)–10(d)), the glacier runoff simulated with a separate debris treatment exceeds that without a debris treatment scenario. The melt factor for the without debris treatment scenario is slightly lower (by 0.49 and 0.17 mm/day/°C for Sust and Taqchakhana, respectively, Table 4(a) and 4(b)) and is within the range of model uncertainties, but the focus here is not on the melt factor itself. This analysis illustrates that without debris treatment, glacier contribution to streamflow might be lost sooner, with a higher initial contribution in the early years, followed by a rapid decline in the later period.

This paper has revealed substantial differences in glacio-hydrological regimes for three Afghan catchments with different climates and geographical settings. The driest catchment, the Sust in eastern Afghanistan, had the largest contribution from glacier runoff, estimated at 76% under current conditions. The Hindu Kush Karakoram and the eastern Afghanistan region get most winter precipitation through westerly and southwesterly air masses (Lutz et al. 2014), so some climatic differences might be expected compared with the other two catchments.

The streamflow was dominated by rain and snowmelt-driven runoff in the Taqchakhana catchment in the northwestern part of Afghanistan (50% contribution to total annual streamflow). This catchment is at the margin between cold and temperate and dry and hot summer climates, mainly receiving precipitation as snow during winter. Glacier melt runoff is lower (33% contribution to the streamflow) but still important.

The Bamyan catchment in the center-west had a completely different hydrological regime. This catchment presented the highest contribution from baseflow by 61%, is located in an arid, steppe, cold climate, and receives much lower annual precipitation.

Overall, the relative contributions of snowmelt and glacier melt to streamflow are variable between catchments. It is not possible to generalize that either glacier melt or snowmelt (cf. Kraaijenbrink et al. 2021) is more dominant, and future water resource projections need to be sensitive to the geographical location of catchments within the HKH mountain area. Generalizations also need to be sensitive to scale effects: as the spatial extent of a catchment increases (with distance downstream), the proportion of the catchment that is glaciated decreases, and so will the degree of influence of glacier melt to water supply (Quincey et al. 2018).

For Afghanistan, very few studies have applied bias corrections when simulating future hydrological responses (Shokory et al. 2023). With recently available observed data, we provide one of the first studies to do so (data are available as supplementary material). Temperature is projected to increase for all catchments and seasons by 2050 but then to reduce until the end of the 21st century under RCP 2.6, while constantly increasing under RCP 8.5. Precipitation was predicted to decrease in spring and autumn for the Bamyan catchment under RCP 2.6 by the 2060s, but fewer changes were significant for the Sust and Taqchakhana catchments. In comparison, precipitation is predicted to reduce significantly for the spring season under RCP 8.5 for the Bamyan and Taqchakhana catchments by 2060 and 2100s.

Such climate forcing greatly affects the glacio-hydrological regime of the catchments. Warming will reduce winter snow accumulation and increase winter runoff; reduced snow accumulation will then reduce summer runoff. Where the warming affects snow accumulation on glaciers, a function of glacier hypsometry, reduced snow accumulation should lead (1) to earlier exposure of glacier ice during summer melt, earlier albedo reduction, and enhanced melt and (2) albeit only if accumulation zones are at a low enough elevation, a reduction in mass accumulation. These together should enhance glacier loss (Azam et al. 2021), which in turn should impact the hydrological regime where glacier melt provides an important contribution. Our study shows that the Taqchakhana, in particular, is likely to see much stronger seasonality, with wetter winters and drier summers (Table 6). In terms of runoff, the Taqchakhana does not have a sufficient glacier cover (Table 1) for rising temperatures (Table 5) to cause a compensatory increase in runoff from melting ice, and so summer streamflow is predicted to decline significantly (Table 7).

In contrast, the Sust catchment is projected to have an increase in precipitation under RCP 8.5 for all seasons, which will increase high-elevation snow accumulation, leading to less negative or even positive mass balance (Azam et al. 2021). Data on glacier change in Afghanistan confirms this is currently the case for glaciers in this region (Shokory & Lane 2024). This is likely to be countered in the longer term if temperature continues to rise and ablation increases (the Sust is too high for winter temperature rise to cause a strong shift from snowfall to rainfall and so modify accumulation). Our simulations of streamflow under RCP 8.5 do suggest an increase in glacier runoff for the Sust catchment (reflecting higher glacier melt rates) as well as increases in baseflow and in rain/snowmelt runoff due to increasing precipitation.

For the glacier retreat scenarios considered here, it is difficult to see clear evidence of peak water except for the Sust under RCP 2.6, where there is a peak in annual glacier runoff around the 2050s (Figure 8(a)). Otherwise, glacier runoff consistently rises for the Sust and Taqchakhana under RCP 8.5 (Figure 8(d) and 8(e)) and declines for Bamyan under both scenarios (Figure 8(c) and 7(f)) and does not change appreciably for Taqchakhana under RCP 2.6 (Figure 8(b)). There is little research that has identified unequivocally the percentage to which glacier cover in a river catchment must decline for peak water to occur. That is not surprising because it is likely also influenced by other parameters such as basin hypsometry, debris cover, glacier aspect, as well as local climate forcing and its evolution in the future. For instance, if we take the warmest scenario (RCP 8.5) for the Sust catchment, increases in winter precipitation should, on average, sustain winter accumulation. There is a progressive increase in glacier melt (Figure 8(d)), but the ice cover is so high, and its rate of decline is likely reduced via increases in winter accumulation, that it is not yet possible to see peak water. Of course, this result depends on how we model future glacier change, and it suggests the need to look further at what is happening in some of these regions of the HKH with high elevations, high glacier cover, high debris cover, and likely future wetter conditions.

By the 2050s, under RCP 2.6, temperatures in the Sust may have stabilized or even been declining (Figure 6(a)); a similar trend is observed in the Hindu Kush Himalayan region as well (Chaturvedi et al. 2014). This is, of course, a highly optimistic scenario, but it is likely that, given there is still extensive ice cover in the Sust in the 2050s (Figure 3), the peak in glacier runoff is not a result of ice cover falling below a critical level but of stabilizing or even decreasing temperature.

Our study is one of the first glacio-hydrological modeling studies that has considered debris cover effects on streamflow at larger scales in this region despite the need to do this (Huss & Hock 2015; Compagno et al. 2022). More complex treatments require extensive input data relating to local surface mass balance parameters, glacier thickness (ice and debris-covered), and so on, which are not available at the spatial extents of the catchments used in this study, especially given the data-poor setting of Afghanistan.

We found different effects of separating debris cover representation with catchment characteristics. No significant effects were found for the Sust catchment because of the lower percentage of debris cover. More important impacts were noted for the Taqchakhana catchment; glacier runoff contribution without debris representation resulted in higher annual glacier melt than with a debris representation under both RCPs (2.6 and 8.5), but this difference declined over time. This effect was not manifest in Sust because whilst it has a high ice cover, it has generally lower proportions of debris-covered ice. The effect was also not present for Bamyan because its glacier extent is so small that any change in glacier representation has little overall impact. Thus, the conclusion is that debris cover matters for ice melt and streamflow simulations at the catchment scale and has to be moderated by the relative importance of bare ice and debris-covered ice in a catchment and not just by the presence of debris-covered ice. If the extent of a glacier upstream of a selected point in a catchment is extensive and if this is a catchment that can produce significant amounts of debris, which is itself spatially variable (Shokory & Lane 2023), then debris may insure against rapid ice loss due to global warming. In such situations, the correct representation of debris is likely necessary.

As for future streamflow and partitioning of streamflow into different sources, these results are sensitive to the way we parameterized changing bare- and debris-covered ice (Figure 3). For instance, with a greater rate of bare ice loss and debris-covered ice gain, the Sust might see lower rates of glacier runoff in the 2100s. Equally, our representation of bare and debris-covered ice is not temperature-dependent (as it is an empirical extrapolation), so it is not sensitive to the RCP we used. Given differences in future temperature under these scenarios (Table 5), this is a particular weakness. This emphasizes the need for a more physically based representation of debris cover formation (e.g. Nicholson et al. 2021).

In this study, we assumed a constant, linear rate of glacier evolution based on current observations primarily due to the limited observations of physical glacier parameters (Shokory et al. 2023; Shokory & Lane 2023; 2024). This simplifies the representation of glacier contributions to streamflow under future climate conditions, which is a common constraint in regions with limited glacier monitoring data (Huss & Hock 2015). While this approach is consistent with observed glacier dynamics, it may lead to an underestimation under warmer future conditions, where glacier melt rates are likely to increase non-linearly as a response to accelerated melting (Huss & Hock 2015). As a result, our model may overestimate future glacier melt contributions in the short term and underestimate them in the long term as glacier volume reduces.

Additionally, our model parameters were calibrated based on current climate conditions, which introduces potential limitations when extrapolating to significantly different future climates (as in all similar climate change impact studies). While our conceptual, process-based hydrological model effectively captures present-day relationships between climate variables and streamflow, the reliability of calibrated models under scenarios involving substantial shifts in temperature, precipitation patterns, or glacier dynamics remains uncertain (Andréassian et al. 2004). However, there is likely no fundamental shift in dominant hydrological processes since the simulated system is not transitioning into a completely new hydroclimatic regime (which would, e.g. be the case if there was no seasonal snow cover in the future). The parameters of temperature-index-based approaches are known to vary depending on the dominant components of the energy balance (Huss et al. 2009). The calibration of the melt parameters on additional data (e.g. remotely sensed snow data, He et al. 2015) covering long time periods could help to address this limitation by assessing the robustness of the parameter sets for different climatic conditions. This will become possible in the future. In exchange, the hydrological parameters describing the filtering effect of the catchment (how liquid water is transformed into runoff) are likely to also hold in the future unless there is a significant shift in the dominant geomorphological and geological features (e.g. in the stream network) that explain this filtering effect. Furthermore, the contribution of snow and glaciers to streamflow in the modeling setup could not be directly verified due to the lack of observed data specific to these components. This highlights the need for future research to incorporate measured data for validation purposes.

Hydrological model evaluation requires aligning performance metrics with study objectives. While the NSE is widely used, it often prioritizes peak flow accuracy at the expense of baseflow representation due to its emphasis on variance, leading to the underestimation of lower flows (Gupta et al. 2009; Mizukami et al. 2017). In contrast, the KGE provides a broader evaluation by incorporating correlation, variability, and bias, making it more suitable for objectives like streamflow contribution estimation rather than flood peaks.

In this study, we assessed the glacio-hydrological regime of three catchments in Afghanistan that represent different geographical contexts in terms of climate, glacier cover, and scale (catchments ranging from 264 to 4609 km²). We used a parsimonious conceptual precipitation and ice melt-runoff model combined with a new representation of debris-covered ice effects to assess current and future streamflow components. Climate change scenarios are based on bias-corrected RCM outputs for RCP 2.6 and RCP 8.5 of the recent CORDEX simulations. First, we found that the pre-processing of streamflow and climate data is a key step to setting up a hydrological model in comparably data-sparse areas. Our results show that particular attention should be paid to elevation precipitation gradients and temperature lapse rates, which differ spatially and temporally. Second, as shown for other regions, bias correction is of prime importance for reliable hydrological modeling of the regions of Afghanistan impacted by snow and glacier melt. Third, whilst all studied catchments had seasonal streamflow patterns, with runoff peaks in June and July, these result from very different contributions of rain and snow runoff versus glacier runoff. The latter varied significantly between catchments and emphasizes the need to contextualize future water resource changes due to rapid climate change in terms of both glacier cover and local meteorological drivers.

Fourth, the primary climate change signal under both RCP 2.6 and 8.5 was changing temperature, with precipitation changes being less important. As expected, the rate of change under RCP 8.5 is greater and continuous through to the 2100s, but under RCP 2.6, temperatures stabilize or begin to decline. The changes in precipitation were found to be drying in the westerly catchments, but only significantly in spring; and under RCP 8.5, on the other hand, for the easterly catchment, an increase in precipitation is predicted for all seasons to the 2100s. These results emphasize that the impact of climate change on precipitation in Afghanistan is likely to be regionally variable, potentially over small distances.

Fifth, streamflow projections associated with these possible climate changes suggest differences in likely basin response along a gradient from the most westerly catchment considered, where there is a tendency toward reduced annual streamflow, to the northern catchment, where there is no apparent change in annual streamflow but a clear change toward increased seasonality to the eastern catchment where annual runoff is projected to increase. These patterns are reinforced under RCP 8.5 as compared with RCP 2.6. Projected total precipitation decreases are evident in both baseflow decreases and rain and snow runoff over Taqchakhana and Bamyan catchments (Tables 5 and 6). Glacier runoff increases substantially in the most glaciated eastern catchment under both scenarios and in the northern catchment under RCP 8.5. The extent to which runoff is increased by glacier melt under climate warming then depends on both glacier cover and the distribution of ice with elevation.

Finally, we have shown the importance of including an explicit representation of debris-covered ice in modeling the contribution of glacier runoff to streamflow for some catchments. The explicit representation of debris-covered ice can reduce estimated runoff from glacier melt in situations where debris cover is important and there is still a reasonably large glacier cover (relative to the total catchment area), as in the northern catchment we considered here. More process-based glacier evolution models for large catchments need to be implemented to refine this result, which remains challenging for debris-covered glaciers.

We greatly appreciate and thank Afghanistan government offices that supported this research, special thanks to the Ministry of Energy and Water (MEW) of Afghanistan for providing us with climate and streamflow datasets. This research was supported financially by Swiss Government Excellence Scholarship Program and the University of Lausanne.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.