2.1 Hydrological model

The large-scale hydrological model Community Water Model (CWatM) is an open-source model developed by the International Institute for Applied Systems Analysis (IIASA) (Burek et al., 2020). The model includes all dominant hydrological storages and processes, such as snow, soil, groundwater and river routing. Additionally, it can simulate irrigation, industrial, domestic and livestock human water use and reservoir regulations. For an in-depth model description, we refer to Burek et al. (2020). CWatM is used at 5 and 30 arcmin resolution but can also be used at 1 km resolution (Guillaumot et al., 2022). At 30 arcmin resolution, CWatM was used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) rounds 2b and 3, which compares large-scale hydrological model outputs (Frieler et al., 2017; Warszawski et al., 2014). CWatM has extensive and publicly available documentation of the source code, the model structure, and model training and tutorials (https://cwatm.iiasa.ac.at/, last access: 26 June 2024). Together with its modular structure, this makes CWatM a highly suitable candidate for the model coupling. CWatM has a simple glacier representation: at high-elevation zones where snow has not melted until the beginning of summer, the melting rate is increased, and the average temperature of the grid cell is applied. This implicitly simulates the downward movement of glaciers (Burek et al., 2020). However, this glacier representation has several shortcomings. First, it assumes that glaciers are in equilibrium, meaning that only snow previously accumulated is melted. This representation neglects the change in glacier-sourced runoff when glaciers retreat, a phenomenon that is ongoing (Hugonnet et al., 2021) and projected to continue throughout this century (Rounce et al., 2023). Second, it assumes the presence of glaciers wherever snow accumulates in the model. We argue that this glacier representation is more akin to a snow redistribution function used to avoid excessive snow accumulation (Freudiger et al., 2017; Burek et al., 2020).

2.2 Glacier model

The Open Global Glacier Model (OGGM) framework is a community-driven, well-documented and open-source global glacier model (Maussion et al., 2019). This makes it stand out from a number of global glacier models which are used for global glacier mass change projections (Marzeion et al., 2020). OGGM has been used in global studies (Marzeion et al., 2020; Rounce et al., 2023) but also on regional and basin scales (e.g. Tang et al., 2023; Furian et al., 2022; Yang et al., 2022). OGGM has a modular structure, which is beneficial for new developments and a variety of research questions requiring different model setups.

Each glacier is simulated individually by OGGM. For this study, OGGM used the glacier outlines from the Randolph Glacier Inventory (RGIv6.0; Pfeffer et al., 2014) and elevation-band flowlines (e.g.  Huss and Farinotti, 2012; Werder et al., 2020) derived from a digital elevation model. A one-dimensional shallow-ice flowline model of OGGM approximates the influence of ice flow on glacier dynamics. A temperature index mass-balance model computing mass balances at each elevation band was calibrated for every glacier by glacier-wide geodetic observations (2000–2019, Hugonnet et al., 2021). The ice volume consensus estimates of Farinotti et al. (2019) at the RGI year (often near the year 2000) were matched regionally by calibrating the ice creep parameter A. The model version does not consider refreezing.

For a detailed model description, we refer to Maussion et al. (2019) and the online documentation (https://docs.oggm.org/, last access: 26 June 2024). The model framework is under continuous development and now possesses the capability to simulate at a daily resolution using a daily mass-balance model and daily climate data, a novel development for large-scale glacier models (Schuster et al., 2023b). In the context of hydrological applications, it is relevant to mention that the utilized OGGM version (v1.5.3) does not differentiate between snowmelt and ice melt on glaciers.

2.3 Model differences in the representation of elevation, precipitation and snowmelt

Both models were slightly adapted to facilitate coupling: in CWatM, we incorporated the rain–snow partitioning function from OGGM. OGGM was modified to provide daily outputs for melt and rainfall on glaciers. Nevertheless, differences between the two models in the representation of snow accumulation and melt processes remain (Table 1).

An important difference is the spatial resolution of snow accumulation and melt processes. CWatM has a maximum of 10 evenly large elevation zones per grid cell (at 5 arcmin resolution, one elevation zone per ∼ 10 km²) for which snow processes are modelled individually. In OGGM, the snow and ice processes are modelled on a much finer scale. Each glacier is split into elevation zones of 30 m height difference. Moreover, the spatial resolution of the model output is different, as OGGM outputs are given per glacier, whereas CWatM inputs and outputs are given per simulation grid cell.

Another relevant difference between the models is the correction of precipitation input data. Precipitation amounts are often underestimated in mountain areas, but the magnitude of this error remains difficult to assess (Beck et al., 2020; Azam et al., 2021). Glacier models use a precipitation factor to correct for measurement errors and additional errors in snow input, such as those caused by avalanches, thereby allowing for better alignment with mass-balance data (Hock et al., 2019; Rounce et al., 2020). In CWatM, a multiplicative factor for snowfall correction is implemented as an optional calibration parameter following the rationale that snowfall has a larger error range than rainfall (Beck et al., 2020). Following the same rationale, we applied the precipitation factor in OGGM only to snow and not to rain. However, efforts to further harmonize the two models by using the same precipitation correction factor resulted in unrealistic calibration and a deterioration in model performance. Therefore, we decided to retain the snow correction procedures of each model, although this leads to disparities between the two models. Thus, the precipitation factor in OGGM was calibrated to a fixed value across all glaciers in the studied region to best represent the observed interannual mass-balance variability from WGMS (2020), which is the standard procedure in OGGM v1.5.3. In CWatM, the snowfall correction factor was only calibrated when there was clear evidence of precipitation underestimation in the river basin. In other cases and for global simulations, no snowfall correction was applied. This leads to larger precipitation inputs in the coupled model compared to the baseline model, a topic further discussed in Sect. 6.2.1.

Both models use a temperature-index model to represent snowmelt and ice melt, with CWatM including melt seasonality and the effect of increasing melt during rainfall. Temperature-index models, which use a degree-day factor, i.e. snowmelt coefficient, are widely used in hydrological models to calculate snowmelt. They have low data requirements, which is beneficial for large-scale studies (Girons Lopez et al., 2020). The utilized OGGM version (v1.5.3) uses one combined degree-day factor for snowmelt and ice melt. This parameter was calibrated per glacier using 20-year average geodetic mass-balance data from Hugonnet et al. (2021). However, in CwatM the degree-day factor was calibrated jointly with other model parameters using discharge data (Table 1), resulting in one parameter value for the whole basin area upstream of the gauge used for calibration.

Table 1Overview of model inputs, resolution and outputs of the applied CWatM and OGGM versions and their configurations.

Download Print Version | Download XLSX

Figure 1Schematization of the steps to sequentially couple the glacier model (OGGM) and the hydrological model (CWatM).

2.4 Model coupling approach

The main rationale behind the coupling was to adhere to the modular philosophy of the models and maintain an open framework to enable their continuous development. Therefore, it was decided to implement a one-way sequential coupling where OGGM was run first, and the output was translated via a pipeline consisting of a set of functions to input into CWatM (Fig. 1). In this manner, the glacier output could potentially also be used with other hydrological models that use NetCDF files as input files. OGGM handled the glacierized portions of the modelling domain, while CWatM handled the non-glacierized portions. The glacier area, which was updated each year, determined the extent of the model domain addressed by OGGM and therefore excluded from CWatM simulations. In the modelling chain, the following sequence was followed: first, the glaciers within the modelled river basin were identified. Second, OGGM was calibrated and run for these glaciers. Third, the pipeline translated output from OGGM into a format readable as input for CWatM. Finally, CWatM was run with glacier area and glacier-sourced melt as input data. The coupling is explained in more detail in the following sections.

2.4.2 Translate OGGM output into CWatM input

CWatM requires input data per grid cell in NetCDF format (latitude, longitude, time dimensions). Therefore, OGGM output per glacier had to be translated into CWatM input per grid cell (Fig. 1 step 3). The following steps were used to generate raster NetCDF files of yearly glacier area and daily glacier-sourced melt and rain on glaciers.

The rain and melt amounts over a glacier were routed to the grid cell in which the terminus of the glacier is located at an RGI date. The glacier-sourced melt includes snowmelt and ice melt on the glacier area. The glacier-sourced melt and rain on all glaciers with a terminus location within the same grid cell were summed to obtain the total glacier-sourced melt (M_t,i) and rain (R_t,i) on a glacier per grid cell (Eqs. 1 and 2).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{M}_{t,i}=\sum _{g=\mathrm{1}}^G{M}_{g,t,i}\text{(2)}& {\displaystyle }& {\displaystyle }{R}_{t,i}=\sum _{g=\mathrm{1}}^G{R}_{g,t,i}\end{array}$$

One glacier can cover several grid cells, and one grid cell can contain several glaciers. Therefore, the geometric glacier outlines of RGIv6.0 were intersected with the model grid at WGS84 and then reprojected to Eckert IV to calculate the glacier area per grid cell at an RGI date similarly to Li et al. (2021). To derive the glacier area of each glacier g at time t in grid cell i ($A_{g,t,i}$), the relative change in glacier area was multiplied by the glacier area contained in the grid cell (Eq. 3). Afterwards, the total glacier area within one grid cell was calculated for every year (Eq. 4). In the case that it exceeded the area of the grid cell due to glacier growth, the surplus area was redistributed evenly to the four neighbouring grid cells.

In reality, the glacier area decreases the most close to the terminus of the glacier, when the glacier melts, whereas during glacier growth, the volume at the top of the glacier increases and subsequently the area at the terminus. However, the coupling was designed for large-scale modelling. Therefore, these small-scale dynamics were neglected when adapting the CWatM modelling domain. We assumed the same relative area reduction or growth for all grid cells that a glacier covers (Fig. S1 in the Supplement).

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{A}_{g,t,i}=\left({\displaystyle \frac{A_{g,t}}{A_{g,\text{RGI}}}}\right)\cdot A_{g,\text{RGI},i}\text{(4)}& {\displaystyle }& {\displaystyle }{A}_{t,i}=\sum _{g=\mathrm{1}}^G{A}_{g,t,i}\end{array}$$

3.1 Calibration of selected river basins (5 arcmin)

To evaluate the model coupling on a river basin scale, four large river basins with glacial influence in Europe (Rhine, Rhône, Glomma) and North America (Fraser) were used (Fig. 2). These river basins are well monitored and have long-term observed discharge time series available at an upstream and downstream gauge. The percentage of glacierized area is the highest for the Fraser (1.2 %), followed by the Rhône (1 %), Glomma (0.7 %) and Rhine (0.2 %) for the region upstream of the downstream gauge.

A genetic algorithm (NSGA-II from the Python package DEAP, Fortin et al., 2012) was used to calibrate 12 parameters in CWatM that govern major hydrological fluxes (Burek et al., 2020, Table 2). These include factors to adjust the evapotranspiration, soil depth, preferential flow to groundwater, soil infiltration and groundwater contribution to streamflow. The parameters were calibrated separately for the upstream and downstream regions of each basin using two gauges (see Fig. 2). Discharge data were obtained from GRDC (GRDC, 2022) and HydroPortail (HydroPortail, 2022). As the objective function, the non-parametric version of the Kling–Gupta efficiency (KGE) was used, which jointly evaluates the errors in the mean, the variability and the dynamics of the discharge (NPE; Pool et al., 2018). It was combined with a penalty for the snow cover error to reduce multi-year snow accumulation (weight 0.8 KGE and 0.2 penalty; see Sect. 2.1 in the Supplement).

CWatM_base and CWatM_glacier were both calibrated for each study basin for the 2000–2009 time period to evaluate the effect of explicit glacier inclusion in the model. Each setup was calibrated five times to capture parameter uncertainty. The performance of the five calibrated parameter sets per river basin was evaluated for the time period prior to calibration (1990–1999).

3.2 Global parameter sets (30 arcmin)

To evaluate the model coupling on a global scale, global parameter sets were used. Most large-scale hydrological models employ one single parameter set for the entire globe and do not undergo an explicit calibration (Telteu et al., 2021). To consider parameter uncertainty we used six parameter sets obtained by calibration of CWatM_base, including the parameter set used in ISIMIP 3 simulations of CWatM. For calibration, we used 601 GRDC stations which complied with criteria laid out in Burek and Smilovic (2022) and had at least 5 years of discharge data in the calibration period. The period for calibration was 1985 to 1994 to maximize the number of discharge stations available. The 12 parameters were calibrated simultaneously to achieve the best possible performance over all stations, similar to Greve et al. (2023). The performance was evaluated for the period of 2004–2013 to maximize the discharge observation length available per river basin. To limit discharge differences between the model setups to the effect of glacier inclusion, the global calibration procedure was not repeated for CWatM_glacier, and the parameters of CWatM_base were used. Repeating the global calibration for CWatM_glacier would only have a minor influence on our results because the weight of glacier-influenced discharge gauges in the global calibration is small.

3.3 Simulations and analysis

CWatM was run from 1990–2019 and for future projections from 2020–2099 using the calibrated parameter sets. A total of 5 past simulations and 50 future simulations (using 5 GCMs and 2 SSPs) were performed for each CWatM_base and CWatM_glacier for the individual river basins. For global simulations, the same strategy was used. For the analysis, the periods of 1990–2019 and 2070–2099 were used because future changes were more pronounced at the end of the century.

Differences in past and future discharge between the baseline and the coupled model were evaluated close to the outlet of the river basins to assess the effect of explicit glacier inclusion on past and future discharge simulations using the ensemble mean values of the calibrations. In addition to projecting future discharge, climate impact studies often assess future changes in discharge by comparing future mean discharges to mean past discharges (e.g. Schewe et al., 2014; Gosling et al., 2017). This quantification of future change is essential for facilitating adaptation. Consequently, we also conducted comparisons between the model setups to assess how the representation of glaciers influences our projections of future changes in discharge.

Figure 3Daily mean hydrographs of the evaluation period (1990–1999) for the downstream station of the study basins. As inset the NPE and its error components over the whole evaluation period are given for CWatM_base (index b) and CWatM_glacier (index g) (α_NP – variability, β_NP – mean, r_s – dynamics; Pool et al., 2018). Values closer to 1 indicate a better match.

Download

4.1 Past model evaluation

The magnitudes and seasonality of flow are generally well represented during the evaluation period (1990–1999) by the calibrated CWatM_base and CWatM_glacier (Figs. 3 and S8). For individual years, the discharge dynamics are also well matched by the simulations (Fig. S5).

The Fraser and Glomma basins have a similar discharge seasonality characterized by a melt regime with peak discharge during the early summer (May–July). Simulated discharge peaks later and declines earlier in the year than observed discharge for both basins, and early winter discharge (October–December) is overestimated. The discharge regimes of the Rhône and Rhine river basins, on the other hand, are characterized by higher discharge during the winter months and relatively low discharge during late summer (August–September) (Fig. 3). The simulations overestimate discharge for the Rhine River basin throughout most of the year, with a mean difference of 8 %. In contrast, the simulations closely match observed discharge for the Rhône River basin. The distribution of flow magnitudes is well represented, but annual maximum discharges are underestimated for the Fraser River basin and overestimated for the Rhine and Rhône river basins. The annual minimum discharge is underestimated for both model setups for all river basins except the Fraser (Figs. S7 and S8).

Observations are better matched by CWatM_glacier for the Rhône River basin during August and to some extent for the Fraser River basin throughout most of the year. No noticeable differences between the model setups are evident for the Glomma and Rhine river basins. For the Rhône River basin, the differences between the two model setups are more pronounced at the upstream station (see Fig. S4). Simulations of CWatM_base show a bias towards too low discharge in August for the Rhône River basin, which suggests that a relevant process during this time of the year is not represented well in this model setup. Using CWatM_glacier, summer discharge is simulated well (see Figs. 3 and S4). This indicates that CWatM_glacier can capture summer discharge better because ice melt is explicitly included.

Figure 4Relative difference in projected mean annual and mean monthly discharge of CWatM_glacier compared to CWatM_base at the downstream gauges for 30 years from 2070–2099 using five GCMs (150 data points per boxplot). Differences are shown per SSP scenario, which translate into global median warming levels of +2.0 and +4.3 °C compared to pre-industrial times. Note that each subplot has a different y-axis scale.

Download

In general, discharge differences between CWatM_base and CWatM_glacier are marginal, which is also reflected by similar performance in the evaluation period. Model calibration can, to some extent, compensate for inadequate process representations within models (Duethmann et al., 2020; Knoben et al., 2020), and consequently, it can also mitigate the limitations arising from the absence of a glacier representation. During calibration of CWatM_base, a parameter set that exhibits a reasonable ability to simulate summer discharge by leveraging other processes within the model is chosen. Parameters influence the partitioning between evapotranspiration and runoff and the timing of runoff. Moreover, adjustments to the precipitation input (and thus the water volume) of the river basin can be made if a snowfall correction factor is applied. For example, the calibrated parameters of CWatM_base reduce evapotranspiration and enhance groundwater contributions to streamflow in the Rhône River basin to compensate for missing runoff from the glaciers in summer (Figs. S2 and S3). Therefore, model setups that do not include all relevant processes can perform fairly well but do not get the right answers for the right reasons (Kirchner, 2006). This presents a challenge when assessing whether an enhanced glacier representation improves model simulations. Additionally, it becomes problematic in a changing environment, a topic discussed in the following section.

4.2 Future projections

The impact of enhancing the glacier representation in CWatM by coupling it with the glacier model OGGM becomes more apparent at the outlets of our study basins when simulating future projections. Future simulated discharge of CWatM_glacier is lower compared to CWatM_base (Fig. 4), with annual discharge differences of 1 % to 4 %. Monthly differences between the model setups are more pronounced and reach 20 % for the Fraser and 35 % for the Rhône but only 3 % for the Rhine and 6 % for the Glomma River basin (Fig. 4). The reason for more pronounced differences in the future compared to the past is that CWatM_glacier captures glacier retreat and therefore changing glacier-sourced runoff during this century, whereas CWatM_base cannot capture this. The calibration of CWatM_base compensates for lacking glacier representation by adjusting parameters to fit the observed discharge in the past. Therefore, it implicitly includes the ice melt during the calibration period in the parameter sets. If these parameter sets are applied to future simulation periods with less glacier-sourced melt, future discharge is likely overestimated in the river basins due to the stationarity of parameters. The differences between the model setups are more pronounced for the higher-emission scenario SSP5–8.5 and in months where ice melt occurs (June to October). However, differences can also be seen during the snowmelt period.

Figure 5Absolute change in mean annual and mean monthly discharge at the downstream stations for the 2070–2099 period compared to 1990–2019 for CWatM_base and CWatM_glacier, shown per SSP scenario, which translates into global median warming levels of +2.0 and +4.3 °C compared to pre-industrial times. The height of each bar indicates the median change based on the GCMs, and the grey lines indicate the range between the maximum and minimum change in individual GCMs.

Download

Comparing future discharge projections to past discharge simulations yields valuable information about the future change in discharge, which is an important metric for adaptation planning. Future change projections of discharge are more negative or less positive for CWatM_glacier compared to CWatM_base (Fig. 5) because CWatM_base slightly underestimates discharge in the past during summer months and likely overestimates future summer discharge, which results in lower future changes.

Differences between the model setups are the largest for August where simulations with CWatM_glacier project a change of −29 % (−58 %) for SSP1–2.6 (SSP5–8.5) for the Rhône River basin towards the end of the century, whereas the projected change with CWatM_base is lower (−2 % (−30 %)). On an annual scale, the differences are lower than for the summer months, with +2 % (−8 %) change in discharge for CWatM_glacier and +6 % (−6 %) changes for CWatM_base (Fig. S12).

For both model setups and all selected river basins, future change in discharge is positive in winter months and negative in summer months, with a larger change for the higher-emission scenario (Fig. 5). These future changes are consistent with previous studies (Stahl et al., 2022; Rottler et al., 2020; Schneider et al., 2013). The range of projected change is large due to different regional projected changes in the forcing data of the GCMs.

The effect of coupling is larger for the summer months, where glaciers contribute the most to discharge. Therefore, it is advisable to include a realistic glacier representation in the hydrological modelling of large river basins, especially when looking beyond the mean annual cycle.

Figure 6Relative mean glacier contribution to annual and monthly discharge (glacier share) at the downstream gauge for the 1990–2019 and 2070–2099 periods for two SSP scenarios, which translate into global median warming levels of +2.0 and +4.3 °C compared to pre-industrial times. The height of the bar indicates the median glacier share of GCMs, and the grey lines indicate the minimum and maximum glacier share of GCMs.

Download

4.3 Glacier-sourced melt contribution

The previous section has revealed the effect of coupling a large-scale hydrological model to a glacier model and its added value in climate change studies for specific large river basins. This section exemplifies the possible use of the coupling to identify glacier-sourced melt contribution to discharge. The coupled model setup was run excluding glacier-sourced melt and the glacier areas, which are updated annually (CWatM_glacier,bare). This simulation was subtracted from CWatM_glacier to derive the glacier contribution to discharge including melt and rain on glaciers. Simulated annual glacier contributions to discharge show a drastic decline even under a moderate warming level of +2.0 °C (SSP1–2.6 at the end of the century). By the end of the century, these contributions are projected to be less than 1.5 % for the selected basins in Europe and North America and thus negligible (Fig. 6). The monthly relative glacier contribution to discharge is the largest in the Rhône River basin due to the substantial percentage of basin area covered by glaciers and low discharge from non-glacierized areas during the months when glacier-sourced melt occurs (see Figs. 3 and S14). The simulated average glacier contribution in August diminishes from 37 % to 14 % (10 %) (340 to 82 m³ s⁻¹ (38 m³ s⁻¹)) at the end of the century for a warming level of +2.0 °C (+4.3 °C) (Figs. 6 and S15). Note that we only show 30-year averaged results here (Fig. 6), but a similar approach could be taken to estimate glacier-sourced melt contributions in individual months and years.

5.1 Past model evaluation

We used six parameter sets within CWatM for the global evaluation to include parameter uncertainty in this study, which is rarely done in large-scale hydrological modelling. The performance of these parameter sets for CWatM_base was comparable to other large-scale hydrological models at 30 arcmin resolution (Müller Schmied et al., 2021; Sutanudjaja et al., 2018). Overall, simulations tend to overestimate observed discharge, which is in line with other studies (Müller Schmied et al., 2021; Wiersma et al., 2022) (see Figs. S9 and S10).

Figure 7(a) Performance comparison using the same discharge stations as presented in Wiersma et al. (2022) between CWatM_base and CWatM_glacier for individual calibrations (grey dots) and the mean of all calibrations (coloured dots) for the 10-year period from 2004 to 2013. The performance metric used is NPE (Pool et al., 2018). The Santa Cruz River basin lies outside the figure boundaries (NPE_glacier is −3.2; NPE_base is −195). The dots with grey outlines show basins smaller than 10 000 km². The median performance increase was 0.05 with an interquartile range from 0 to 0.2. Performance increased for 23 out of 30 basins. (b, c) Comparison of mean discharge of observations and simulations by CWatM_base and CWatM_glacier. The individual lines show the simulations of the different parameter sets, and the bold line shows the mean of the simulations.

Download

Figure 8Relative difference in average discharge for September between CWatM_glacier and CWatM_base during the period of 1990–2019. Positive values indicate a larger discharge of CWatM_glacier. The outlines of large-scale glacierized river basins according to Huss and Hock (2018) (basin area > 5000 km², glacier coverage > 0.01 %) are shown in black. Note the difference in degree intervals.

Looking at the 56 large-scale glacierized river basins used by Huss and Hock (2018), observed discharge data are available for 30 of these basins, mainly located in Europe and North America (Wiersma et al., 2022), with glacier coverage ranging from 0.02 % to 35 %. The performance in simulating discharge improved or is at least similar for CWatM_glacier compared to CWatM_base for all these river basins (Fig. 7a), except for the Skagit River basin. This river basin is comparatively small (7850 km²) and therefore likely not well represented at 30 arcmin resolution, where it covers only five grid cells.

The improvement in performance when explicitly including glaciers, as compared to CWatM_base, was the largest for highly glacierized river basins (Fig. 7a). For example, the Copper River basin, which is 22 % glacierized, exhibited substantial improvement in simulation performance when glaciers were explicitly included. In this case, the simulation of CWatM_glacier matched the observations very well. The analysis of the discharge regime revealed that for CWatM_base, ice melt was missing as the relevant process in the Copper River basin during the summer months. Consequently, simulations with CWatM_base failed to capture the observed behaviour (Fig. 7b).

In general, performance was very variable among the river basins for both model setups and all parameter sets. This is likely due to using global parameter sets which were not adjusted for each basin individually and inaccurate precipitation input. For example, the glacier inclusion improved the simulations of the Fraser River basin, but simulations were still unsatisfactory because discharge was highly underestimated (Fig. 7c). Thus, further effort in calibration is likely to enhance the advantages of CWatM_glacier for global simulations. However, improving the global performance of large-scale hydrological models is a research field in itself and beyond the scope of this study.

Figure 9Comparison of relative future discharge change for the month with the largest absolute glacier-sourced melt contribution in the past, at the end of the 21st century, considering both CWatM_base and CWatM_glacier, across 56 glacierized river basins for SSP1–2.6. (a) The coloured dots show the median of all GCMs, and the grey dots show individual GCMs. (b) The boxplots show the relative future change in all basins for CWatM_base and CWatM_glacier and the difference between the two.

Download

5.2 Effect of coupling (past and future differences)

Differences between past simulations of CWatM_glacier and CWatM_base can be seen worldwide in large glacierized river basins (Fig. 8). The relative difference in monthly discharge between the two model setups can exceed 50 % close to the glaciers and decreases further downstream. Nevertheless, differences between the model setups can also be detected downstream. For example, in September, the CWatM_glacier simulations show 15 % higher discharge in the Rhône Basin, 26 % higher discharge in the Glomma Basin and 17 % higher discharge in the Amu Darya Basin. This shows the influence of glaciers on discharge in large river basins. The timing of the largest differences is delayed further downstream due to travel time through the river network. In winter, no glacier-sourced melt is expected; thus the change in model setup should not influence the winter discharge. However, a decrease in discharge is simulated for winter and spring months in CWatM_glacier compared to CWatM_base (see Fig. S21). The reason behind this is likely the lack of glacier representation in CWatM_base where other runoff processes occur on the glacier area. Another explanation is the finer spatial resolution of elevations in OGGM (see Table 1) and therefore a later snowmelt occurrence in OGGM than in the same areas in CWatM.

On an annual average, the discharge simulated by CWatM_glacier is larger than or equal to that of CWatM_base. This is because the glaciers are diminishing, resulting in a negative change in storage that contributes to discharge. This additional water source simulated by OGGM is included in CWatM_glacier but not in CWatM_base. The increase in annual discharge is larger for highly glacierized basins (Fig. S19). The annual difference between the model setups exceeds 10 % for 15 basins and 5 % for 22 basins.

Regarding future changes in river discharge, changes can be positive or negative depending on the regional changes in temperature and precipitation. In many river basins, the discharge in the month with the largest past glacier-sourced melt contributions is projected to decrease at the end of the century compared to the past (Fig. 9), whereas the annual discharge is projected to increase in many basins (Fig. S22). At the outlets of the 56 glacierized river basins, negative changes are exacerbated when glaciers are included explicitly in large-scale hydrological modelling, whereas positive changes are attenuated (Fig. 9). The reason is that simulated glacier-sourced melt is smaller at the end of the century than in the past, and this negative change is included in CWatM_glacier, whereas precipitation projections are identical for both model setups. This pattern holds for the majority of river basins, with the exception of a few basins where simulated glacier-sourced runoff is greater or similar at the end of the century compared to the past because glacier peak water is reached later, most notably Santa Cruz, Tarim, Amu Darya and Jokulsa (Fig. S17). The difference in the change between the model setups is more pronounced for basins with a higher degree of glaciation. Overall, changes are more substantial for the higher-emission scenario (Fig. S23).

Figure 10(a) Simulated glacier contribution to total basin discharge in the past (1990–2019) and at the end of the century (2070–2099) for the month with the highest relative glacier-sourced melt contribution to discharge in the past. (b) The difference between the two periods, i.e. the simulated future change in glacier contribution. The outlined circles in (a) show the Amu Darya and Tarim basins (black and grey, respectively).

Download

5.3 Glacier-sourced melt contribution

The simulated relative glacier-sourced melt contributions to discharge decreased for future periods at the outlets of all glacierized river basins, except for the Amu Darya and Tarim basins, both on an annual average and during the month with the largest glacier-sourced melt contribution (Fig. 10). In the Amu Darya and Tarim basins, the projected future glacier-sourced melt is either larger than or similar to the past. This, combined with projected decreased discharge, increases relative glacier-sourced melt contributions. The decrease in glacier-sourced melt contribution to discharge is large in many basins. The highest monthly glacier contribution decreased by more than 10 % in 22 river basins for SSP1–2.6 and 26 river basins for SSP5–8.5. It decreased by more than 20 % in 12 river basins for SSP1–2.6 and 10 river basins for SSP5–8.5.

However, it should be noted that deriving glacier-sourced melt contributions can be highly uncertain if observed discharge is not accurately represented in discharge simulations. An underestimation of discharge would result in an overestimation of relative glacier-sourced melt contributions, whereas an overestimation of discharge would result in an underestimation of relative glacier-sourced melt contributions.

6.1 Benefits of coupling

Glaciers were explicitly included in the simulations of the large-scale hydrological model CWatM, which simulates human water use in addition to natural hydrological processes. Incorporating human water use is essential for studying large river basins. No less important is improving the representation of mountainous areas within large-scale hydrological models.

Therefore, including glacier-sourced runoff in large-scale hydrological modelling is an important step towards bridging the gap between assessments of natural environmental changes and their potential impacts on societies on larger scales. It is also a step towards better understanding the relevance of glaciers for water availability beyond small headwater catchments because the coupling opens the possibility of estimating glacier contributions to discharge by running the model in different setups. Thus, the inclusion of glacier-sourced runoff offers the opportunity for further studies to delve into the impacts of glaciers on water availability.

To the best of our knowledge, our study is the first to explicitly incorporate glacier-sourced runoff generated on a daily resolution and account for ice dynamics in a large-scale hydrological model (see Cáceres et al., 2020; Wiersma et al., 2022). Our analysis demonstrated that explicitly including glaciers in both past and future simulations leads to notable changes in simulated discharge. Simulation differences were mainly present in the summer months, a period where water availability is especially relevant for irrigated agriculture. In specific years or months, the effect of including glaciers is likely larger than shown in this analysis, for which 30-year average values were used. This is potentially relevant for climate impact assessments with a focus on summer months or extreme years.

Climate change leads to an increase and subsequent decline in glacier-sourced melt volumes, occurring once when the point of glacier peak water is surpassed (Huss and Hock, 2018). Therefore, the impact of improved glacier representation on future discharge changes depends on the period of comparison relative to glacier peak water. If the past period coincides with glacier peak water, the inclusion of glaciers in CWatM simulations results in a better representation of negative future discharge changes. This pattern holds for many river basins worldwide in this study when comparing 1990–2019 to 2070–2099, including the four selected river basins (see Figs. 9 and S17). The rapid changes in glacier-sourced melt volumes throughout this century and the dependence of results on the selected period emphasize the relevance of explicitly including glacier representation in large-scale hydrological modelling studies that focus on future discharge changes.

6.2 Caveats of coupling

The results of coupling a glacier model (OGGM) to a large-scale hydrological model (CWatM) are influenced by the limitations of each model. A major limitation in both models is uncertain precipitation data. Global glacier models such as OGGM are further limited by the scarcity of long-term in situ mass-balance observations (WGMS, 2020) and thus mostly rely on one observation per glacier (i.e. the 2000–2020 geodetic mass-balance average estimate of Hugonnet et al., 2021). Nevertheless, we recommend including glacier runoff simulations for future discharge projections in large river basins with many glaciers because glaciers are a rapidly changing water resource and influence future discharge changes, especially in the summer. The absence of globally available and temporally and spatially resolved mass-balance estimates poses a major challenge for the robust calibration of more enhanced mass-balance models (such as energy-balance models). Yet, such mass-balance models would be valuable for improving the representation of tropical glaciers, particularly since their current representation is limited by temperature-index models (Fernández and Mark, 2016). Moreover, model parameter choices affect glacier volume and runoff projections (Rounce et al., 2020; Schuster et al., 2023b).

Large-scale hydrological models, including CWatM, are limited by the coarse representation of mountainous areas when typical resolutions of 10–50 km (5–30 arcmin) are used. We did not conduct a comparison between the simulations at different modelling scales for the selected river basins because simulation differences could not solely be attributed to scale differences but were also a result of different calibration procedures. At the global scale, the calibration and evaluation of models are limited by the availability of discharge data, which are biased towards mid-latitudes and temperate climates (Krabbenhoft et al., 2022). Additionally, large-scale hydrological models often lack a detailed calibration, resulting in poor model performance, especially in arid regions and regions with few discharge observations (Burek et al., 2020). Regionalized and calibrated parameter sets that encompass all major hydrological processes would certainly improve simulations but are beyond the scope of this study. In the following section, the most important limitations of the models and their coupling are discussed in more detail.

6.3 Comparison to other studies

The performance of global hydrological simulations increases in particular in highly glacierized river basins when glaciers are explicitly included in modelling (Wiersma et al., 2022). While Wiersma et al. (2022) demonstrated a decrease in model performance in specific months for weakly glacierized basins when glaciers are explicitly included, our study did not find an overall decrease in performance. This might be explained by the underestimation of discharge in many basins in the present study, whereas discharge was mostly overestimated in Wiersma et al. (2022). Like Wiersma et al. (2022), we detected a decrease in winter and spring discharge for the coupled model. Comparing glacier-sourced melt contributions to other studies is challenging due to differences in modelling time periods and varying definitions of what constitutes glacier melt (Gascoin, 2024). Nevertheless, we compare our results to other studies to see whether trends and magnitudes match.

The derived glacier-sourced contribution to discharge in the Rhine River basin at 5 arcmin resolution fits the results of a detailed study on discharge components of the Rhine River basin (Stahl et al., 2016), e.g. 2.6 % (4.2 %) ice melt at Lobith in August (September) (1901–2006) vs. 7.5 % (6 %) glacier-sourced melt in this study (1990–2019). Differences between the studies can be explained by different periods and by comparing ice melt to glacier-sourced melt. While our study suggests a glacier-sourced melt contribution of 40 % for the upper Indus, previous studies reveal a contribution of 60 % (Lutz et al., 2014) and 45 % (Su et al., 2016). Nevertheless, all studies show a similar pattern of high glacier contributions for the upper Indus and rather low glacier contributions for the upper Brahmaputra. The glacier contributions we derived for the Himalayan basins are greater than the results reported by Kraaijenbrink et al. (2021). Differences might stem from differences in evaporation estimates or other discharge components. While Huss and Hock (2018) show a strong decrease in relative glacier-sourced runoff for basins in Central Asia, our results show an increase in glacier contribution to discharge in the Amu Darya and Tarim river basins. This projected increase in glacier contributions is consistent with only a slight decrease in projected glacier volumes in these basins compared to all other regions in High Mountain Asia at the end of the 21st century, as reported by Miles et al. (2021).