Increased rainfall-runoff drives flood hazard intensification in Central Himalayan river systems

The Central Himalayan floodplain is one of the most flood-affected regions of the world¹ with the September 2024 floods causing damages of 1% of Nepal’s GDP with 236 deaths and 8,400 people displaced² and by 2050, flood damages are projected to account for 2.2% of annual GPD³. Heavy rainfall and pluvial floods are projected to increase in this region, with changes occurring at faster rates than the global average^4,5. Fluvial flood response to climate change is more uncertain, due to the complex interactions between environmental, climatic and hydrological conditions⁶. Understanding this response and its driving processes is critical for the development of flood risk management strategies targeted at the catchment-scale.

The current understanding of the climate-change impact on Central Himalayan river floods is largely based on global-scale projections of Global Circulation Model (GCM) runoff projections⁷, hydrodynamic routing of GCM runoff projections^8,9, or hydrological modelling using GCM climate projections¹⁰. While global studies are valuable for understanding regional trends, their validity on the catchment-scale is limited for several reasons: (i) the global application leads to a loss of performance at the catchment-scale (e.g. the observed direction of change is not captured by GCM runoff predictions for 40% of the catchments⁷; (ii) the projections are averaged over a coarse modelling grid, which restricts change impact assessments to large catchments⁷; (iii) the small-scale variability of the hydro-climatic system of mountainous catchments is not captured by the coarse modelling grid of GCMs^5,11. Regional-scale assessments provide an improved representation of this variability due to finer modelling grids, but the studies investigating the Central Himalayas have used older generations of climate projections^12,13, and/or do not capture appropriate design floods used for flood risk management purposes (e.g. 1% Annual Exceedance Probability (AEP)). Furthermore, the studies that use Flood Frequency Analysis (FFA) to determine design flood magnitudes do not estimate the uncertainties introduced by the FFA, although this is also needed to correctly design flood mitigation strategies¹⁴.

Here, we present the first, comprehensive regional climate change impact assessment of design flood hazards in the Central Himalayan Karnali River in Nepal and China (Fig. 1) derived from high-resolution, large ensemble modelling using the latest generation of climate projections, and capturing climatic, hydrological and statistical uncertainties. The applied framework couples the climate projections with a continuous, fully distributed (500 × 500 m) hydrological model and a FFA. We apply: probabilistic climate projections (12 bias-corrected and downscaled CMIP6 models¹⁵ for medium (SSP245) and high emission (SSP585) scenarios for the near-future (2020–2059) and far future (2060–2099) to estimate the climatic uncertainty; a Generalised Likelihood Uncertainty Estimation framework (GLUE)¹⁶ with 64 parameter sets to capture the hydrological uncertainty, and utilise the GLUE predictions to capture the measurement uncertainty; and a bootstrapping framework¹⁷ to capture the sampling uncertainty in the FFA. We apply the ‘Spatial Processes in Hydrology’ (SPHY) hydrological model¹⁸ which represents cryospheric and hydrological processes, including the simulation of glacier retreat¹⁹, which improves predictions for long-term climate change assessments. However, the current land cover is maintained throughout the modelling period as no high-resolution projections of potential future land cover exist for the Karnali catchment, and estimating the potential future land cover change would increase the uncertainty in the model predictions. Thus, the flood hazard projections do not account for potential land use changes.

Our simulations find, with high confidence, that flood magnitudes for the 1% AEP will increase by 23% and 26% in the near-future (2020–2059), and increase further in the far-future (2060–2099) by 40% and 79% for the medium and high emissions scenarios respectively. We also show that rainfall-runoff is the largest source of flood water (> 90%), but that snowmelt is an important source for individual medium-discharge flood events. We find that the FFA uncertainty is strongly affected by the timing and occurrence of rare precipitation events, which dominate flood magnitude response. We therefore argue that a large ensemble of climate models is needed to capture a large enough sample size of rare events to enable the analysis of flood frequencies under future climates. Both Atmosphere Only General Circulation Models (AOGCMs) and Earth System Models (ESMs) should be included in the ensemble to capture climate model uncertainty, as AOGCMs tend to project a wider range of changes with respect to the baseline, with lower minimum bounds, whilst ESMs project a narrower range of and higher median changes. These results build on existing studies to provide an advance in understanding of fluvial flood uncertainty and drivers, and propose a modelling framework for predicting AEP flood magnitudes under climate change, necessary for flood risk mitigation.

Projected changes in flood magnitudes

Flood discharge projections obtained from hydrological modelling (Fig. 2a) are used in the FFA to estimate flood magnitudes for probabilistic climate projections. Our results show that in the near-future (2020–2059), the P₅₀ flood magnitudes of events with a 1% AEP are projected to increase from the baseline (1975–2014) scenario by + 23% (P₁₀ + 21%, P₉₀ + 17%) for the medium-emission scenario, SSP245, and + 26% (P₁₀ + 25%, P₉₀ + 27%) for the high-emission scenario, SSP585 (Fig. 2b), with 10 (SSP245) and 11 (SSP585) of the 12 climate models projecting increases. However, the ensemble-median difference between both scenarios is 3%, with higher increases for SSP245 than for SSP585 for 6 of the 12 climate models.

In the far-future (2060–2099), the P₅₀ 1% AEP flood magnitude increases from the baseline by 40% (P₁₀ + 33%, P₉₀ + 48%) for SSP245 and 79% (P₁₀ + 54%, P₉₀ + 82%) for SSP585 (Fig. 2b) with 12 (SSP245) and 11 (SSP585) of the 12 climate models projecting increases. Similarly, the frequency of a given flood magnitude increases, with the P₅₀ 1% AEP event projected to occur once every 11 years for SSP245 and every 3 years for SSP585, with 12 (SSP245) and 11 (SSP585) of the 12 climate models projecting increases. An intensification of magnitude occurs with higher emissions, with 11 of the 12 climate models projecting higher increases for SSP585 than for SSP245.

The timing when both emission scenarios diverge varies between the percentiles of the ensemble. The P₅₀ 30-year-mean flood discharge increases at similar rates for SSP245 (+ 20%) and SSP585 (+ 21%) between the years 2020 and 2059 (Fig. 2c). In the year 2060, the P₅₀ 30-year-mean flood discharge of both scenarios starts diverging, and it further increases by 7% (SS245) and 35% (SSP585) between 2060 and 2099. The timing when higher flood discharges start to be predicted for SSP585 varies with the percentile, the P₉₀ starts diverging between 2035 and 2040, the P₅₀ around the year 2060, and the P₁₀ between 2065 and 2070.

We do not identify systematic patterns in the projected floods for related climate models²⁰, despite models in the same family usually having similar climate properties; models of the same family can project both small and large changes in 1% AEP flood magnitudes (Fig. 3). Some models (e.g. HadAM, CCM) project smaller spreads, whilst others project larger spreads (e.g. ECMWF). Nevertheless, we find that ESMs generally project larger increases and provide a smaller range of change than AOGCMs, possibly related to their inclusion of the carbon cycle, vegetation, atmospheric chemistry and other processes, although other comparisons accounting for differences between AOGCMs and Earth System Models (ESMs) are currently lacking²⁰, and further testing is necessary to confirm this result.

The sources of floodwater

We analyse the floodwater sources simulated by the hydrological model, divided between direct rainfall-runoff, glacier melt, groundwater-driven baseflow and snowmelt, to understand the drivers of the projected flood intensification. In the baseline period, rainfall-runoff is the most important source, with a mean floodwater contribution of 78 ± 11% (± standard deviation), followed by glacier melt (10 ± 5%), baseflow (8 ± 5%) and snowmelt (3 ± 8%) (Fig. 4). The importance of rainfall-runoff increases with the flood discharge, with rainfall-runoff contributing ≥ 90% to events with discharges ≥ 17,500 m³/s (high flows) (Fig. 5a). Consequently, the importance of the other sources decreases with increasing discharge. However, other sources are important for events with lower (< 10,000 m³/s) and medium (10,000–17,500 m³/s) discharges with contributions ≥ 25%, and 10–25% respectively (Fig. 5a). High rainfall-runoff contributions are predicted without clear seasonal patterns, but event timing is important: higher snowmelt contributions are predicted in the pre- and early-monsoon seasons (typically up to the 27th week of the calendar year); glacier meltwater contributions increase towards the mid-monsoon season (weeks 27–35); and higher baseflow contributions are predicted in the late-monsoon and post-monsoon seasons (week 35 onwards). These results provide a new contribution to the literature.

The drivers of change

In the near-future and far-future, the projected flood intensification is driven by increasing rainfall-runoff for both scenarios. The mean flood discharge increases by 1,380–4,581 m³/s compared to the baseline, whereas the mean rainfall-runoff increases between 1,202 and 4,255 m³/s, and thus accounts for 91–93% of the projected increase of the mean flood discharge (Fig. 4). Furthermore, increasing rainfall-runoff also causes the intensification of high flood discharges, which remain composed of ≥ 90% rainfall-runoff (Fig. 5a). The baseflow and glacier melt contributions to the mean flood discharge increase by 77–215 m³/s and 82–215 m³/s, respectively (Fig. 4). Their relative importance decreases due to the higher rainfall-runoff increase, and we project no changes in their seasonality. The flood discharge intensification is dominated by increasing rainfall-runoff in both emissions scenarios, and thus, differences in future emissions affect the rate of change rather than their driving processes.

Temperature-related changes in melt processes lead to a decrease in snowmelt runoff for most flood events. The mean snowmelt contribution decreases from 3 ± 8% in the baseline to 2 ± 7% in the near-future (SSP245 and SSP585) and 1 ± 6% (SSP245) and 0 ± 2% (SSP585) in the far-future, which is caused by an earlier onset of the melting season (Fig. 5b). However, the occurrence probability of medium-discharge flood events with 10–40% snowmelt contributions is projected to increase from 0.30% in the baseline to 0.65% in the near-future, after which it decreases to 0.17% for SSP585. For SSP245, an increase to 0.45% is delayed to the far-future. The changes in the frequency of these events indicate that rising temperatures cause an initial increase in melt contributions for individual pre-monsoon and early monsoon season flood events. The warming also shifts the melting season earlier, which leads to lower flood discharge contributions as the temperature continues to increase, as indicated by the lower frequency in the far-future of SSP585 (Fig. 5a). While the occurrence probability is below 1% at the mountain outlet, the projected changes may indicate a shift in the flood seasonality, frequency and magnitude in smaller upstream subbasins.

Performance of the hydrological model

The hydrological ensemble is a behavioural representation of the hydrological system of the Karnali River catchment, indicated by the hydrograph inspection, the model efficiencies and the streamflow composition. The ensemble replicates the flow seasonality with high flows during the monsoon and low flows during the winter and times the transitions in the Pre-monsoon and Post-monsoon seasons well (Fig. 6). The seasonality of the runoff composition (Supplementary Figure S4) with snowmelt contributions in the Pre- and early Monsoon season (Mar – Jun), rainfall-runoff and glacier melt during the Monsoon season (Jun – Sep), and baseflow in the Post-monsoon season (Oct – Nov) reflects our understanding of Central Himalayan catchment hydrology^11,21,22,23. The good model performance is, furthermore, illustrated by the high ensemble median efficiencies of: (i) the Nash-Sutcliffe-Efficiency (NSE) of 0.85 (range: 0.75–0.87) for the calibration period and 0.82 (range: 0.70–0.85) for the validation period, and (ii) the coefficient of determination (R²) of 0.85 (range: 0.77–0.87) and 0.84 (0.80–0.86) for the calibration and validation periods, respectively. These efficiencies are at the higher level of reported efficiencies in Central Himalayan catchments, ranging between 0.56 and 0.88 (NSE) and 0.63–0.89 (R²)^{23,24,25,26,27}.

The evapotranspiration module was calibrated and validated against satellite-based annual actual evapotranspiration (ET_a) estimates (Fig. 7a). The ensemble-median PBIAS is -11% (-20% to + 10%) in the calibration period and − 13% (-22% to + 7%) in the validation period, indicating a moderate underestimation of the catchment-mean ET_a rates. However, the performance of each ensemble member is satisfactory or better²⁸, and the correlation between observed and simulated land-cover-specific ET_a rates is very high (R² > 0.95), indicating a very good representation of the spatial ET_a variability (Fig. 7b). However, the model’s ability to reproduce ET_a rates varies between the years, and the percentage difference is within − 3% to -15%, but reaches up to -22% to -32% in the years 2005, 2014 and 2015. However, given the large uncertainty in the climate input datasets, these biases are likely caused by precipitation input deficits, as precipitation products are known to underestimate the high-mountain precipitation^29,30. If Et_a observations are not considered in the model calibration, precipitation input deficits may be compensated by an ET_a underestimation^31,32. Thus, we argue that it is important to consider ET_a observations during the calibration of hydrological models in Himalayan catchments.

The uncertainty in the modelling ensemble

We analyse the uncertainty composition of the ensemble predictions using the 3-way analysis of variance (ANOVA) framework³³. However, instead of decomposing the uncertainty contribution of climate models, climate scenarios and the internal variability, we assess the contribution of the ensembles of climate models, hydrological parametersets, and flood frequency curves. It is worth noting that the FFA ensemble estimates the uncertainty associated with the internal variability from a bootstrapping approach¹⁷.

The hydrological ensemble is the largest uncertainty source in the baseline, contributing 52% to the overall variance of the 1% AEP event (Fig. 8). For the projected scenarios, the hydrological uncertainty varies between the climate scenarios, contributing 53% (SSP245) and 43% (SSP585) to the 1% AEP in the near-future, and 37% (SSP245) and 35% (SSP585) in the far-future. However, the mean standard deviation (σ) of the hydrological ensemble averaged over the climate ensemble remains with 20–22% stable for all climate scenarios and return periods (Supplementary Figure S 5). This indicates that the hydrological parametersets introduce a large, but relatively constant uncertainty, and that the decreasing relative variance contribution reflects an increase in the overall variance.

The variance introduced by the climate models increases with time and emissions, but this trend is superimposed by an increasing FFA variance for higher frequencies. The climate model ensemble contributes 30% to the overall variance of the 1% AEP in the baseline. For SSP245, this fraction decreases to 26% in the near-future, after which it increases to 36% in the far-future. For SSP585, the contribution is with 41% (near-future) and 49% (far-future) notably higher than for SSP245, which indicates that the flood response becomes increasingly uncertain with time and emissions. The initial decrease in the climate ensemble variance for the near-future of SSP245 is caused by superimposing effects of the FFA variance, which is indicated by consistent climate variance increases with time and emissions for lower frequencies 10% AEP to 2% AEP).

The uncertainty associated with the FFA ensemble is the lowest of the three components, but it interacts with the climate model uncertainty. Generally, the FFA uncertainty increases with the frequency, because the FFA is is sensitive to the slope of the flood frequency curve above the 5% AEP, which is fitted to a low sample size and extrapolated beyond the 40-year record length to the 1% AEP^34,35. In the baseline, the FFA variance increases with the event frequency from 8% (10% AEP) to 11% (1% AEP). For the projections, the FFA variance patterns vary between the scenarios, whereas the 1% AEP variance decreases to 8–9% for SSP585, and increases to 12% (near-future) and 15% (far-future) for SSP245. Generally, the variance of component interactions is less than 3%, but in the case of SSP245 (far-future), the interaction effect of climate and FFA ensembles is 7% (p-value < 0.1).

We attribute this CM-FFA interaction effect, and the higher FFA variance to the flood discharge time series projected for individual climate members. The Flood Percentile (FP_n) aggregates the flood events by the n^th frequency percentile (see Supplementary Information 1). In the case of the far-future of SSP245, 5 out of 12 models predict > 40% higher FP₁₀₀ discharge (highest discharge in the record) relative to the FP_92.5 discharge (4th highest discharge) (Fig. 9). This large difference between the rarest and more frequent events increases the sampling uncertainty because the rarest events are removed or duplicated in the bootstrapping process, which increases the difference between the slopes of the individual flood frequency curves. The case of the far-future of SSP245 indicates that the FFA uncertainty is affected by projections of rare extreme precipitation events and is thus sensitive towards the selection of climate models and the time frame in which the most extreme events may or may not occur.

Current climate projections indicate with high confidence that heavy precipitation and pluvial floods increase with increasing Greenhouse Gas (GHG) emissions in the Central Himalayas, but are less certain how these increases affect fluvial floods^4,6. Previous studies suggest an increase in fluvial flood discharges, and while direct comparisons are hampered by differences in the research design (e.g. spatial and temporal extent, datasets, models, catchments, event frequencies), our projections fall within the reported range of + 20 to + 108%^9,10,12,13. We can attribute these changes to increasing rainfall-runoff, consistently predicted for a large ensemble of climate projections and hydrological model configurations, suggesting that projected increases in heavy precipitation will cause an intensification of river floods in the foreland of Central Himalayan catchments.

The uncertainty in predicting the hydrological response to climate change can exceed the uncertainty in predicting climate response to GHG emissions in mountainous catchments in which melt processes are important³⁶. In the Himalayas, the uncertainty in hydrological modelling applications is increased by the coarse resolution of physical and climatic datasets and a lack of hydro-meteorological observations at higher elevations (> 4,000 masl)^5,23,29,37. Here, the hydrological modelling is the largest source of uncertainty in the baseline. However, the ‘good’ to ‘very good’ modelling efficiencies²⁸, the agreement in the seasonality of the runoff composition with other studies^11,21,22, and the similarities in the projected changes in the hydrological behaviour, such as an earlier onset of the snow-melt season^11,38,39, indicate that the hydrological ensemble reflects the catchment hydrology well. The relative uncertainty in the hydrological predictions decrease for future climates, which may be attributed to the dominance of rainfall-runoff in floodwater generation, as the rainfall-to-runoff conversion is less complex than other processes that include freezing, melting, and percolation.

This study does not consider projected land use changes in the experimental design, although the land use is an important component of the hydrological cycle and affects flood hazard characteristics⁴⁰. In the case of Nepal, two drivers of land use change may be of particular relevance: (i) climatic changes shift climate zones and, thus, the natural vegetation cover^41,42, and; (ii) the expansion or abandonment of agricultural terraces (e.g. caused by migration), which affect the surface and subsurface runoff generation and the hydrological connectivity and, thus, alter the contribution to the flood discharge^43,44.

The Himalayan foreland in Nepal is, under both past and current climatic conditions, prone to flooding, with high mortalities and impacted populations, economic damages, and other repercussions, including increasing food insecurity and the outbreak of epidemics^1,45,46,47. Future GHG emissions are likely to exacerbate these impacts. For the near-future, the relationship between flood magnitude and emissions is less pronounced due to the small differences between the SSP245 and SSP585 flood magnitudes from the relatively small differences in GHG emissions between the two climate scenarios (also possibly due to internal variability (i.e. natural climatic fluctuations) obscuring responses^48,49. For the far-future however, our simulations suggest that GHG emissions cause a sustained flood hazard intensification – even for the SSP245 scenario, where emissions are higher for each year of the near-future than for the far-future, while the 30-year mean flood discharge continuously increases until the year 2099. The flood projections imply that: (i) Central Himalayan floods will intensify with time and emissions; (ii) this intensification is likely to continue for decades after the peak emissions and; (iii) flood magnitudes are likely to remain above current levels until the end of the century.

Flood magnitudes in the far-future are characterised by an increased spread of flood magnitude predictions between the individual CMIP6 members, ranging from + 9% to + 69% (SSP245) and − 12% to + 144% (SSP585) for the median 1%-AEP projections from the hydrological model and FFA ensemble. Furthermore, the FFA projections are sensitive to the slope of the frequency curves, and thus to the rarest, most extreme events (i.e., the FP₁₀₀) and the discharge difference between the rarest and more common events (i.e., the FP_92.5 and FP₁₀₀). Both the climate model differences and the FFA sensitivity to rare events make the projection of design floods (i.e. the 1%-AEP event) susceptible to the unknown true frequency of the rarest projected events. Therefore, it is important to maintain a large ensemble of climate models to guarantee a larger sample size of the rarest events (i.e. the FP₁₀₀ events of the climate models). If necessary, a new selection criterion should be developed for sub-setting climate models, which account for the frequency of extremes rather than the commonly used temperature and annual precipitation subset (e.g. dry-cold, dry-warm, wet-cold, wet-warm). Furthermore, a subset should include a combination of AOGCMs and ESMs, as these tend to project different rates of change.

The near-future flood projections indicate that the Karnali River is already in a phase of flood hazard intensification, so existing flood risk management strategies need to be adapted to maintain current flood protection levels. Flood defences, such as embankments and outlet structures, are typically designed and implemented on decadal time scales and, therefore, new projects should take the projected future changes in flood hazards into consideration^50,51. The inherent uncertainties of future emissions and their impact on the climate and flood response create conflicts between designing protection levels for the worst-case scenario versus the economic feasibility of such structures. New, flexible, and anticipatory flood risk management strategies are therefore needed that can account for emerging new knowledge and understanding. Such strategies could include flexible designs for future flood risk infrastructure, which include the potential for later adaptations and upgrades, such as increases in storage space for floodwater, or increased dam walls levels. Additional strategies could include allocating space now for flood water retention areas within urban-zoning development plans, and policies and institutional workflows that reduce the period between planning and operation.

This study presents an environmental modelling framework that uses high-resolution, large ensemble modelling to predict future design floods at the catchment scale for the Central Himalayan Region. We provide design floods (i.e. 1%-AEP) for future climate change scenarios, and the hydrological processes driving these potential changes, thereby providing practical information for climate change adaptation for flood risk management. The framework extends current uncertainty analyses by incorporating sampling uncertainty into the FFA. The uncertainty analysis highlights the sensitivity of the FFA to the flood discharge projections of individual climate models, emphasising the necessity to maintain large climate ensembles for flood impact analyses. These ensembles should contain both AOGCMs and ESMs, as these families differ in their flood projections.

The application of this framework projects the intensification of flood hazards over time and with increasing emissions. The hydrological and climate ensemble consistently predicts that the flood hazard intensification is driven by rainfall-runoff increases, which increases the evidence that increases in heavy precipitation will lead to increased river floods in the foreland of Central Himalayan catchments.

Research design

We combine hydrological and statistical modelling with probabilistic climate projections to predict the potential changes in future flood hazards in the Karnali River in Nepal and China. In the first stage, we implement the hydrological model Spatial Processes in Hydrology (SPHY) for historical conditions (past and observed hydro-meteorological conditions) and identify an ensemble of behavioural parameter sets. This ensemble is then forced with probabilistic climate projections for three scenarios: the Baseline (1975–2014), the medium-emission scenario SSP245 (2020–2099) and the high-emission scenario SSP585 (2020–2099). In the last stage, a Flood Frequency Analysis (FFA) is conducted with the simulated flood discharges, which are classified by the Annual Maximum flow (AMAX). We compare the differences between the Baseline and the projected scenarios (SSP245 and SSP585) to quantify the changes in flood hazard magnitude and frequency between the past (1975–2014) and the near-future (2020–2059), as well as the far-future (2060–2099). We chose 40-year time frames over the commonly used 30-year time frames to increase the length of the flood record, which is the main source of uncertainty in the FFA³⁴.

Data processing

The utilised datasets comprise gridded and point-based datasets that determine the environmental (static) and climatic (dynamic) boundary conditions, plus calibration and validation datasets for the hydrological model (Supplementary Table 1). All gridded datasets are resampled to the modelling resolution (500 × 500 m) and reprojected to the UTM44N coordinate system (EPSG: 32644). The preprocessing of hydro-climatic data includes:

• Precipitation (historical): The network of precipitation gauges is concentrated in the floodplains and valleys, and precipitation at high elevations is poorly observed²⁹. We therefore use monthly gridded satellite-based precipitation estimates (GPM IMERG Final Precipitation L3 1 Month V006⁵²) and disaggregate them to the daily resolution using precipitation measurements of the Department of Hydrology and Meteorology, Nepal (DHM), following the approach of Arias-Hidalgo et al. (2013)⁵³.

• Temperature (historical): We use the WATCH forcing Data methodology applied to ERA-Interim dataset (WFDEI)^54,55 to determine the historical temperature boundary conditions. This dataset has a spatial resolution of 0.5 × 0.5° and is downscaled to the modelling resolution using a lapse rate derived from DHM temperature observations.

• Climate projections: We use the temperature and precipitation predictions of an ensemble of 12 CMIP6 Global Circulation Models (GCM), which have been bias-corrected using Empirical Quantile Mapping and downscaled to 0.25 × 0.25º by Mishra et al. (2020)¹⁵. The CanESM5 member was removed from the ensemble because it does not depict the monsoon seasonality in the catchment (Supplementary Fig. 3). The temperature data is downscaled to the modelling resolution using a lapse rate, and the precipitation data is downscaled using bilinear interpolation.

• Discharge: We use daily discharge data and stage-discharge observations (DHM) at the mountain outlet for the model calibration. The uncertainty in the discharge is estimated using the BaRatin approach⁵⁶.

The hydrological model

The catchment hydrology is simulated by the Spatial Processes in Hydrology (SPHY) model version 3.0¹⁸. SPHY is a distributed (raster-based) model that simulates soil processes, snow processes and evapotranspiration on the grid scale, and glacial processes on a sub-grid scale. The model distinguishes four types of runoff; rainfall-runoff Q_RR, baseflow Q_BF, snowmelt Q_SM, and glacier melt Q_GM which compose the total runoff Q_TOT of each cell:

The total runoff is then routed to the downstream cell based on a flow direction raster. The model implements a flow recession coefficient to account for channel friction (Supplementary Table 2).

The soil processes are calculated as in the SWAT and SWAP models⁵⁷. The soil is divided into a rootzone layer SW₁, a subzone layer SW₂ and the groundwater layer SW₃. The upper two layers SW₁ and SW₂ contribute to rainfall-runoff whereas the model simulates surface runoff (saturation excess runoff and infiltration excess runoff), and slower lateral flow from the soil layers. Percolation to the groundwater storage and baseflow are calculated using a recession coefficient to relate the baseflow to the groundwater recharge.

The reference evapotranspiration is estimated from the modified Hargreaves evapotranspiration equation⁵⁸ and crop coefficients are used to estimate the potential evapotranspiration rates for different land covers⁵⁹. Snow processes are simulated by a dynamic snow storage model⁶⁰ using a degree-day-factor approach to quantify the snow melt⁶¹. Glacier melt is also calculated using the degree-day factor approach whereas separate factors are defined for clean-ice and debris-covered glaciers. The evolution of glaciers (the retreat or advance) is simulated from the mass balance (snow accumulation – glacier melt) using a mass conserving ice redistribution approach which is applied once per year¹⁹.

Hydrological modelling setup

We use the SPHY model version 3.0 with the dynamic glaciers, infiltration, groundwater and snow modules. The model is applied at 500 × 500 m grid resolution with daily time steps. It is calibrated and validated for 2002–2006 and 2007–2015, respectively.

A stepwise calibration approach is used to identify behavioural parameter sets. In the first stage, a Regional Sensitivity Analysis with 1,500 Latin Hypercube Samples (LHS)⁶² of 21 parameters is conducted to identify sensitive parameters. In the second stage, the 15 sensitive parameters are calibrated within a GLUE framework using 10,000 parameter sets sampled with LHS. We implemented two scaling factors as calibration parameters: (i) a precipitation correction factor to account for the bias in the precipitation data (similar to Lutz et al. (2014)¹¹, and (ii) a multiplication factor for the crop coefficients which control the evapotranspiration for different land cover classes (similar to Khanal et al. (2021)¹⁹. The precipitation correction factor is used for the model calibration but is not transferred to the climate projections because the CMIP6 data is bias-corrected.

We focus on the model performance during high flows for selecting the behavioural parameter sets using a multi-criteria approach to increase the robustness of the ensemble predictions⁶³. However, we exclude parameter sets with percentage BIAS (PBIAS) exceeding ± 30% for the annual actual evapotranspiration from MODIS, the 8-day snow extent from MODIS, and the discharge during the winter months to prevent the selection of parameter sets that do not depict critical aspects of the catchment hydrological behaviour. Of the remaining parameter sets, we select the 25 members with the highest Nash-Sutcliffe Efficiency (NSE), PBIAS of high flows (≥ 5,000 m³/s), and a modified version of the extended GLUE (eGLUE)⁶⁴ to account for the uncertainty in the discharge observations. The eGLUE quantifies the number of days during which the simulated discharge is within the uncertainty interval of the observed discharge. The modified eGLUE is the sum of the simulated discharge for all days that fall within the observed confidence interval divided by the total number of days. This adjustment is made to prevent lower flows, which occur throughout most of the year, from dominating the classification. After the removal of duplicates, 64 parameter sets are maintained.

Flood frequency analysis

We apply a Flood Frequency Analysis (FFA) to quantify the changes in flood magnitudes for the projected climate scenarios. The Flood Frequency Analysis (FFA) is a statistical method that utilises an Extreme Value Distribution (EVD) to estimate flood magnitude based on flood frequency (occurrence probability) from a record of flood events⁶⁵. The FFA is fitted to each AMAX record of the modelling ensemble (12 CMIP6 members X 64 hydrological models) using the L-Moments approach for parameter estimation⁶⁶. We evaluated the performance of eight EVD because it is generally impossible to determine the most suitable distribution before the application⁶⁵, and selected the Wakeby distribution for the Flood frequency analysis (Supplementary Table 4). We employ a bootstrapping approach (n = 1,000) based on Burn (2003)¹⁷ to estimate the sampling uncertainty, which relates to the difference between the distribution used in the FFA and the true distribution of the flood discharge. The flood discharges are sorted before the bootstrapping, and the same bootstrap samples are taken for all combination of climate models and hydrological parametersets to ensure the comparability of FFA samples in the uncertainty decomposition.