Journal cover Journal topic
Hydrology and Earth System Sciences An interactive open-access journal of the European Geosciences Union
Journal topic
Hydrol. Earth Syst. Sci., 23, 1833-1865, 2019
https://doi.org/10.5194/hess-23-1833-2019
Hydrol. Earth Syst. Sci., 23, 1833-1865, 2019
https://doi.org/10.5194/hess-23-1833-2019

Research article 03 Apr 2019

Research article | 03 Apr 2019

# Future evolution and uncertainty of river flow regime change in a deglaciating river basin

Future evolution and uncertainty of river flow regime change
Jonathan D. Mackay1,2, Nicholas E. Barrand1, David M. Hannah1, Stefan Krause1, Christopher R. Jackson2, Jez Everest3, Guðfinna Aðalgeirsdóttir4, and Andrew R. Black5 Jonathan D. Mackay et al.
• 1School of Geography, Earth and Environmental Sciences, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK
• 2British Geological Survey, Environmental Science Centre, Keyworth, Nottingham, NG12 5GG, UK
• 3British Geological Survey, Lyell Centre, Research Avenue South, Edinburgh, EH14 4AS, UK
• 4Institute of Earth Sciences, University of Iceland, 101 Reykjavík, Iceland
• 5Geography and Environmental Science, University of Dundee, Dundee, DD1 4HN, UK
Abstract

The flow regimes of glacier-fed rivers are sensitive to climate change due to strong climate–cryosphere–hydrosphere interactions. Previous modelling studies have projected changes in annual and seasonal flow magnitude but neglect other changes in river flow regime that also have socio-economic and environmental impacts. This study employs a signature-based analysis of climate change impacts on the river flow regime for the deglaciating Virkisá river basin in southern Iceland. Twenty-five metrics (signatures) are derived from 21st century projections of river flow time series to evaluate changes in different characteristics (magnitude, timing and variability) of river flow regime over sub-daily to decadal timescales. The projections are produced by a model chain that links numerical models of climate and glacio-hydrology. Five components of the model chain are perturbed to represent their uncertainty including the emission scenario, numerical climate model, downscaling procedure, snow/ice melt model and runoff-routing model. The results show that the magnitude, timing and variability of glacier-fed river flows over a range of timescales will change in response to climate change. For most signatures there is high confidence in the direction of change, but the magnitude is uncertain. A decomposition of the projection uncertainties using analysis of variance (ANOVA) shows that all five perturbed model chain components contribute to projection uncertainty, but their relative contributions vary across the signatures of river flow. For example, the numerical climate model is the dominant source of uncertainty for projections of high-magnitude, quick-release flows, while the runoff-routing model is most important for signatures related to low-magnitude, slow-release flows. The emission scenario dominates mean monthly flow projection uncertainty, but during the transition from the cold to melt season (April and May) the snow/ice melt model contributes up to 23 % of projection uncertainty. Signature-based decompositions of projection uncertainty can be used to better design impact studies to provide more robust projections.

1 Introduction

Mountain watersheds have been referred to as the world's water towers , partly because they receive large quantities of precipitation relative to adjacent lowlands, but also because they regulate runoff through the seasonal accumulation and melt of snow and ice. The presence of snow and ice profoundly affects characteristics of downstream river flow regime including flow magnitude, timing and variability over a range of timescales . This is partly due to the periodic (diurnal and seasonal) variations and longer-term (decadal) trends in meltwater inputs brought about by fluctuations in glaciological mass balance. In addition, the dynamic water storage and release properties of snow and ice (runoff-routing) control downstream river flow response to runoff over hourly to seasonal timescales (Willis2005). As such, glaciated basins exhibit river flow regimes that differ from their non-glaciated equivalents. For example, the so-called “compensation effect” has been widely observed in the Northern Hemisphere, whereby partially glaciated catchments demonstrate reduced intra-annual flow variability . Indeed, the compensation of runoff from melt inputs can actually serve to increase mean runoff during anomalously dry heatwave events .

Mountain glaciers are retreating at unprecedented rates while snow coverage is receding , resulting in observable changes to downstream river flows . With near-surface air temperature projected to rise over the coming decades future changes in river flow regimes in response to cryosphere change could have wide-ranging socio-economic and environmental impacts. Long-term reductions in meltwater inputs will disrupt the supply of water available for irrigation . Increased inter-annual and intra-annual flow variability will threaten infrastructure projects such as hydroelectric power stations . The loss of the runoff-regulating effects of snow and ice could result in more frequent short-term very high flows putting downstream populations and infrastructure at risk . Changes in flow magnitude and variability from annual to sub-daily timescales will threaten the sustainability of some of the world's most pristine freshwater ecosystems . Therefore, it is of paramount importance to make reliable projections of changes in downstream river flow regimes from glaciated watersheds so that future impacts can be adapted to and mitigated.

Computational glacio-hydrological models (GHMs) driven by numerical climate model projections allow us to assess how future river flow regimes will change in glaciated river basins. Past studies have focussed on projecting changes in decadal, annual and seasonal variations in runoff magnitude. Decadal changes in runoff are inevitable over the coming century , where enhanced melt will result in increased river discharge to a point in time termed “peak water”, after which the continued loss of snow and ice will result in an overall decrease in river flow. It has been shown that many basins, particularly those with small glaciers, have already reached peak water and face a future of dwindling water supply . Seasonal flow magnitudes are also projected to change as melt cycles evolve and watersheds shift from glacial–nival to pluvial runoff regimes .

Some impact studies show robust changes in the magnitude of the highest and lowest river flows, including , who projected an increase in the magnitude of the 10 % exceedance flow (Q10) for river basins across the Hindu Kush Himalayan region. Other studies for the Rhine , upper Indus and upper Yellow River basin show high-flow magnitudes will increase. projected a decrease in low-flow magnitude (Q90) for the snow-covered Sierra Nevada and upper Colorado River basins due to shifts in the snowmelt season and changes in precipitation type from snow to rain. For the Hindu Kush, found the opposite impact with an increase in the magnitude of low-flow events. The projected trends in Q90 for the upper Yellow River basin by were inconclusive as they showed an even spread of positive and negative trends under the warmest climate scenarios.

Of course, one could go beyond projecting changes in seasonal to decadal mean flow magnitudes and quantiles of the flow duration curve (FDC). A branch of streamflow analysis that has been widely adopted in hydrology is the calculation of river flow “signatures”, which are metrics derived from river discharge time series that represent different characteristics of river flow over specific timescales. These may include mean flows and FDC quantiles as well as metrics to quantify the variability (e.g. coefficient of variation), timing (e.g. peak flow month) and flashiness (e.g. autocorrelation) of flows. Signatures have been used in the past to analyse catchment runoff behaviour and similarity . Furthermore, their ability to localise specific aspects of runoff behaviour make them ideal diagnostic evaluation metrics for model hypothesis testing and calibration . They also offer an opportunity to evaluate past and future river flow regime change. For example, projected changes in 14 different river flow signatures for 14 snow-covered catchments in Sweden and showed daily to annual river flow magnitude, timing and variability were all sensitive to climate change. An analysis like this is yet to be undertaken for any glaciated river basins.

Projections of river flow regime are inherently uncertain due to assumptions made about the formulation, parameterisation, and boundary conditions of the underlying GHM and climate model, be that a general circulation model (GCM) or a combined GCM and regional climate model (GCM–RCM) . Uncertainties may also be introduced by intermediary steps employed to link the two sets of models together, such as downscaling (DS). Quantifying the propagation of uncertainties from all sources in the model chain provides a basis for assigning more robust levels of confidence to river flow projections. Additionally, one can assess the relative contributions of model chain components to the total projection uncertainty, providing empirical evidence for future research needs (e.g. Meresa and Romanowicz2017). Ensemble-based experiments have been used in the past to provide this understanding. Here, different components of the model chain are perturbed, typically using a one-at-a-time (OAT) approach where the spread in projections for each perturbed component is evaluated. perturbed three components of a model chain applied to the Hunza River basin, northern Pakistan, including the GCM, statistical DS model and parameterisation of the GHM. They showed that all three sources contributed to annual runoff projection uncertainty, but for the heavily glaciated subregions of the catchment, the GHM parameter uncertainty exceeded the effect of other sources. investigated uncertainty in seasonal river flow projections over the 21st century for the Findelengletscher catchment, Switzerland, by modifying the GCM–RCM, GHM melt model structure and initial ice volume boundary condition. Of these, they found that the GCM–RCM and initial ice volume were most important, while the melt model structure was of secondary importance. investigated uncertainties in 21st century river flow projections for the Clutha River basin, New Zealand. They evaluated contributions from the emission scenario (ES), GCM–RCM, statistical DS approach and melt model structure. Similarly to , they found that uncertainty in the choice of GCM–RCM dominated total projection uncertainty.

The OAT method provides a first-order approximation of the relative contribution of each component to the total projection uncertainty. However, findings are dependent on how the non-perturbed model components are fixed. Furthermore, this approach cannot resolve interactions between model components which may also contribute to projection uncertainty . The analysis of variance (ANOVA) statistical method addresses these shortcomings and has been adopted in a number of recent regional- and global-scale hydrological modelling studies to compare uncertainties stemming from the ES, climate model, hydrological model structure and DS approach. While uncertainties associated with future climate tend to dominate projections of river flow, glacier-fed river flow projections have shown been to be highly sensitive to hydrological model structure , particularly in relation to high flows . Furthermore, the contribution of projection uncertainty from interactions between model chain components can exceed individual components . Several issues not considered in these studies, however, are yet to be addressed. Firstly, none have investigated a full range of characteristic changes in river flow regime covering decadal to sub-daily timescales. Second, all have incorporated hydrological model uncertainty using multiple model codes, each with their own unique set of process representations, resolution, time step and climate interpolation strategies, making it difficult to determine which model components contribute most to projection uncertainty. Finally, none included a fully integrated mass-conserving, dynamic glacier evolution model component and therefore could not fully account for atmosphere–cryosphere–hydrosphere feedbacks.

This study uses a GCM–RCM–DS–GHM model chain to simulate the impact of 21st century climate change on downstream river flow regime in the deglaciating Virkisá river basin in southern Iceland. Five components of the model chain are perturbed to represent uncertainty of ES, GCM–RCM, statistical DS parameterisation and structure parameterisation of two primary controls on river flow regime in the GHM: melt and runoff-routing processes. The study has two principal aims: (i) to determine how climate change and consequent cryospheric change will impact downstream river flow regime over the 21st century; and (ii) to quantify the relative influence of the five model chain components to projection uncertainty across the different characteristics of river flow regime. This study addresses each of the aforementioned gaps in previous work. Firstly, changes in river flow regime are assessed quantitatively using 25 river discharge signatures which define different characteristics of river flow regime over a range of timescales. Second, a single, consistent, GHM code is used that can incorporate different model structures and parameterisations of melt and runoff-routing processes, allowing for uncertainty stemming from these to be localised using ANOVA. Finally, a fully integrated mass-conserving, dynamic glacier evolution routine is included in the GHM code.

2 Methodology

## 2.1 Study site

The Virkisá river basin covers an area of 22 km2 on the western side of the ice-capped Öræfajökull stratovolcano in south-east Iceland (Fig. 1) and forms a primary drainage channel for accumulating ice at the mountain summit (∼2000 m a.s.l.). The glacier flows in a south-westerly direction (average ice surface slope =0.25) along two distinct glacier arms, Virkisjökull and Falljökull (hereafter referred to as Virkisjökull), around a bedrock ridge before meeting again at the terminus (∼150 m a.s.l.). Virkisjökull currently covers ∼60 % of the river basin area but has been in a phase of retreat since 1990. Between 1988 and 2011 Virkisjökull lost ∼0.3 km3 of ice and retreated 0.5 km. A small proglacial lake at the terminus forms the headwater of the Virkisá river. The Virkisá flows through an 800 m bedrock-controlled section flanked on either side by push moraines and then over the Skeiðarársandur floodplain, typically comprising unconsolidated glacial outwash sediment. The steep-sided valley walls and glacial activity only allow for sporadic development of thin soils with limited vegetation including mosses, grass and shrubs.

Figure 1Location of the Virkisá river basin in Iceland with glaciated areas highlighted in grey (a), location of Öræfajökull (b), and detailed topographical map of the basin including major land surface types as well as meteorological and stream gauging stations (c).

The local climate is characterised by cool summers (∼10C on average at automatic weather station 1, AWS1) and mild winters (∼1C on average at AWS1) with an average temperature lapse rate of −0.44C 100 m−1 (Flett2016). There is a significant lateral precipitation gradient due to the prevailing north-easterly winds and orographic effects with more than 5 times the annual precipitation falling at the Öræfajökull summit (∼8000 mm yr−1) compared to lower down at the catchment outlet to the west (∼1500 mm yr−1) .

## 2.2 Climate data

### 2.2.1 Historical climate

Historical climate data were available from 1981 to 2016 inclusive. A detailed description of these are provided by . For brevity, only a summary of these data is provided here. The historical climate data include continuous hourly near-surface air temperature measurements from two automatic weather stations in the catchment (Fig.  1c) which were installed in 2009 (AWS1) and 2011 (AWS4). Temperature data from the nearby Icelandic Meteorological Office Fagurhólsmýri weather station (12 km south of the study site) were used to extend the AWS1 time series back to 1981 using a linear regression model (R2=0.92) to bias-correct against the AWS data. A seasonally variable hourly lapse rate calculated between AWS1 (156 m a.s.l.) and AWS4 (805 m a.s.l.) was used to extrapolate near-surface air temperature across the study region, and an on-ice temperature correction function was employed to account for katabatic cooling of air in the glacier valleys. Continuous hourly incident solar radiation data were also available from AWS1. A random resampling strategy that accounted for the dependence between intraday solar radiation and temperature variability was employed to generate a continuous time series back to 1981. For precipitation, the 2.5 km gridded hourly total precipitation data produced as part of the ICRA atmospheric reanalysis project were used, which are currently considered the most accurate gridded precipitation product over Iceland . These were bias-corrected against hourly precipitation measurements from AWS1 using equidistant quantile mapping .

### 2.2.2 Regional climate projections

Future climate time series until 2100 were constructed using regional climate projections from the Coordinated Regional Climate Downscaling Experiment (CORDEX) . These are based on an ensemble of RCMs driven by GCM projections from the Coupled Model Intercomparison Project (CMIP5) . Iceland is covered by the EURO-CORDEX and ARCTIC-CORDEX regional model domains. Following the review by , the EURO-CORDEX data were used as these include projections at a higher 0.11 spatial resolution and a larger ensemble of GCM–RCM combinations allowing better exploration of climate model uncertainty.

The 0.11 EURO-CORDEX simulations span the years 1950–2100 with simulations up to 2005 constituting the “recent past” where influences such as atmospheric composition, solar forcing and emissions are imposed based on observations. From 2006, three future ESs or representative concentration pathways (RCPs) were imposed on the models, including RCP2.6, RCP4.5 and RCP8.5, which represent an additional radiative forcing by 2100 relative to pre-industrial values of +2.6, +4.5 and +8.5 W m−2 respectively. All simulations are available at 3-hourly to 3-monthly resolution; however, the 3-hourly simulations were only produced using four GCM–RCM combinations while daily to seasonal simulations were produced using 15. Given the intent of this study to analyse projection uncertainty, it was decided that the daily data were most suitable. RCP2.6 was omitted as only 8 of 15 GCM–RCM combinations within the CORDEX archive used this ES. Furthermore, the probability of achieving the RCP2.6 targets is increasingly unlikely and arguably completely infeasible given the current global emission trajectory.

One of the 15 GCM–RCM combinations (GCM: CNRM-CM5, RCM: CNRM-ALADIN53) was removed from the ensemble given that it showed an extreme negative winter temperature bias and a consistently low skill when compared to daily observed climate data (see Appendix Appendix A). Figure 2 shows the seasonal bias of each of the 14 remaining GCM–RCM combinations when compared to observations between 1981 and 2005. For temperature, the coldest days (T1) typically show a negative bias, particularly in winter, spring and autumn. Biases for T99 are generally positive but smaller in magnitude. The average absolute bias in mean seasonal temperature (Tmean) is 1.4 C, but the majority of GCM–RCM combinations show absolute biases <1.2C. Biases in seasonal incident solar radiation projections are almost exclusively positive, with the largest biases associated with SWmean and SW99, particularly in spring and summer where they can exceed 100 W m−2. Total precipitation biases are typically largest in winter and autumn where proportionally biases in Pmean can exceed the magnitude of the observations (see SON for EC-EARTH–HIRHAM5). The largest biases however are seen in extremes (P99) which range from −86.9 to 77.5 mm d−1. While positive and negative precipitation biases are present throughout the ensemble, the sensitivity of precipitation simulations to the RCM is clear. For example, the CCLM4-8-17 RCM has a systematic negative bias and the HIRHAM5 RCM has a systematic positive bias.

Figure 2Comparison of seasonal catchment-average observed and simulated near-surface air temperature (T), incident solar radiation (SW) and total precipitation (P) between 1981 and 2005 for the 14 GCM–RCM models used in this study. The top row shows the observed value and all subsequent rows indicate the GCM–RCM biases. The 1st percentile, mean and 99th percentile are denoted by the subscripts 1, mean and 99 respectively. All statistics are calculated for the recent past (1981–2005) for winter (DJF), spring (MAM), summer (JJA) and autumn (SON).

### 2.2.3 Downscaling regional climate projections

The statistical delta-change downscaling approach which has been widely applied in hydrological impact studies was employed. While most studies have used monthly mean delta-change values to capture seasonal shifts in climate, several recent investigations have used advanced quantile-based approaches which account for changes in higher-order statistical properties of future climate by evaluating shifts in the empirical cumulative distribution functions (ECDFs) of climate variables. Including these higher-order changes has shown to be important for evaluating shifts in extreme high flows and sub-seasonal metrics of river flow projections . In addition, shifts in the day-to-day variability of temperature impact projections of glacier retreat as these variations control the periodic rising of temperature above the melting point . Accordingly, the advanced delta-change approach was adopted in this study. The approach is summarised in five steps which were applied to each combination of GCM–RCM, climate variable and ES separately:

1. The climate variable time series was divided into four 25-year-long periods including the recent past (1981–2005) and early (2006–2030), mid- (2041–2065) and late (2076–2100) 21st century.

2. For each of the four periods, all daily data points were further divided into 12 subsamples representing each month of the year. An ECDF was constructed for each month of each period.

3. For each month of each future period, 10 deltas were calculated by taking the mean difference between the recent past and future ECDF for each 10 % section (see grey bars in Fig. 3a for example).

4. Given the need for transient climate time series to simulate glacier evolution over the 21st century, a daily delta time series from 2006 to 2100 was constructed for each ECDF section of each month by linearly interpolating between the calculated deltas of each future period (e.g. as implemented by  Farinotti et al.2012), using the midpoints of the future periods as interpolation points (Fig. 3b).

5. The hourly historic observation data for the recent past were randomly sampled (with replacement) on a year-by-year basis to generate an initial unperturbed future climate variable time series (blue dash, Fig. 3c). The daily deltas were applied to this time series for each month and ECDF section separately to generate a future perturbed climate time series at an hourly resolution (orange dash Fig. 3c). It was noted upon visual inspection that the inter-annual variability of the future climate time series was very sensitive to the random sampling of the historic climate data. Accordingly, uncertainty associated with this aspect of the DS parameterisation was considered by using 10 different random historic climate samples.

Figure 3Example of the advanced delta-change approach when applied to near-surface air temperature data based on the RCP8.5 projections using the CNRM-CM5 GCM and CCLM4-8-17 RCM. Deltas (grey bars) derived from ECDFs (black curves) for April in the late 21st century (a); daily delta time series for each section of the April ECDFs (green line represents 40th–50th percentile section) (b); initial and perturbed future temperature time series when deltas for all months and ECDF sections are applied (c).

For temperature, catchment-average daily deltas were applied evenly across the catchment and each daily period of the unperturbed time series. Accordingly, diurnal temperature variability and lapse rates were assumed not to change in the future. For incident solar radiation and precipitation, proportional deltas were used to prevent negative values and preserve the sub-daily proportional distribution of these variables in space and time. A total of 280 future climate time series from 2 ESs with 14 GCM–RCM models and 10 DS parameterisations were generated for this study.

## 2.3 Glacio-hydrological model

The distributed GHM code implemented by was used in this study because it includes a dynamic, mass-conserving glacier evolution component and also allows the user to utilise different model structures of melt and runoff-routing processes. The GHM resolves glacio-hydrological processes over a regular 2-D Cartesian grid of 50 m cells driven by hourly climate data including precipitation, temperature and incident solar radiation. Empirical index-based equations simulate the melt of snow and ice. Snow redistribution by drift and avalanches is calculated using the curvature and slope of the surface while the mass-conserving Δh parameterisation of glacier retreat resolves changes in the glacier geometry. A soil infiltration and evapotranspiration model based on the well-established FAO56 soil moisture accounting procedure solves the water balance for model nodes with no ice or snow coverage. This model has been applied extensively and has been shown to compare favourably to physically based models at the field scale where interception losses are small . Excess soil moisture, rainfall and melt are routed to the catchment outlet via a semi-distributed network of linear reservoir cascades (Ponce1989) which represent the average water storage and release characteristics of the major hydrological pathways in the watershed (see Fig. 2 in Mackay et al.2018).

### 2.3.1 Modification to Δh parameterisation of glacier retreat

Under periods of sustained positive mass balance, simulations from the Δh glacier evolution model may result in an unrealistic build-up of ice at the glacier tongue without any simulated areal advance. Given the potential for periods of glacier advance under a changing climate, such behaviour is likely to result in significant projection biases. Recently, presented an implementation of the original Δh parameterisation that provides more realistic simulations of glacier advance. They propose running the Δh parameterisation a priori outside of the GHM. A small negative mass balance is used to force the Δh model from an initial glacier profile (ideally its maximum observed extent) until the glacier has disappeared completely. At each step, the glacier mass and geometry are stored in the form of a lookup table. On running the GHM, the retreat/advance of the glacier is derived from the lookup table as a function of the simulated glacier mass. One important drawback of using this static lookup table is that it modifies the behaviour of the Δh formulation during periods of retreat. More specifically, this approach neglects the transient annual sequencing of glacier mass balance which influences simulated glacier geometry due to the non-linear structure of the Δh polynomial that defines the relationship between mass balance and glacier geometry. Accordingly, a modified implementation of the approach was used in this study, which behaves like the original Δh formulation during periods of glacier retreat and allows for the simulation of glacier advance while accounting for mass balance sequencing effects on the model behaviour. For periods of negative glacier mass balance the original Δh formulation was used. The GHM was then modified so that for each simulation year, the simulated glacier mass and geometry were stored in memory. If a positive glacier mass balance (ΔM) was simulated, the GHM would log the current glacier mass (Mcurrent) and then look for the most recent historical simulated glacier mass (Mhist) that exceeded McurrentM. The Δh model was then run with a negative mass balance (ΔM*) so that ${M}_{\mathrm{hist}}+\mathrm{\Delta }{M}^{*}={M}_{\mathrm{current}}+\mathrm{\Delta }M$.

### 2.3.2 Melt and runoff-routing model structures

The selection of melt and runoff-routing model structures was based on the findings of . They applied nine different combinations of three melt model structures and three runoff-routing model structures of varied complexity in the GHM and evaluated their ability to capture a range of river discharge signatures derived from observation time series at automatic stream gauge 1 (ASG1 in Fig. 1c), which has been operational since 2012. They also used observation data of ice melt and snow coverage to derive signatures that described different aspects of these data (see Appendix Appendix B). They showed that while introducing model complexity did improve simulations when evaluated against specific signatures, it did not necessarily result in better consistency across all signatures, emphasising model selection uncertainty. The most complex runoff-routing structure, however, was consistently the least efficient when compared to the two simpler alternatives, particularly in relation to capturing signatures representing high-river-flow events. As such, this model structure was discarded, and only the remaining six combinations of melt and runoff-routing model structures were used in this study. These included three melt model structures: (i) the classic temperature-index model (TIM1) where melt increases linearly with near-surface air temperature above a critical threshold (e.g. Braithwaite1995); (ii) the enhanced temperature-index model (TIM2) proposed by Hock (1999) which accounts for topographic effects on incident solar radiation including surface slope, aspect and shading from the surrounding landform; and (iii) the enhanced temperature-index model (TIM3) proposed by which accounts for topographic effects and also includes a dynamic snow albedo parameterisation , which accounts for the drop in snow albedo as it ages. Each melt model structure was combined with the two runoff-routing structures: (i) a single linear reservoir cascade (ROR1) which routes runoff from all sources (ice melt, snowmelt, rainfall and excess soil water) simultaneously; and (ii) two linear reservoir cascades in parallel (ROR2), where the first represents the slow percolation of water through the snow and firn and the second represents faster flow of water through and over bare ice and over land. The simplest ROR1 structure assumes all catchment water stores delay and diffuse downstream river response to runoff in the same way, effectively fixing the runoff-routing behaviour of the catchment over time. The more complex ROR2 structure accounts for temporal variations in the drainage efficiency of the catchment according to changes in snow and ice coverage.

## 2.4 Signatures of river flow regime

Table 1 lists the 25 signatures of river discharge used to evaluate future changes in river flow regime. The majority of the signatures were selected from past studies and were chosen to reflect the types of changes that one might expect to see in snow- and ice-covered catchments. They also broadly follow those used in the model assessment study of . The signatures are grouped into seven different attributes and further categorised by the characteristic(s) of flow regime that they evaluate and their temporal scale. At the decadal timescale, two signatures were selected. These include the “peak water”, which defines the timing (by year) of maximum flow, as well as the inter-annual flow range, which characterises long-term flow variability. Changes in mean annual river flow were also evaluated, while mean monthly flows were used to evaluate changes to the seasonal timing and magnitude of river flow. The range in mean monthly flows was also chosen to evaluate intra-annual flow variability. In addition, eight signatures were selected which broadly describe the magnitude and variability of slow-release low flows (99 %–95 % exceedance flows), moderate flows (52 %–48 % exceedance) and quick-release high flows (5 %–1 % exceedance). For these, the quantiles of the FDC were used to assess changes in the magnitude of these flow types. The standard deviation was also used to define flow variability of each flow type. Finally, the integral scale, which measures the lag time at which the autocorrelation function of the river flow time series falls below $\frac{\mathrm{1}}{e}$, was utilised as an indicator of the response time of the catchment to runoff events (flashiness).

Table 1Summary of 25 river discharge signatures used to evaluate future changes in river flow regime. Those with available limits of acceptability were also used as part of the GHM calibration and evaluation procedure.

## 2.5 GHM calibration

Given the focus on projecting changes in river discharge signatures, these were explicitly included in the GHM calibration procedure as this gives better signature simulations than using traditional global objective functions . Calibrating against river flow data alone can lead to unrealistic snow and glacier melt rates, inhibiting model consistency and increasing projection uncertainties . Accordingly, a novel signature-based calibration of the GHM was undertaken by evaluating the GHM against 20 of the river discharge signatures in Table 1 for which observation data exist, calculated from hourly river discharge measurements at the automatic stream gauge (ASG1 in Fig. 1) in combination with 12 signatures of ice melt and snow coverage (Appendix Appendix B).

For each signature, model simulations were compared to observations using a continuous acceptability score that is analogous to those used in other signature-based hydrological studies . This objective function explicitly accounts for uncertainty in the observation signatures, hereafter termed “limits of acceptability” (LOA), so that decisions about model appropriateness can be made within the uncertainties of observation data. In this study the 95 % confidence bounds were used to define the LOA for the river discharge signatures (Table 1) and the ice melt and snow coverage signatures (Table B1). Details of how these were derived can be found in the study of . The acceptability for signature j is defined as

$\begin{array}{}\text{(1)}& {s}_{j}=\left\{\begin{array}{ll}\mathrm{0}& {\mathrm{low}}_{j}\le {\mathrm{sim}}_{j}\le {\mathrm{upp}}_{j}\\ \frac{{\mathrm{sim}}_{j}-{\mathrm{upp}}_{j}}{{\mathrm{upp}}_{j}-{\mathrm{obs}}_{j}}& {\mathrm{sim}}_{j}>{\mathrm{upp}}_{j}\\ \frac{{\mathrm{sim}}_{j}-{\mathrm{low}}_{j}}{{\mathrm{obs}}_{j}-{\mathrm{low}}_{j}}& {\mathrm{sim}}_{j}<{\mathrm{low}}_{j}\end{array}\right\,\end{array}$

where obsj and simj are the observed and simulated values, and uppj and lowj are the upper and lower LOA. A score of zero indicates that the model captures the signature within the LOA. A non-zero score is given for any simulation that falls outside of the LOA with a sign that indicates the direction of bias and a magnitude that indicates the model's performance relative to the LOA. A score of −3 would indicate that the model underestimates the signature by 3 times the observation uncertainty. This score therefore does not penalise a model if it falls within the observation uncertainty of a signature. It is also tolerant of projections that fall outside of the LOA where observation uncertainty is high – a desirable attribute given the range of signatures the GHM was evaluated against.

The aim of the calibration was to extract an ensemble of GHM compositions (TIM and ROR structure–parameter combinations) that were most acceptable across the river discharge signatures whilst broadly reproducing the snow coverage and ice melt signatures. This was achieved using a two-stage Monte Carlo procedure which was devised so that the resultant GHM ensemble reflected the uncertainty in model selection given the known inconsistencies of the GHM across the signatures.

### 2.5.1 Stage 1: TIM calibration

The first stage aimed to extract the optimal TIM compositions (structure–parameter combinations) by calibrating them against the 12 snow coverage and ice melt signatures. Here, for each of the three TIM structures, 5000 TIM parameter sets were drawn from predefined uniform distributions (Table C1 in the Appendix) using the quasi-random Sobol sampling strategy to sample the parameter space as efficiently as possible. For each parameter set, the GHM was spun up for 3 years from 1985 to 1988 with a static ice geometry fixed to a 1988 ice digital elevation model (DEM) . The GHM was then run from 1988 to the end of 2016 with a freely evolving glacier geometry.

Given the high degree of glaciation in the study catchment, and its recent rapid retreat, an initial emphasis of the calibration was put on the model's ability to capture the long-term glacier volume change signature. Accordingly, only those TIM compositions that captured this signature within the LOA were considered and the rest were discarded. These remaining compositions were then further refined by evaluating them against the remaining 11 snow and ice signatures. First, the TIM compositions were sorted by structure (TIM1, TIM2, TIM3). Then, for a given TIM structure, the following steps were applied:

1. Find the TIM parameter set(s) that capture the signature within the LOA and discard the rest. If more than one parameter set captures the current signature, go to step 2. If none capture the current signature, discard none and go to step 2.

2. Of the remaining models, find that which best captures the 10 remaining snow and ice signatures overall according to the weighted mean scores obtained in Eq. (1). The weights were applied to ensure that equal preference was given to ice melt and snow coverage signatures.

Twenty-four unique TIM compositions were obtained from this calibration stage made up of eight unique parameterisations of each of the three TIM structures. In some cases the same composition was selected more than once, which was accounted for by weighting the simulations in the results presented throughout this study.

### 2.5.2 Stage 2: ROR calibration

The second calibration stage aimed to extract the optimal ROR compositions when used in combination with the 24 preselected TIM compositions by calibrating them against 20 of the river discharge signatures obtained from observations of river discharge for the years 2013 and 2014 (see signatures with calibration LOA in Table 1). Note that the inter-annual flow signatures and the mean December river flow signatures could not be calculated as there was insufficient observation data. Furthermore, the mean annual river flow and mean monthly flow range were not included as this information was already accounted for in the mean monthly flow signatures. Here, 5000 random ROR parameter sets were drawn for each ROR structure. Each was used in combination with the preselected TIM compositions in the GHM. Then, the two steps outlined in calibration stage 1 were applied using the 20 calibration river discharge signatures with two notable differences. Firstly, for each ROR structure and each river discharge signature, rather than selecting a unique ROR parameter set for each of the 24 TIM compositions, a single parameter set was selected based on its mean performance across the 24 TIM compositions. This was done to satisfy the ANOVA requirements so that the TIM and ROR composition uncertainty could be analysed separately. Furthermore, for step 2, the signatures were weighted so that each of the attributes in Table 1 were weighted equally. In total, 14 unique ROR compositions made up of seven unique parameterisations of the ROR1 and ROR2 structures were selected, giving a total of $\mathrm{24}×\mathrm{14}=\mathrm{336}$ unique GHM compositions.

## 2.6 ANOVA uncertainty analysis

For the 21st century runs, all 336 GHM compositions were run to the end of 2016 using the historic observed climate to capture the evolving ice geometry as accurately as possible. From 2017 to 2100, the 280 downscaled future climate time series were used to drive the GHM compositions resulting in 94 080 individual model runs. For each model run, projections of watershed snow and ice coverage and the 25 river discharge signatures were extracted for six 21st century 25-year time slices centred on the 2030s (2023–2047), 2040s (2033–2057), 2050s (2043–2067), 2060s (2053–2077), 2070s (2063–2087) and 2080s (2073–2097). Future changes in these were then calculated relative to a reference 25-year period (1991–2015). This reference period was chosen because ice coverage data (used to initialise the GHM) were only available from 1988 and historic climate data were available up to the end of 2016. ANOVA was used to quantify the effect size of the five components of the model chain, hereafter termed “factors”, on each signature for each 21st century time slice. Note that the peak water (PW) signature can only be calculated taking into account the full projection time series and, as such, it was not possible to apply ANOVA to each time slice for this signature. The five factors include the 2 future ESs, 14 GCM–RCM combinations, 10 DS parameterisations, 24 TIM compositions and 14 ROR compositions. ANOVA offers an intuitive approach to estimate the effect size of each factor on each signature by partitioning the total sum of squares (SStot) in the response variable over all combinations of factor levels:

$\begin{array}{}\text{(2)}& {\mathrm{SS}}_{\mathrm{tot}}={\mathrm{SS}}_{a}+{\mathrm{SS}}_{b}+{\mathrm{SS}}_{c}+{\mathrm{SS}}_{d}+{\mathrm{SS}}_{e}+{\mathrm{SS}}_{\mathrm{I}}+{\mathrm{SS}}_{\mathit{\epsilon }},\end{array}$

where

$\begin{array}{}\text{(3)}& {\mathrm{SS}}_{\mathrm{tot}}=\sum _{i=\mathrm{1}}^{{n}_{a}}\sum _{j=\mathrm{1}}^{{n}_{b}}\sum _{k=\mathrm{1}}^{{n}_{c}}\sum _{l=\mathrm{1}}^{{n}_{d}}\sum _{m=\mathrm{1}}^{{n}_{e}}\left({y}_{i,j,k,l,m}-\stackrel{\mathrm{‾}}{Y}{\right)}^{\mathrm{2}},\end{array}$

where na, nb, nc, nd and ne are the number of levels for each factor, y is the response for a given treatment (i.e. combination of factor levels) and $\stackrel{\mathrm{‾}}{Y}$ is the grand mean of the response variable over all treatments. SSa, SSb, SSc, SSd and SSe in Eq. (2) are the sum of squares due to the main effects, i.e. the variability in the response variable due to varying a given factor on its own. For example,

$\begin{array}{}\text{(4)}& {\mathrm{SS}}_{a}={n}_{b}{n}_{c}{n}_{d}{n}_{e}\sum _{i=\mathrm{1}}^{{n}_{a}}\left({y}_{i,\circ ,\circ ,\circ ,\circ }-\stackrel{\mathrm{‾}}{Y}{\right)}^{\mathrm{2}},\end{array}$

where indicates averaging over an index. SSI includes all nonadditive interaction terms where the combined effect of two or more factors is not the sum of their main effects. For a five-factor ANOVA, one could include all unique n-tuple combinations of factors where $n=\left(\mathrm{2},\mathrm{3},\mathrm{4},\mathrm{5}\right)$. Given the difficulty in interpreting these higher-order interactions, and computational requirements, it was decided to investigate the nine first-order interactions only, so that

$\begin{array}{ll}{\mathrm{SS}}_{\mathrm{I}}& ={\mathrm{SS}}_{ab}+{\mathrm{SS}}_{ac}+{\mathrm{SS}}_{ad}+{\mathrm{SS}}_{ae}+{\mathrm{SS}}_{bc}+{\mathrm{SS}}_{bd}+{\mathrm{SS}}_{be}\\ \text{(5)}& & +{\mathrm{SS}}_{cd}+{\mathrm{SS}}_{ce}+{\mathrm{SS}}_{de}.\end{array}$

The sum of squares for a first-order interaction are calculated as follows using factors a and b as an example:

$\begin{array}{}\text{(6)}& {\mathrm{SS}}_{ab}={n}_{c}{n}_{d}{n}_{e}\sum _{i=\mathrm{1}}^{{n}_{a}}\sum _{j=\mathrm{1}}^{{n}_{b}}\left({y}_{i,j,\circ ,\circ ,\circ }-{y}_{i,\circ ,\circ ,\circ ,\circ }-{y}_{\circ ,j,\circ ,\circ ,\circ }+\stackrel{\mathrm{‾}}{Y}{\right)}^{\mathrm{2}}.\end{array}$

Finally, the SSε term includes all unexplained variance, i.e. error in the ANOVA model.

Having partitioned the sum of squares, the effect size, η2, for any term in Eq. (3) can be taken as the proportion of the total sum of squares:

$\begin{array}{}\text{(7)}& {\mathit{\eta }}_{*}^{\mathrm{2}}={\mathrm{SS}}_{*}/{\mathrm{SS}}_{\mathrm{tot}},\end{array}$

where * can be any of the main effects, interactions or error term.

showed that because ANOVA is based on a biased variance estimator that underestimates the variance in small sample sizes, the calculated effect sizes are biased if a different number of levels are used for each factor. Given that the number of factor levels ranges from 2 to 24, a pure application of ANOVA using all possible treatments would lead to biased results. outlined a method to correct for this which involves subsampling the factor levels down to the smallest number levels across all factors. The procedure is repeated using every possible combination of factor levels with unbiased effect size taken as the mean across all subsamples. However, given that there are >108 unique combinations of factor levels when subsampled down to two (and discarding factor level repetitions), it would have been infeasible to account for every possible combination. Instead, it was decided to calculate the effect sizes in this manner using five different subsample sizes (101,102 … 105). The results were then analysed to see if the effect sizes converged. It was found that 103 subsamples were sufficient to converge the effect sizes for all river discharge signatures and projections of snow and ice coverage. Accordingly, this subsampling strategy was adopted in this study.

3 Results

## 3.1 Evaluation of calibrated GHM compositions

The simulated river discharge time series and signatures using the calibrated GHM compositions were evaluated against river discharge observations covering the years 2015 and 2016. Note that no data for mean January and February flows were available for these years. Figure 4a and b show the simulated “capture ratio” (the ratio of the 336 GHM compositions that capture the observation data within their 95 % uncertainty bounds) time series projected onto the mean observed river discharge for the years 2015 and 2016 respectively. Also shown is the ensemble mean simulated river discharge (dashed black line) which, while not indicative of a single GHM simulation, does provide an indication of overall projection bias.

Figure 4Capture ratio projected onto observed river discharge data during evaluation period for 2015 (a), 2016 (b), and over the FDC (c). The weighted ensemble mean simulation is shown as a dashed black line. Also shown are the range of acceptability scores for each of the available river discharge signatures over the evaluation period (d). Acceptable simulations in (d) are those contained within the dashed black lines.

A total of 56 % of the observation time series were captured by at least half of the GHM compositions, while 41 % and 28 % of the observations were captured by at least 75 % and 90 % of the GHM compositions. A total of 12 % of the observations were not captured by any of the GHM compositions. These included some of the low flows observed at the beginning of the year outside of the melt season, particularly in 2015, where the GHM showed consistent negative biases. Some rainfall-induced summer peak flows were also not captured, particularly during the late summer months of August and September. Furthermore, the sustained summer melt runoff discharge in between rainfall-induced peak flows tended to be overestimated (for example during July and August 2016). Even so, the flow duration curve in Fig. 4c shows that almost the entire FDC was captured by all of the GHM simulations except for some of the lowest flows on record. Indeed, Fig. 4d reveals that GHMs were least efficient at capturing the low-flow signatures, particularly the variability signature (σ99–95), where simulations were positively biased by almost 4 times the observation uncertainty. For the remaining signatures, though, the ensemble of GHMs was remarkably efficient, with the majority of simulations (and in most cases all of them) capturing these signatures within their 95 % observation uncertainty bounds.

Figure 5Seasonal average projected changes in ECDFs for near-surface air temperature (a, d, g, j), incident solar radiation (b, e, h, k), and total precipitation (c, f, i, l) for the late 21st century (2076–2100) relative to the recent past (1981–2005). Changes are plotted for each 10 % section of the ECDFs. For each section, blue and yellow dots represent each of the 140 downscaled future climate time series for the RCP4.5 and RCP8.5 ES respectively (280 in total). Winter: Dec, Jan, Feb; spring: Mar, Apr, May; summer: Jun, Jul, Aug; autumn: Sep, Oct, Nov.

## 3.2 Future climate projections

Projections of temperature for the late 21st century (2076–2100) consistently show an increase relative to the recent past (1981–2005). The largest increases are projected for the coldest days of the year during the winter (Fig. 5a), spring (Fig. 5d) and autumn (Fig. 5j) months as shown by the positive skew in the lower sections of the ECDFs. However, these changes are also typically associated with the greatest projection uncertainty. RCP4.5 projects annual mean near-surface air temperature to rise by between 1.1 and 3.6 C by the late 21st century relative to the recent past with an ensemble mean projection of +2.0C. RCP8.5 projects an equivalent rise of between 2.3 and 4.9 C with an ensemble mean projection of +3.3C.

Projected changes in incident solar radiation span positive and negative values, but the median projections are consistently negative, indicating reductions in incident solar radiation are most likely. Uncertainty in the magnitude of change is highest during the spring and summer months (Fig. 5e and h) when incident solar radiation peaks. Under RCP4.5 annual mean incident solar radiation is projected to change by between −10.7 and +0.8 % by the late 21st century with an ensemble mean projection of −4.4 %. Under RCP8.5, changes of between −15.3 and 0.4 % are projected with an ensemble mean projection of −7.7 %.

Projected changes in total precipitation are negligible for the four lowest 10 % sections of the precipitation ECDFs but significant for the two highest sections. In the winter (Fig. 5c) and autumn (Fig. 5l) months, absolute changes exceed 40 mm d−1. The direction of change is uncertain apart from autumn where median projections are consistently positive for the upper sections of the ECDF. The magnitude of change is also uncertain. RCP4.5 projects annual mean precipitation will change by between −13.5 and +21.6 % relative to the recent past by the late 21st century with an ensemble mean projection of +1.7 %. Under RCP8.5, changes of between −25.7 and 25.1 % are projected with an ensemble mean projection of +1.4 %.

Figure 6 shows the correlation matrix calculated between seasonal average climate variables for the late 21st century. For all climate variables, between-season changes (scores within green borders in Fig. 6) are positively correlated, indicating that an increase in summer temperature typically corresponds with an increase in winter temperature for example. Temperature has the greatest between-season correlation while precipitation is the least well correlated. Within-season, between-variable correlation scores (within purple border in Fig. 6) show that precipitation and incident solar radiation are negatively correlated and that the correlation between precipitation and temperature depends on the time of year. For the cooler winter, spring and autumn months, temperature and precipitation are positively correlated, but there is a weak negative correlation for the summer months. Temperature and incident solar radiation are negatively correlated, most strongly for the cooler winter, spring and autumn months.

Figure 6Correlation matrix between seasonal average climate variables calculated for the late 21st century (2076–2100) using the 280 downscaled future climate time series. Within-variable, between-season correlation scores are contained within the green borders and within-season, between-variable correlation scores are contained within the purple borders. Those regions of the correlation matrix that do not cover these two groups are shaded in black.

## 3.3 Future evolution of snow and ice coverage

The ensemble mean projections of annual mean watershed snow coverage (Fig. 7a) show that it will decrease from 12.2 km2 in 2016 to 9.2 km2 in 2100 (25 % reduction) under RCP4.5 and 6.0 km2 (51 % reduction) under RCP8.5. The 95 % projection confidence intervals indicate that by 2100 the watershed could be almost entirely free of snow (2.5 km2 remaining) under RCP8.5 or could have a coverage exceeding present levels (13.3 km2) under RCP4.5.

Figure 7Projected annual mean watershed snow coverage (a) and ice coverage (b) including the projection confidence intervals (bands) and ensemble mean projections (thick solid lines) for the RCP4.5 (blue) and RCP8.5 (yellow) projections. Also shown are projection confidence levels for a reduction in coverage relative to 2016 (thin solid lines, right-hand axis).

Beyond 2050, there is high confidence (≥95 %) that snow coverage will reduce relative to 2016 levels under RCP8.5 (thin yellow line in Fig.  7a) and equally high levels of confidence apply to projected reductions in snow coverage beyond 2066 under the cooler RCP4.5 (thin blue line in Fig. 7a).

The ensemble mean projection of ice coverage (Fig. 7b) projects a 31 % reduction relative to 2016 by 2100 under RCP4.5 and a more severe 63 % reduction under RCP8.5. There is high confidence (≥95 %) that ice coverage will be less than 2016 levels from 2037 onwards under RCP4.5 and from 2030 onwards under RCP8.5, but the magnitude of change is uncertain under both emission scenarios. By 2100, the 95 % confidence intervals for both RCP4.5 and RCP8.5 are 6.5 km2 wide (more than half the 2016 watershed ice coverage).

The simulation that projected the minimum ice coverage by 2100 under the RCP8.5 emission scenario shows sustained retreat of the glacier between 2000 and 2100 resulting in a watershed that is almost entirely ice-free by the end of the century (Fig. 8). In contrast, the maximum ice coverage simulation under the RCP4.5 emission scenario projects two periods of glacier advance between 2010 and 2030 and between 2060 and 2100. By the end of the century, this simulation projects ice coverage will be similar to that in 2000.

Figure 8Simulated ice thickness between 2000 and 2100 based on simulations that projected the maximum (RCP4.5) and minimum (RCP8.5) ice coverage by 2100. Watershed outline shown in magenta.

Figure 9a, b and c show the climate projection time series that produced the minimum (dotted lines) and maximum (dashed lines) snow (blue lines) and ice (red lines) coverage by 2100. The minimum coverage simulations were forced with some of the highest temperature time series while the maximum coverage simulations were forced with some of the lowest. The maximum coverage simulations show higher-than-average incident solar radiation inputs (Fig. 9b) and lower precipitation volumes than the minimum coverage simulations. Indeed, correlation scores calculated between seasonal average climate variables and the simulated snow and ice coverage by 2100 (Fig. 9d) show that there is a strong negative correlation between mean temperature and projected snow and ice coverage and a weaker positive correlation between snow and ice coverage and incident solar radiation. An even weaker negative correlation exists between autumn and winter precipitation and snow and ice coverage.

Figure 9Relationship between driving climate data and projected snow and ice coverage including annual mean downscaled climate time series of temperature (a), incident solar radiation (b), and total precipitation (c) with time series that produced the minimum (dotted lines) and maximum (dashed lines) snow and ice coverage by the end of 2100. Also included are correlation scores calculated between seasonal average climate variables over the entire future period (2017–2100) and simulated snow and ice coverage by the end of 2100 (d).

## 3.4 Sources of uncertainty in snow and ice coverage projections

The effect size of the main, interaction and error terms calculated using ANOVA for projected changes in snow and ice coverage is shown in Fig. 10. Note that ROR effects are not included here as this model chain component has no influence on cryospheric processes in the GHM. The effect size of each ANOVA term changes through the decades and also varies between snow and ice coverage. Throughout the 21st century, TIM uncertainty contributes <3 % to the total projection uncertainty of snow coverage. For projections of ice coverage, ${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}$ is greater than 0.11 up to and including the 2060s, but then gradually falls to 0.07 by the 2080s. ${\mathit{\eta }}_{\mathrm{DS}}^{\mathrm{2}}$ and ${\mathit{\eta }}_{\mathrm{I}}^{\mathrm{2}}$ never exceed 0.1 for snow and ice coverage, and as with ${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}$, they gradually reduce through the latter half of the 21st century. GCM–RCM uncertainty is the largest contributor to ice coverage projection uncertainty in the 2030s with an effect size of 0.47. For snow coverage, ES and GCM–RCM have similar effect sizes of 0.45 and 0.4 respectively. However, for the mid- and latter half of the 21st century ES uncertainty dominates, contributing 73 % and 72 % of snow and ice coverage total projection uncertainty by the 2080s respectively.

Figure 10Effect size (η2) of main effects (ES, GCM–RCM, DS and TIM), interactions (I), and remaining error (ϵ) on projected changes in snow and ice coverage calculated using ANOVA for the six 21st century time slices. Note that the ROR main effect is not included here as it does not influence cryospheric processes in the GHM.

## 3.5 Future evolution of primary runoff components

As an initial indication of the potential for downstream river flow regime change, Fig. 11 shows the 21st century evolution of changes in the four primary runoff components relative to the reference period. The ensemble means (solid coloured lines) indicate that by the end of the century rainfall will increase for all months under both emission scenarios except for August, where RCP8.5 shows a small decrease in rainfall on average (Fig. 11a). The largest increases are shown during the autumn (SON) and winter (DJF) months under RCP8.5. The confidence in the direction of change by the end of the century is ≥90 % for 6 months under RCP8.5 (as indicated by the coloured bands) but only for 2 months (March and April) under RCP4.5. However, ≥75 % of the projections from both RCPs project an increase in rainfall between October and April (as indicated by the markers in Fig. 11a). A comparison of the reference and 2080s monthly ensemble means (inset in Fig. 11a) indicates that the peak rainfall month will shift from September to October.

Figure 11Projections of monthly mean runoff components including rainfall (a), snowmelt (b), ice melt (c), and evapotranspiration (d) for RCP4.5 (blue) and RCP8.5 (yellow). For each month, the trajectory of the ensemble mean change over the 21st century time slices (2030s to 2080s) relative to the reference period (1991–2015) is shown by the solid coloured lines. These lines are marked for each time slice where there is ≥75 % confidence in the direction of change. They are bounded by the 10th and 90th percentiles of the projections (bands). Inset in each plot are ensemble mean monthly runoff volumes averaged over the reference period (black solid line) and 2080s (dashed lines).

For snowmelt, the greatest changes are projected to occur in the summer months of July and August under RCP8.5 where there is ≥90 % confidence that melt will reduce relative to the reference period from the 2040s onwards (Fig. 11b). RCP4.5 also projects decreases in summer melt, but the magnitude of change is smaller. In the winter months, both RCPs project a small increase in melt on average by the end of the century. The ensemble means project that total summer (JJA) melt will reduce by 19 % under RCP4.5 and 37 % under RCP8.5 by the 2080s (inset in Fig. 11b). Annual melt will decrease by 12 % under RCP4.5 and 26 % under RCP8.5. A similar pattern of change is projected for ice melt (Fig. 11c) where total summer (JJA) melt will reduce by 33 % under RCP4.5 and 58 % under RCP8.5 by the 2080s. There is high confidence (≥90 %) that mean monthly ice melt will reduce for all months except December under RCP8.5. Under RCP4.5 a small increase in winter ice melt is projected for the early and mid-21st century, but by the 2080s winter melt is projected to reduce near to reference levels on average. Under RCP8.5, winter ice melt is projected to reduce relative to reference levels for the latter half of the 21st century.

Projections consistently (≥90 %) show an increase in evapotranspiration for all months of the year (Fig. 11d) with the largest increases projected under RCP8.5 towards the end of the 21st century. However, the volume of evapotranspiration will remain a small component of the overall water balance.

## 3.6 Future evolution of river flow regime

Figure 12 shows the projected changes in river discharge signatures relative to the reference period (except peak water for which the raw projections are shown). Under RCP4.5 the ensemble mean projection of peak water is 2045, while under RCP8.5 peak water is projected to occur 17 years earlier in 2028. Indeed, the RCP8.5 projections of the mean annual flow signature ($\stackrel{\mathrm{‾}}{Q}$) show a consistent decline through the 21st century with ≥90 % confidence that flows will reduce by the end of the century by 19 % on average. In contrast, under RCP4.5 the magnitude of the decline is smaller (ensemble mean projects a 7 % decrease for 2090s) and the direction of change is more uncertain. Both RCPs project an increase in inter-annual flow range (RANN) throughout the 21st century (≥75 % under RCP8.5). Under RCP4.5 the ensemble mean projects a 47 % increase in RANN by the 2080s while RCP8.5 shows a 71 % increase.

Figure 12Projected changes in river discharge signatures. For each signature, the trajectory of the ensemble mean change over the 21st century time slices (2030s to 2080s) relative to the reference period (1991–2015) is shown by the solid coloured lines. These lines are marked for each time slice where there is ≥75 % confidence in the direction of change. They are bounded by the 10th and 90th percentiles of the projections (bands). Also shown are 2080s ensemble mean change expressed as a percentage of simulated signatures for the reference period (text). Note that the peak water (PW) signature is not expressed as a change but as the overall raw projections.

Seasonally, monthly winter (DJF) flows are projected to increase while ≥90 % of the ensemble projects a decrease in summer (JJA) flows by the 2090s under both RCPs. The absolute change in mean monthly flows is larger for summer flows on average, but proportionally the winter flows are projected to change most, particularly in February where the ensemble mean projects a 60 % and 67 % increase under RCP4.5 and RCP8.5 respectively by the end of the century. The combined effect of increased winter flows and decreased summer flows results in decreased intra-annual flow variability. Under both RCPs, more than 90 % of the ensemble projects a decrease in Rmnth relative to the reference period from the 2050s onwards. The ensemble mean projections under RCP8.5 show a decade-on-decade reduction in Rmnth with time and a 41 % reduction by the end of the century.

Of those signatures with units m3 s−1, the high-flow Q01 signature shows the largest ensemble mean increase of 2.8 and 2.5 m3 s−1 for RCP4.5 and RCP8.5 respectively by the end of the century. There is high confidence (≥75 %) that Q01 will increase relative to the reference period under RCP8.5, but the magnitude of change is uncertain under both RCPs. For Q05, the ensemble means from both RCPs both show a reduction throughout the 21st century; however, the 10th and 90th percentile span positive and negative values of change for all decades. The ensemble mean projections of changes to high-flow variability (σ05–01) are positive throughout the 21st century under both RCPs. In the latter half of the century, ≥75 % of the projections under RCP4.5 show an increase in σ05–01 while ≥90 % of the projections under RCP8.5 show an increase.

For moderate flows, the ensemble mean of the RCP4.5 projections shows a small reduction in Q50 of approximately 0.15 m3 s−1 throughout the 21st century while the RCP8.5 ensemble mean projects a decade-on-decade reduction in Q50, and by the end of the century there is high confidence (≥90 %) that moderate flows will reduce under this emission scenario. Moderate-flow variability (σ52–48) is projected to reduce with high confidence under both RCPs, albeit by only 0.03 and 0.06 m3 s−1 by the 2080s under RCP4.5 and RCP8.5 respectively.

For the slow-release low-flow signatures, ≥90 % of the projections are positive throughout the 21st century under both RCPs, indicating an increase in the magnitude of low-flow events (or equivalently a reduction in the frequency of these flow events) and variability of low flows. The absolute changes in the ensemble means never exceed 0.1 m3 s−1 for these signatures, although proportionally they show the largest degree of change, particularly for Q99 where the proportional change exceeds 2000 % under RCP4.5.

Finally, the response time to runoff (τ) is projected to decrease throughout the 21st century under both RCPs (≥90 % confidence), indicating the catchment will likely become more flashy. The magnitude of change is small where the ensemble mean projects a small reduction in τ of 3.9 h under RCP4.5 and a slightly greater reduction of 4.7 h under RCP8.5.

## 3.7 Sources of uncertainty in river flow regime projections

Figure 13 shows the ANOVA effect sizes calculated for the 2030s and 2080s for each river discharge signature. The error term (${\mathit{\eta }}_{\mathit{ϵ}}^{\mathrm{2}}$) never exceeds 0.09 and for 21 of the 25 signatures it is <0.03, indicating that the main effects and first-order interaction terms explain the majority of projection uncertainty. For the 2030s, ES uncertainty contributes 4 %–27 % of the total projection uncertainty across the signatures. By the 2080s, ES contributes up to 65 % of the total projection uncertainty. In fact, for all but four signatures, ES contributes proportionally more to total projection uncertainty in the 2080s than the 2030s. By the 2080s the five signatures with the highest ${\mathit{\eta }}_{\mathrm{ES}}^{\mathrm{2}}$ include the mean monthly flows from May to August and the mean monthly flow range (Rmnth) signature (Table 2) for which the effect sizes are at least 0.47. GCM–RCM uncertainty is the largest contributor to total projection uncertainty for 21 of the 25 river discharge signatures for the 2030s, and it still remains a significant contributor to projection uncertainty by the 2080s with a mean effect size across the signatures of 0.3. Four of the five most sensitive signatures to GCM–RCM uncertainty for the 2030s remain in this top five for the 2080s (Table 2), and these include the mean monthly winter flows in January and February and two of the quick-release high-flow signatures (Q01 and Q05).

Figure 13Effect size (η2) of all main effects (ES, GCM–RCM, DS, TIM and ROR), interactions (I) and remaining error (ϵ) on projected changes in the 25 river discharge signatures at the start (2030s, a) and end (2080s, b) of the 21st century.

Table 2Top five river discharge signatures ranked according to the average effect size for each of the main effects, interactions, and remaining error on projected changes for the 2030s and 2080s. Effect sizes are included in brackets.

On average, the DS parameterisation contributes 18 % of the total projection uncertainty across the signatures for the 2030s. In fact, ${\mathit{\eta }}_{\mathrm{DS}}^{\mathrm{2}}$ is relatively consistent across the signatures, ranging from 0.1 to 0.2 for 18 of the 25 signatures. For the 2080s, ${\mathit{\eta }}_{\mathrm{DS}}^{\mathrm{2}}$ reduces for all signatures except mean November and December flows and the inter-annual flow range (RANN). For RANN, DS has the largest effect size, contributing 43 % of the total projection uncertainty. Autumn and winter monthly mean flows for September, November, December and February make up the remainder of the top five signatures most affected by DS uncertainty for the 2080s. On average TIM uncertainty contributes 9 % of the total projection uncertainty across the different signatures for the 2030s. For this period it is the largest contributor to RANN uncertainty (${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}=\mathrm{0.35}$) and it also shows significant contributions to mean monthly flow projection uncertainty for April (${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}=\mathrm{0.17}$) and May (${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}=\mathrm{0.23}$) at the beginning of the melt season. It is also the largest contributor to uncertainty of projections of response time to runoff (τ) where ${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}=\mathrm{0.20}$. For the 2080s the average ${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}$ falls slightly to 7 %, but TIM uncertainty remains an important contributor to total projection uncertainty for τ, April and May flows, and two of the low-flow signatures (Q95 and σ99–95) where ${\mathit{\eta }}_{\mathrm{TIM}}^{\mathrm{2}}\ge \mathrm{0.12}$. Uncertainty stemming from the ROR model chain component has a negligible influence on the decadal signatures (PW and RANN) or those signatures characterising annual and monthly mean flows for the 2030s and 2080s. For the 2030s, ROR uncertainty is important for projections of low-flow magnitude (Q99 and Q95, ${\mathit{\eta }}_{\mathrm{ROR}}^{\mathrm{2}}=\mathrm{0.43}$ and 0.20 respectively) and variability (σ99–95, ${\mathit{\eta }}_{\mathrm{ROR}}^{\mathrm{2}}=\mathrm{0.13}$). In fact, for Q99, ROR is the single largest contributor to total projection uncertainty. For the 2080s, its influence on low-flow quantiles remains significant and it is the single largest contributor to both the Q99 and τ projection uncertainty. It also remains a significant contributor to the high-flow variability signature, σ05–01, where ${\mathit{\eta }}_{\mathrm{ROR}}^{\mathrm{2}}=\mathrm{0.12}$.

Unlike ice and snow coverage, interactions between model components significantly contribute to the total projection uncertainty across the signatures where ${\mathit{\eta }}_{\mathrm{I}}^{\mathrm{2}}$ ranges between 0.07 and 0.27 for the 2030s and between 0.07 and 0.32 for the 2080s. Figure 14 shows the decomposition of the five interaction terms with the largest effect sizes on average for the 2030s (a) and 2080s (b). The interactions between the ES and GCM–RCM model chain components dominate the contribution to projection uncertainty. However, interactions between the climate model chain components and the GHM (e.g. DS–TIM) may also contribute to the projection uncertainty. For RANN, DS–TIM interaction contributes 7 % of the total projection uncertainty for the 2030s and 2080s. Furthermore interactions between the TIM and ROR in the GHM contribute some (albeit small) amounts to the total projection uncertainty. For 16 of the signatures, the contribution from interactions between model chain components increases from the 2030s to the 2080s. These include all of the signatures that characterise, high-, moderate- and low-flow magnitude and variability, but the largest increases are shown for March and October mean monthly flows.

Figure 14Effect size (η2) of the five most significant interactions on projected changes in the 25 river discharge signatures at the start (2030s, a) and end (2080s, b) of the 21st century.

4 Discussion

## 4.1 Future evolution of river flow regime

There is high confidence that near-surface air temperature will increase by the late 21st century (2076–2100) relative to conditions in the recent past (1981–2005). Precipitation and incident solar radiation were projected to slightly increase and decrease respectively on average – a finding that is consistent with other analyses of the EURO-CORDEX projections for northern Europe . The primary driver of changes in snow and ice is near-surface air temperature, while precipitation and incident solar radiation are of secondary importance. Because of this, there is high confidence that glacier ice and snow will continue to retreat as near-surface air temperature rises throughout the 21st century, which would leave the river basin almost free of snow and ice by 2100 under the warmest climate projections.

The signature-based analysis undertaken in this study has revealed how climate change will impact the magnitude, timing and variability of downstream river flows over a range of timescales in the Virkisá river basin. Projected changes in flow seasonality broadly follow those shown for other mid-latitude alpine river basins where the loss of snow and ice will reduce meltwater inputs in the summer months and a phase shift of precipitation from snowfall to rainfall combined with enhanced melt during the colder months will increase winter runoff . Summer runoff is projected to decrease by 24 % under RCP4.5 and 40 % under RCP8.5 by the 2080s while winter runoff is projected to increase by 59 % under RCP4.5 and 57 % under RCP8.5 by the 2080s. The consequence of these seasonal shifts in runoff is that intra-annual (monthly) flow variability will reduce by 25 % (RCP4.5) and 41 % (RCP8.5) by the 2080s. Furthermore, the magnitude of very low flow events (Q99), which typically occur during the winter months, is likely to increase.

On average, the projections indicated that the seasonal redistribution of runoff will have little influence on mean annual flows under RCP4.5 (−7 % by the 2080s) as changes in summer and winter flows approximately compensate for one another. Under RCP8.5, however, the more pronounced reduction in summer melt inputs results in a 19 % reduction by the 2080s. The loss of a consistent melt input to the river basin and its evolution to a hydrological regime governed by rainfall-runoff processes means inter-annual flow variability (RANN) will increase by 47 % (RCP4.5) and 71 % (RCP8.5) by the 2080s. The increase in rainfall inputs, particularly during the storm-prone autumn and winter months, likely explains the projected increased magnitude of very high flow events (Q01), a finding that is in agreement with other studies that have investigated changes in high-flow magnitudes in glaciated river basins . It is likely that the intensification of peak flow magnitudes will be further exacerbated by the projected decrease in river flow response time to runoff (τ), which is an artefact of losing the runoff-regulating ice and snow water stores. Accordingly, the river basin will become more flashy and flood-prone in the future.

Increased flood frequency has major implications for local infrastructure in the vicinity of the Virkisá river basin. In particular, the southern section of the Route 1 highway which passes over the Skeiðarársandur floodplain navigates a large number of glacier-fed rivers including the Virkisá. Due to the unconsolidated nature of the floodplain lithology, the morphology of these rivers can change rapidly, particularly during high-flow events (Marren2005), and often at considerable expense to the road authority . Accordingly, the projected increase in frequency and severity of high-flow events will likely incur further expenses to maintain this transport link in the future.

Beyond local implications, one should be cautious in extrapolating the findings from this study to other glaciated catchments in Iceland or beyond, but it is clear that the timing, magnitude and variability of glacier-fed river flows over a range of timescales are sensitive to climate change. For Iceland, these changes could impact glacier-fed hydroelectric dams, which are a primary source of electricity for the country. Increased frequency and magnitude of high-flow events could render current dams unsafe if their designed flood capacity can no longer meet regulation requirements . The redistribution and levelling out of seasonal flows, however, could actually have a beneficial effect on the running costs and capacity to produce electricity from such projects .

## 4.2 Uncertainties in projections of river flow regime

Projections of the direction of change relative to the reference period were well constrained for the majority of river discharge signatures, particularly towards the end of the 21st century and for the warmer RCP8.5 emission scenario. Even so, there was considerable spread in the projected magnitude of these changes due to uncertainties in the driving climate data (ES, GCM–RCM, DS) and representation of glacio-hydrological processes (TIM, ROR) in the model chain. Uncertainty in future snow and ice coverage primarily stemmed from the ES due to its control on future near-surface air temperature. In fact, the proportional contribution of the ES to projection uncertainties increased throughout the 21st century and, consequently, the ES was also found to be the dominant source of uncertainty for projections of mean monthly flows during the melt season by the 2080s. The growing influence of the ES over time was also shown by for six alpine catchments in Switzerland and by for two mountain river basins in the Tian Shan. Interestingly though, these studies along with the recent study of found that climate model uncertainty was still the dominant source for projections of monthly river flows. postulated that this was likely because of the high uncertainty in future precipitation across the climate models. Indeed, others have also attributed future runoff uncertainty in glaciated river basins to variability in precipitation projections , a finding which is compounded by an increasingly warm and thus rainfall-dominated precipitation input. In this study, however, the GCM–RCM model chain component only dominated river flow projection uncertainty during the winter months while summer flow uncertainty was dominated by the ES. There are two key reasons that could explain this. Firstly, precipitation uncertainty across the GCM–RCM combinations showed to be especially high during winter (Fig. 5), which coupled with the fact that rainfall is the primary source of runoff during winter likely explains the dominant role GCM–RCM plays in projection uncertainty during the winter months. Furthermore, it should be noted that the Virkisá river basin has a much higher proportional glacier coverage (60 %) compared to the aforementioned studies (1.8 %–22.3 %). Therefore, it is postulated that the influence of the ES in the summer is related to the relatively high proportion of melt runoff that the Virkisá river receives during these months and the fact that the ES showed to be the dominant contributor to future ice coverage uncertainty. Importantly, this finding also serves to highlight the need to represent atmosphere–cryosphere–hydrosphere feedbacks adequately in future studies, particularly where glacier coverage is high, through the inclusion of a dynamic glacier evolution model in the model chain like that implemented in this study.

For projections of the inter-annual flow range, the DS procedure was the largest contributor to projection uncertainty by the end of the 21st century, which should be expected given that the perturbation of this procedure accounted for uncertainty in the random year-by-year sampling of the historic climate data. Uncertainty in the TIM structure parameterisation was the dominant contributor to the spread in projections of moderate monthly flows during the transition from the cold to melt season, which corroborates the model comparison study of , who found that the structural representation of melt was important for controlling the initiation of the melt season due to the contrasting sensitivity of the models to temperature and incident solar radiation. also concluded that signatures derived from the flow duration curve as well as those representing flashiness were most sensitive to the configuration of the ROR component of the GHM. Indeed, here it was found that uncertainty in the ROR structure parameterisation significantly contributed to the total projection uncertainty of slow-release low-flow signatures as well as the response time (flashiness) of the catchment to runoff. Similar sensitivities in low-flow metrics to the choice of hydrological model have been shown for non-glaciated river basins by , and they postulated that these might stem from differences in water-storage-release processes in the models. However, a key drawback of this study and other studies that have investigated the role of hydrological model uncertainty in glaciated river basins (e.g. Giuntoli et al.2015; Vetter et al.2017) is that they have implemented multiple model codes and therefore cannot make any definite conclusion about the source of the projection uncertainties. For example, concludes that the sensitivity to the choice of hydrological model could stem from any number of differences between model codes including the structure, climate interpolation method and calibration strategy. In this study, it has been demonstrated that, by using a single but flexible model code, it is possible to separate out the sources of projection uncertainties down to the process level. Such insights can be used to help prioritise those aspects of the GHM that require (i) additional refinement (e.g. through model development) and (ii) adequate representation of their uncertainty to improve projection robustness.

Furthermore, the signature-based analysis undertaken in this study has shown that the importance of these different sources, be it from the GHM or the climate projections, is dependent on which signature of river discharge is being evaluated. It is clear, therefore, that signature-based analyses could be used to help prioritise uncertainty sources based on the characteristic of flow one is interested in. For example, the results from this study indicate that for evaluating changes in monthly melt season runoff only, it may be beneficial to ignore ROR uncertainty and focus time and computational resources on quantifying uncertainties stemming from the remaining model components. In this respect, the time frame of the projections should also be considered, given the apparent change in effect sizes with time demonstrated for projections of snow and ice coverage and river flow signatures (see Appendix Appendix D).

More broadly, the results from this study emphasise the need for impact studies to represent uncertainties stemming from model chain components that control future climate and glacio-hydrological behaviour, the second of which has been widely neglected. The need for this is compounded by the fact that interactions between model chain components exceeded individual main effects for some river discharge signatures. Accordingly, an ensemble that includes perturbations of multiple components of the model chain simultaneously will provide the most rigorous quantification of projection uncertainty.

## 4.3 Limitations

While some characteristics of projected river flow regime change are broadly in agreement with other studies in similar mid-latitude alpine settings (e.g. changes in flow seasonality and projected increase in high-flow magnitude), it is important to emphasise that the projected river flow regime shifts should not be generalised across glaciated river basins. Indeed, recent regional and global studies have shown that local catchment characteristics such as climate and glacier hypsometry largely influence seasonal river flow response to 21st century climate change. In this study a small absolute increase in low-flow magnitude was projected, indicating that climate change and deglaciation could help to limit periods of water scarcity. However, in more arid regions, where rainfall cannot compensate for reductions in melt, the opposite effect has been shown . One might also expect to see much greater changes in the river flow response time to runoff as snow and ice retreat in other river basins. For the Virkisá river basin, a relatively small reduction in response time (τ) was projected on average by the end of the 21st century. This, perhaps, should not be surprising given the small size of the river basin and the fact that previous investigations have shown that Virkisjökull has a well developed conduit drainage system that routes runoff efficiently year-round . For larger river basins with more expansive cryospheric water stores, changes in the response time to runoff could be much greater, substantially increasing the risk of pluvial flooding.

Similar inter-catchment variability should also be expected with regards to the dominant sources of projection uncertainty. Indeed, as already noted in this discussion, some of the results from this study contrast the limited number of studies that have investigated uncertainty sources in other glaciated basins. suggests that catchment elevation influences the importance of the ES on projection uncertainty whereby runoff from higher elevation catchments with more snow and ice is more sensitive to the ES. It is therefore vital that signature-based evaluations like the one undertaken in this study are applied to other glaciated river basins in the future so that regional variations in river flow regime change and uncertainty sources can be evaluated.

It is also important to consider potential deficiencies in the calibrated GHMs. In fact, the model evaluation demonstrated that they were able to capture the majority of the observed river discharge signatures within their observation uncertainty bounds. Even so, it should be noted that there are several limitations in the calibration approach that could have hindered the efficiency of the calibrated models. Firstly, given the distributed structure of the GHM and the fact that it runs on an hourly time step, running the GHM over multiple years required considerable computation time, which limited the number of runs that could be undertaken in the Monte Carlo calibration procedure. A total of 5000 runs was adopted as an appropriate compromise, balancing the density of parameter sampling with available computational resources. Even so, it is recognised that, particularly for the more complex model structures which employ more calibration parameters, a denser parameter sampling could help to find more efficient model parameterisations. It should also be noted that the models were calibrated and evaluated on 4 years of river flow data only. This detail is particularly important given the conceptual nature of the GHM and thus the potential for the calibration parameters to become less applicable when applied to periods outside of the calibration data. Additionally, it is important to highlight possible model deficiencies brought about by the two-step GHM calibration procedure employed, in which the TIM and ROR model chain components were calibrated independently. This was necessary so that the main effects (Eq. 4) and interaction terms (Eq. 6) for both components could be calculated separately (thus achieving the second aim of the study). However, the drawback of implementing this stepwise calibration procedure over one that calibrates both model components simultaneously is that it neglects any interactions between the TIM and ROR models. Of course, its should be noted that the ANOVA results showed that TIM and ROR interactions are negligible except for two of the 25 signatures evaluated.

In the previous model evaluation study undertaken by , they highlighted the historic observed precipitation data as source of model deficiencies. They noted the lack of available precipitation data at higher elevations, making the gridded dataset employed in this study less reliable near the basin summit. They also analysed the effectiveness of the bias-correction procedure applied to the precipitation dataset and showed that it resulted in time series that were well correlated to the AWS data over a 3-day time step but that this correlation degraded at shorter daily and hourly time steps, which could have contributed to the model's inability to capture snow coverage observations higher up in the catchment and river discharge signatures relating to the timing of flows.

Indeed, uncertainties in the historic precipitation data were not included as part of this study, partly because there was almost no information that would have allowed one to quantify these uncertainties (e.g. rain gauge errors), particularly higher up in the catchment where measurements are least reliable. Additionally, though, it would have meant further increasing the size of the model chain ensemble which was already at the very limit of what was computationally feasible. This, however, raises an important broader limitation of the study in that the total projection uncertainties reported are not indicative of the “true” uncertainty. Further insights could undoubtedly be gained by perturbing other model chain components including the historic climate time series and components related to key glacio-hydrological processes such as the snow redistribution routine and glacier evolution model. Certainly, calculated that the bias correction of precipitation contributed up to 22 % of seasonal streamflow projection uncertainty.

Furthermore, the representations of uncertainty in the five components evaluated in this study are themselves not exhaustive. It is well established that uncertainties in climate model ensembles are under-represented , and steps were taken in this study to limit the total ensemble size so that the experiments were computationally feasible. For example, only 10 random DS sequences were generated, and indeed other aspects of the downscaling procedure could have also been modified (e.g. replacing the linear interpolation of change factors with a moving-average model). Additionally, the melt and runoff-routing model structures implemented represent a subset of a much larger population of available models. For example, we adopted simplified energy balance models and the concept of linear reservoirs to route runoff. However, other model structures that employ more complex physically based energy balance approaches and hydraulic models that simulate discrete flow pathways through the glacier (e.g. Arnold et al. 1998) could also be implemented to provide a more accurate representation of the true projection uncertainty.

5 Conclusions

Twenty-first century climate change is projected to alter the magnitude, timing and variability of river flows over decadal to sub-daily timescales in the Virkisá river basin. Relative to the 1990s reference period, there was high confidence in the direction of change for the majority of the 25 river discharge signatures over the 21st century. The magnitude of change, however, was more uncertain. The application of ANOVA demonstrated that the climate model chain components (ES, GCM–RCM, DS) were the main sources of this uncertainty. However, uncertainty relating to glacio-hydrological process representation in the model chain (TIM, ROR) was the dominant source of projection uncertainty for some river discharge signatures. Furthermore, interactions between model chain components can exceed individual main effects. Based on these findings, we make several recommendations for future studies that aim to assess climate change impact on glacier-fed river flows:

1. Studies should seek to evaluate multiple characteristics of river flow regime change (magnitude, timing and variability) over different timescales where possible so that a more thorough understanding of potential environmental and socio-economic impacts can be deduced from projections. Signatures of river discharge provide the ideal tool to evaluate these changes quantitatively. Changes in the magnitude of river flows over decadal to seasonal timescales are already known to be highly site-specific and therefore we should expect that other signatures of regime change will also show considerable inter-catchment variation.

2. Studies should account for uncertainties stemming from both the climate projections and glacio-hydrological process representations so that more robust projections of river flow regime change are produced. The latter has largely been neglected in studies to date.

3. Careful consideration of which model chain components are the dominant sources of projection uncertainty (through the use of methods such as ANOVA) would help to prioritise resources (e.g. computational) to further enhance projection robustness. The results from this study indicate that such decisions will depend on the signatures of river flow regime change that one is interested in projecting.

Data availability
Data availability.

Appendix A: EURO-CORDEX models

A total of 15 unique GCM–RCM combinations using six GCMs and seven RCMs were available to use in this study (Table A1). Figure A1 shows the EURO-CORDEX 0.11 RCM grids. After comparing monthly average simulations from each GCM–RCM over the recent past (1981–2005) against the observed climate data, it was found that the CNRM-CM5–ALADIN53 GCM–RCM has anomalously large negative temperature biases, particularly during the winter months of the year (see red line in Fig. A2d). In addition to this, a root mean squared error (RMSE) score was calculated for each climate variable by comparing monthly observed and simulated empirical distribution functions constructed from catchment-average daily climate data (Fig. A2a–c). When ranked according to their RMSE scores, the CNRM-CM5–ALADIN53 GCM–RCM ranked 14, 13 and 15 out of 15. Given the anomalously high biases in temperature and the importance of temperature for driving hydrological change in the catchment (both in terms of melt rate and the proportion of precipitation falling as rainfall), coupled with the fact that the model was relatively poor across all three climate variables, it was deemed appropriate to remove this model from the ensemble.

Table A1List of GCMs and RCMs used in this study.

Figure A1EURO-CORDEX 0.11 RCM grid lines. RCM nodes are situated at grid line intersects. All RCMs utilise the green grid except for REMO2009 which uses the blue grid.

Figure A2Root mean squared error scores calculated by comparing monthly empirical distribution functions constructed from catchment-average daily observed and simulated (GCM–RCM) total precipitation (a), incident solar radiation (b), and near-surface air temperature (c) data over the recent past (1981–2005). Also shown are the observed and simulated monthly mean near-surface air temperatures for the recent past (d).

Appendix B: Ice melt and snow coverage signatures used for model calibration

Twelve signatures of ice melt and snow coverage which were previously derived by were used for model calibration and are shown in Table B1. These signatures include (i) measurements of winter and summer ice melt in the main ablation zone between 2012 and 2014 which were derived from ablation stake data; (ii) an estimate of long-term glacier volume change calculated using two DEMs of the ice for 1988 and 2011; and (iii) estimates of the average seasonal snow coverage for spring (March and April), early summer (May and June) and late summer (July and August) in the lower (77–587 m a.s.l.), middle (587–776 m a.s.l.) and upper (776–1123 m a.s.l.) sections of the glacier-free basin area. These were calculated from the remotely sensed MOD10A1 MODIS product for the years 2001 to 2015 inclusive .

Table B1Summary of 12 ice melt and snow coverage signatures used to calibrate the GHM with their limits of acceptability. Note that snow coverage is expressed as a proportion of the glacier-free basin area.

Appendix C: GHM calibration parameters

Table C1 lists all of the calibration parameters and their predefined calibration ranges for the melt and runoff-routing model structures used during the GHM calibration procedure. The three melt model structures include the classic temperature-index model (TIM1), the enhanced temperature-index model proposed by Hock (1999) (TIM2) and the enhanced temperature-index model proposed by (TIM3). The two runoff-routing model structures include the single linear reservoir cascade (ROR1) and two linear reservoir cascades in parallel (ROR2).

Table C1Calibration parameters for the melt and runoff-routing model structures.

Appendix D: Decadal changes in effect size for river discharge signatures

Figure D1Effect size of all main effects, interactions and remaining error on projected decadal changes in the 25 river discharge signatures for all future time slices centred on the 2030s to the 2080s.

Author contributions
Author contributions.

JDM undertook all practical elements of this study, including regional climate projection downscaling, GHM calibration, 21st century projections, ANOVA, analysis of results and production of figures. He also led the writing of the manuscript. JE and ARB managed the design, commissioning and operation of the hydro-meteorological monitoring used in this study. All co-authors contributed to the formulation and discussion of methods used as well as writing of the manuscript.

Competing interests
Competing interests.

The authors declare they have no competing interests.

Acknowledgements
Acknowledgements.

This work was supported by a NERC studentship awarded to Jonathan D. Mackay via the Central England NERC Training Alliance (CENTA). Jonathan D. Mackay, Christopher R. Jackson and Jez Everest publish with permission of the Executive Director of the British Geological Survey.

Review statement
Review statement.

This paper was edited by Markus Hrachowitz and reviewed by Massimiliano Zappa and one anonymous referee.

References

Addor, N., Rössler, O., Köplin, N., Huss, M., Weingartner, R., and Seibert, J.: Robust changes and sources of uncertainty in the projected hydrological regimes of Swiss catchments, Water Resour. Res., 50, 7541–7562, https://doi.org/10.1002/2014WR015549, 2014. a, b, c, d, e, f, g

Ali, G., Tetzlaff, D., Soulsby, C., McDonnell, J. J., and Capell, R.: A comparison of similarity indices for catchment classification using a cross-regional dataset, Adv. Water Resour., 40, 11–22, https://doi.org/10.1016/j.advwatres.2012.01.008, 2012. a

Allen, R., Pereira, L., Raes, D., and Smith, M.: Crop evapotranspiration – Guidelines for computing crop water requirements – FAO Irrigation and drainage paper 56, Tech. rep., Food and Agriculture Organization of the United Nations, Rome, Italy, 1998. a

Baraer, M., Mckenzie, J., Mark, B. G., Gordon, R., Bury, J., Condom, T., Gomez, J., Knox, S., and Fortner, S. K.: Contribution of groundwater to the outflow from ungauged glacierized catchments: A multi-site study in the tropical Cordillera Blanca, Peru, Hydrol. Process., 29, 2561–2581, https://doi.org/10.1002/hyp.10386, 2015. a

Bartók, B., Wild, M., Folini, D., Lüthi, D., Kotlarski, S., Schär, C., Vautard, R., Jerez, S., and Imecs, Z.: Projected changes in surface solar radiation in CMIP5 global climate models and in EURO-CORDEX regional climate models for Europe, Clim. Dynam., 49, 2665–2683, https://doi.org/10.1007/s00382-016-3471-2, 2017. a

Beamer, J. P., Hill, D. F., McGrath, D., Arendt, A., and Kienholz, C.: Hydrologic impacts of changes in climate and glacier extent in the Gulf of Alaska watershed, Water Resour. Res., 53, 7502–7520, https://doi.org/10.1002/2016WR020033, 2017. a

Beer, C., Porada, P., Ekici, A., and Brakebusch, M.: Effects of short-term variability of meteorological variables on soil temperature in permafrost regions, The Cryosphere, 12, 741–757, https://doi.org/10.5194/tc-12-741-2018, 2018. a

Björnsson, H. and Pálsson, F.: Icelandic glaciers, Jökull, 58, 365–386, 2008. a

Bliss, A., Hock, R., and Radić, V.: Global response of glacier runoff to twenty-first century climate change, J. Geophys. Res., 119, 717–730, https://doi.org/10.1002/2013JF002931, 2014. a

Bosshard, T., Carambia, M., Goergen, K., Kotlarski, S., Krahe, P., Zappa, M., and Schär, C.: Quantifying uncertainty sources in an ensemble of hydrological climate-impact projections, Water Resour. Res., 49, 1523–1536, https://doi.org/10.1029/2011WR011533, 2013. a, b, c, d, e

Braithwaite, R. J.: Positive degree-day factors for ablation on the Greenland Ice-sheet studied by energy balance modeling, J. Glaciol., 41, 153–160, 1995. a

Brately, P. and Fox, B. L.: Algorithm 659: Implementing Sobol's quasirandom sequence generator, ACM T. Math. Software, 14, 88–100, 1988. a

Brock, B. W., Willis, I. C., and Sharp, M. J.: Measurement and parameterisation of albedo variations at Haut Glacier d'Arolla, Switzerland, J. Glaciol., 46, 675–688, https://doi.org/10.3189/172756500781832675, 2000. a, b

Bunn, S. E. and Arthington, A. H.: Basic principles and ecological consequences of altered flow regimes for aquatic biodiversity, Environ. Manage., 30, 492–507, https://doi.org/10.1007/s00267-002-2737-0, 2002. a

Carey, M., Baraer, M., Mark, B. G., French, A., Bury, J., Young, K. R., and McKenzie, J. M.: Toward hydro-social modeling: Merging human variables and the social sciences with climate-glacier runoff models (Santa River, Peru), J. Hydrol., 518, 60–70, https://doi.org/10.1016/j.jhydrol.2013.11.006, 2014. a

Carvajal, P. E., Anandarajah, G., Mulugetta, Y., and Dessens, O.: Assessing uncertainty of climate change impacts on long-term hydropower generation using the CMIP5 ensemble – the case of Ecuador, Climatic Change, 144, 611–624, https://doi.org/10.1007/s10584-017-2055-4, 2017. a

Casper, M. C., Grigoryan, G., Gronz, O., Gutjahr, O., Heinemann, G., Ley, R., and Rock, A.: Analysis of projected hydrological behavior of catchments based on signature indices, Hydrol. Earth Syst. Sci., 16, 409–421, https://doi.org/10.5194/hess-16-409-2012, 2012. a

Chen, J. and Ohmura, A.: On the influence of Alpine glaciers on runoff, IAHS-AISH P., 193, 117–126, 1990. a

Collins, M., Knutti, R., Arblaster, J., Dufresne, J.-L., Fichefet, T., Friedlingstein, P., Gao, X., Gutowski, W., Johns, T., Krinner, G., Shongwe, M., Tebaldi, C., Weaver, A., and Wehner, M.: Long-term Climate Change: Projections, Commitments and Irreversibility, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., Cambridge University Press, Cambridge, United Kingdom and New York, USA, 1029–1136, 2013. a

Coxon, G., Freer, J., Wagener, T., Odoni, N. A., and Clark, M.: Diagnostic evaluation of multiple hypotheses of hydrological behaviour in a limits-of-acceptability framework for 24 UK catchments, Hydrol. Process., 28, 6135–6150, https://doi.org/10.1002/hyp.10096, 2014. a, b

Daron, J. D. and Stainforth, D. A.: On predicting climate under climate change, Environ. Res. Lett., 8, 034021, https://doi.org/10.1088/1748-9326/8/3/034021, 2013. a

Duethmann, D., Menz, C., Jiang, T., and Vorogushyn, S.: Projections for headwater catchments of the Tarim River reveal glacier retreat and decreasing surface water availability but uncertainties are large, Environ. Res. Lett., 11, 054024, https://doi.org/10.1088/1748-9326/11/5/054024, 2016. a, b

Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., and Savenije, H. H. G.: A framework to assess the realism of model structures using hydrological signatures, Hydrol. Earth Syst. Sci., 17, 1893–1912, https://doi.org/10.5194/hess-17-1893-2013, 2013. a

Farinotti, D., Usselmann, S., Huss, M., Bauder, A., and Funk, M.: Runoff evolution in the Swiss Alps: projections for selected high-alpine catchments based on ENSEMBLES scenarios, Hydrol. Process., 26, 1909–1924, https://doi.org/10.1002/hyp.8276, 2012. a, b

Finger, D., Pellicciotti, F., Konz, M., Rimkus, S., and Burlando, P.: The value of glacier mass balance, satellite snow cover images, and hourly discharge for improving the performance of a physically based distributed hydrological model, Water Resour. Res., 47, W07519, https://doi.org/10.1029/2010WR009824, 2011. a

Flett, V., Maurice, L., Finlayson, A., Black, A. R., MacDonald, A. M., Everest, J., and Kirkbride, M. P.: Meltwater flow through a rapidly deglaciating glacier and foreland catchment system: Virkisjökull, SE Iceland, Hydrol. Res., 1, nh2017205, https://doi.org/10.2166/nh.2017.205, 2017. a

Flett, V. T.: Glacier retreat and projected river regime changes in the hydrologically highly-coupled Virkisjökull catchment, SE Iceland, Doctor of philosophy, University of Dundee, Dundee, Scotland, 2016. a

Fountain, A. G. and Tangborn, W. V.: The Effect of Glaciers on Streamflow Variations, Water Resour. Res., 21, 579–586, https://doi.org/10.1029/WR021i004p00579, 1985. a

Fyke, J. and Matthews, H. D.: A probabilistic analysis of cumulative carbon emissions and long-term planetary warming, Environ. Res. Lett., 10, 115007, https://doi.org/10.1088/1748-9326/10/11/115007, 2015. a

Garee, K., Chen, X., Bao, A., Wang, Y., and Meng, F.: Hydrological Modeling of the Upper Indus Basin: A Case Study from a High-Altitude Glacierized Catchment Hunza, Water, 9, 1–20, https://doi.org/10.3390/w9010017, 2017. a

Gaudard, L., Romerio, F., Dalla Valle, F., Gorret, R., Maran, S., Ravazzani, G., Stoffel, M., and Volonterio, M.: Climate change impacts on hydropower in the Swiss and Italian Alps, Sci. Total Environ., 493, 1211–1221, https://doi.org/10.1016/j.scitotenv.2013.10.012, 2014. a

Giorgi, F., Jones, C., and Asrar, G. R.: Addressing climate information needs at the regional level: The CORDEX framework, World Meteorological Organization Bulletin, 58, 175–183, 2009. a, b

Giuntoli, I., Vidal, J.-P., Prudhomme, C., and Hannah, D. M.: Future hydrological extremes: the uncertainty from multiple global climate and global hydrological models, Earth Syst. Dynam., 6, 267–285, https://doi.org/10.5194/esd-6-267-2015, 2015. a, b, c

Gosseling, M.: CORDEX climate trends for Iceland in the 21st century, Tech. rep., Icelandic Meteorological Office, Reykjavik, Iceland, 2017. a

Griffiths, J., Keller, V., Morris, D., and Young, A. R.: Continuous Estimation of River Flows (CERF), Tech. rep., Environment Agency, Bristol, UK, 2008. a

Hanzer, F., Helfricht, K., Marke, T., and Strasser, U.: Multilevel spatiotemporal validation of snow/ice mass balance and runoff modeling in glacierized catchments, The Cryosphere, 10, 1859–1881, https://doi.org/10.5194/tc-10-1859-2016, 2016. a

Hernández-Henríquez, M. A., Sharma, A. R., and Déry, S. J.: Variability and trends in runoff in the rivers of British Columbia's Coast and Insular Mountains, Hydrol. Process., 31, 3269–3282, https://doi.org/10.1002/hyp.11257, 2017. a

Hingray, B., Schaefli, B., Mezghani, A., and Hamdi, Y.: Signature-based model calibration for hydrological prediction in mesoscale Alpine catchments, Hydrolog. Sci. J., 55, 1002–1016, https://doi.org/10.1080/02626667.2010.505572, 2010. a

Hock, R.: A distributed temperature-index ice- and snowmelt model including potential direct solar radiation, J. Glaciol., 45, 101–111, 1999. a, b

Hrachowitz, M., Fovet, O., Ruiz, L., Euser, T., Gharari, S., Nijzink, R., Freer, J., Savenije, H., and Gascuel-Odoux, C.: Process consistency in models: The importance of system signatures, expert knowledge, and process complexity, Water Resour. Res., 50, 7445–7469, 2014. a

Huss, M. and Hock, R.: Global-scale hydrological response to future glacier mass loss, Nat. Clim. Change, 8, 135–140, https://doi.org/10.1038/s41558-017-0049-x, 2018. a, b

Huss, M., Bauder, A., Funk, M., and Hock, R.: Determination of the seasonal mass balance of four Alpine glaciers since 1865, J. Geophys. Res.-Earth, 113, F01015, https://doi.org/10.1029/2007JF000803, 2008. a

Huss, M., Jouvet, G., Farinotti, D., and Bauder, A.: Future high-mountain hydrology: a new parameterization of glacier retreat, Hydrol. Earth Syst. Sci., 14, 815–829, https://doi.org/10.5194/hess-14-815-2010, 2010. a

Huss, M., Zemp, M., Joerg, P. C., and Salzmann, N.: High uncertainty in 21st century runoff projections from glacierized basins, J. Hydrol., 510, 35–48, https://doi.org/10.1016/j.jhydrol.2013.12.017, 2014. a, b, c, d, e

Immerzeel, W. W., Pellicciotti, F., and Bierkens, M. F. P.: Rising river flows throughout the twenty-first century in two Himalayan glacierized watersheds, Nat. Geosci., 6, 742–745, https://doi.org/10.1038/ngeo1896, 2013. a, b

Jackson, C., Wang, L., Pachocka, M., Mackay, J., and Bloomfield, J.: Reconstruction of multi-decadal groundwater level time-series using a lumped conceptual model, Hydrol. Process., 30, 3107–3125, https://doi.org/10.1002/hyp.10850, 2016. a

Jakob Themeßl, M., Gobiet, A., and Leuprecht, A.: Empirical-statistical downscaling and error correction of daily precipitation from regional climate models, Int. J. Climatol., 31, 1530–1544, https://doi.org/10.1002/joc.2168, 2011. a

Jansson, P., Hock, R., and Schneider, T.: The concept of glacier storage: a review, J. Hydrol., 282, 116–129, https://doi.org/10.1016/S0022-1694(03)00258-0, 2003. a

Jobst, A. M., Kingston, D. G., Cullen, N. J., and Schmid, J.: Intercomparison of different uncertainty sources in hydrological climate change projections for an alpine catchment (upper Clutha River, New Zealand), Hydrol. Earth Syst. Sci., 22, 3125–3142, https://doi.org/10.5194/hess-22-3125-2018, 2018. a, b, c, d, e, f

Jóhannesson, T., Aðlgeirsdóttir, G., Björnsson, H., and Crochet, P.: Effect of climate change on hydrology and hydro-resources in Iceland, Tech. rep., National Energy Authority, Orkugarður, 2007. a

Kelleher, C., McGlynn, B., and Wagener, T.: Characterizing and reducing equifinality by constraining a distributed catchment model with regional signatures, local observations, and process understanding, Hydrol. Earth Syst. Sci., 21, 3325–3352, https://doi.org/10.5194/hess-21-3325-2017, 2017. a

Kiesel, J., Guse, B., Pfannerstill, M., Kakouei, K., Jähnig, S. C., and Fohrer, N.: Improving hydrological model optimization for riverine species, Ecol. Indic., 80, 376–385, https://doi.org/10.1016/j.ecolind.2017.04.032, 2017. a

Kobierska, F., Jonas, T., Zappa, M., Bavay, M., Magnusson, J., and Bernasconi, S. M.: Future runoff from a partly glacierized watershed in Central Switzerland: A two-model approach, Adv. Water Resour., 55, 204–214, https://doi.org/10.1016/j.advwatres.2012.07.024, 2013. a, b

Konz, M. and Seibert, J.: On the value of glacier mass balances for hydrological model calibration, J. Hydrol., 385, 238–246, https://doi.org/10.1016/j.jhydrol.2010.02.025, 2010. a

Laghari., J. R.: Melting glaciers bring energy uncertainty, Nature, 502, 617–618, https://doi.org/10.1038/502617a, 2013. a, b

Li, H., Sheffield, J., and Wood, E. F.: Bias correction of monthly precipitation and temperature fields from Intergovernmental Panel on Climate Change AR4 models using equidistant quantile matching, J. Geophys. Res.-Atmos., 115, D10101, https://doi.org/10.1029/2009JD012882, 2010. a

Luce, C. H. and Holden, Z. A.: Declining annual streamflow distributions in the Pacific Northwest United States, 1948–2006, Geophys. Res. Lett., 36, L16401, https://doi.org/10.1029/2009GL039407, 2009. a

Lutz, a. F., Immerzeel, W. W., Shrestha, A. B., and Bierkens, M. F. P.: Consistent increase in High Asia's runoff due to increasing glacier melt and precipitation, Nat. Clim. Change, 4, 587–592, https://doi.org/10.1038/nclimate2237, 2014. a

Lutz, A. F., Immerzeel, W. W., Kraaijenbrink, P. D. A., Shrestha, A. B., and Bierkens, M. F. P.: Climate change impacts on the upper indus hydrology: Sources, shifts and extremes, PLoS ONE, 11, e0165630, https://doi.org/10.1371/journal.pone.0165630, 2016. a, b, c, d, e

Macdonald, A. M., Black, A. R., Dochartaigh, B. É. Ó., Everest, J., Darling, W. G., Flett, V., and Peach, D. W.: Using stable isotopes and continuous meltwater river monitoring to investigate the hydrology of a rapidly retreating Icelandic outlet glacier, Ann. Glaciol., 57, 1–8, https://doi.org/10.1017/aog.2016.22, 2016. a

Mackay, J., Jackson, C., and Wang, L.: A lumped conceptual model to simulate groundwater level time-series, Environ. Modell. Softw., 61, 229–245, https://doi.org/10.1016/j.envsoft.2014.06.003, 2014. a

Mackay, J., Jackson, C., Brookshaw, A., Scaife, A., Cook, J., and Ward, R.: Seasonal forecasting of groundwater levels in principal aquifers of the United Kingdom, J. Hydrol., 530, 815–828, https://doi.org/10.1016/j.jhydrol.2015.10.018, 2015. a

Mackay, J. D., Barrand, N. E., Hannah, D. M., Krause, S., Jackson, C. R., Everest, J., and Aðalgeirsdóttir, G.: Glacio-hydrological melt and run-off modelling: application of a limits of acceptability framework for model comparison and selection, The Cryosphere, 12, 2175–2210, https://doi.org/10.5194/tc-12-2175-2018, 2018. a, b, c, d, e, f, g, h, i, j

Magnússon, E., Muñoz-Cobo Belart, J., Pálsson, F., Ágústsson, H., and Crochet, P.: Geodetic mass balance record with rigorous uncertainty estimates deduced from aerial photographs and lidar data – Case study from Drangajökull ice cap, NW Iceland, The Cryosphere, 10, 159–177, https://doi.org/10.5194/tc-10-159-2016, 2016. a

Mandal, S. and Simonovic, S. P.: Quantification of uncertainty in the assessment of future streamflow under changing climate conditions, Hydrol. Process., 31, 2076–2094, https://doi.org/10.1002/hyp.11174, 2017. a

Mankin, J. S., Viviroli, D., Singh, D., Hoekstra, A. Y., and Diffenbaugh, N. S.: The potential for snow to supply human water demand in the present and future, Environ. Res. Lett., 10, 114016, https://doi.org/10.1088/1748-9326/10/11/114016, 2015. a

Mansour, M. M., Wang, L., Whiteman, M., and Hughes, A. G.: Estimation of spatially distributed groundwater potential recharge for the United Kingdom, Q. J. Eng. Geol., 51, 247–263, 2018. a

Marren, P. M.: Magnitude and frequency in proglacial rivers: a geomorphological and sedimentological perspective, Earth-Sci. Rev., 70, 203–251, https://doi.org/10.1016/j.earscirev.2004.12.002, 2005. a

Matti, B., Dahlke, H. E., Dieppois, B., Lawler, D. M., and Lyon, S. W.: Flood seasonality across Scandinavia – Evidence of a shifting hydrograph?, Hydrol. Process., 31, 4354–4370, https://doi.org/10.1002/hyp.11365, 2017. a

McDowell, J. Z. and Hess, J. J.: Accessing adaptation: Multiple stressors on livelihoods in the Bolivian highlands under a changing climate, Global Environ. Chang., 22, 342–352, https://doi.org/10.1016/j.gloenvcha.2011.11.002, 2012. a

Meresa, H. K. and Romanowicz, R. J.: The critical role of uncertainty in projections of hydrological extremes, Hydrol. Earth Syst. Sci., 21, 4245–4258, https://doi.org/10.5194/hess-21-4245-2017, 2017. a

Mora, C., Frazier, A. G., Longman, R. J., Dacks, R. S., Walton, M. M., Tong, E. J., Sanchez, J. J., Kaiser, L. R., Stender, Y. O., Anderson, J. M., Ambrosino, C. M., Fernandez-Silva, I., Giuseffi, L. M., and Giambelluca, T. W.: The projected timing of climate departure from recent variability, Nature, 502, 183–187, https://doi.org/10.1038/nature12540, 2013. a

Naiman, R. J., Latterell, J. J., Pettit, N. E., and Olden, J. D.: Flow variability and the biophysical vitality of river systems, C. R. Geosci., 340, 629–643, https://doi.org/10.1016/j.crte.2008.01.002, 2008. a

Nawri, N., Pálmason, B., Petersen, G. N., Björnsson, H., and Þorsteinsson, S.: The ICRA atmospheric reanalysis project for Iceland, Tech. rep., Icelandic Meteorological Office, Reykjavík, Iceland, 2017. a, b

Nolin, A. W., Phillippe, J., Jefferson, A., and Lewis, S. L.: Present-day and future contributions of glacier runoff to summertime flows in a Pacific Northwest watershed: Implications for water resources, Water Resour. Res., 46, W12509, https://doi.org/10.1029/2009WR008968, 2010. a

Pellicciotti, F., Brock, B., Strasser, U., Burlando, P., Funk, M., and Corripio, J.: An enhanced temperature-index glacier melt model including the shortwave radiation balance: development and testing for Haut Glacier d'Arolla, Switzerland, 51, 573–587, 2005. a, b

Phillips, E., Finlayson, A., Bradwell, T., Everest, J., and Jones, L.: Structural evolution triggers a dynamic reduction in active glacier length during rapid retreat: Evidence from Falljökull, SE Iceland, J. Geophys. Res.-Earth, 119, 2194–2208, https://doi.org/10.1002/2014JF003165, 2014. a

Pianosi, F., Beven, K., Freer, J., Hall, J. W., Rougier, J., Stephenson, D. B., and Wagener, T.: Sensitivity analysis of environmental models: A systematic review with practical workflow, Environ. Modell. Softw., 79, 214–232, https://doi.org/10.1016/j.envsoft.2016.02.008, 2016. a

Ponce, V. M.: Engineering hydrology: Principles and practices, Prentice-Hall, Englewood Cliffs, New Jersey, 1989. a

Pool, S., Vis, M. J. P., Knight, R. R., and Seibert, J.: Streamflow characteristics from modeled runoff time series – importance of calibration criteria selection, Hydrol. Earth Syst. Sci., 21, 5443–5457, https://doi.org/10.5194/hess-21-5443-2017, 2017. a

Ragettli, S., Pellicciotti, F., Bordoy, R., and Immerzeel, W. W.: Sources of uncertainty in modeling the glaciohydrological response of a Karakoram watershed to climate change, Water Resour. Res., 49, 6048–6066, https://doi.org/10.1002/wrcr.20450, 2013. a, b

Ragettli, S., Immerzeel, W. W., and Pellicciotti, F.: Contrasting climate change impact on river flows from high-altitude catchments in the Himalayan and Andes Mountains., P. Natl. Acad. Sci. USA, 113, 9222–9227, https://doi.org/10.1073/pnas.1606526113, 2016. a, b

Riggs, G. and Hall, D.: MODIS Snow Products Collection 6 User Guide, Tech. rep., available at: https://nsidc.org/sites/nsidc.org/files/files/MODIS-snow-user-guide-C6.pdf (last access: 6 October 2016), 2015. a

Samaniego, L., Kumar, R., Breuer, L., Chamorro, A., Flörke, M., Pechlivanidis, I. G., Schäfer, D., Shah, H., Vetter, T., Wortmann, M., and Zeng, X.: Propagation of forcing and model uncertainties on to hydrological drought characteristics in a multi-model century-long experiment in large river basins, Climatic Change, 141, 435–449, https://doi.org/10.1007/s10584-016-1778-y, 2017. a

Sanford, T., Frumhoff, P. C., Luers, A., and Gulledge, J.: The climate policy narrative for a dangerously warming world, Nat. Clim. Change, 4, 164–166, https://doi.org/10.1038/nclimate2148, 2014. a

Sawicz, K. A., Kelleher, C., Wagener, T., Troch, P., Sivapalan, M., and Carrillo, G.: Characterizing hydrologic change through catchment classification, Hydrol. Earth Syst. Sci., 18, 273–285, https://doi.org/10.5194/hess-18-273-2014, 2014. a

Schaefli, B.: Snow hydrology signatures for model identification within a limits-of-acceptability approach, Hydrol. Process., 30, 4019–4035, https://doi.org/10.1002/hyp.10972, 2016. a, b

Schaefli, B. and Huss, M.: Integrating point glacier mass balance observations into hydrologic model identification, Hydrol. Earth Syst. Sci., 15, 1227–1241, https://doi.org/10.5194/hess-15-1227-2011, 2011. a

Seibert, J., Vis, M. J. P., Kohn, I., Weiler, M., and Stahl, K.: Technical note: Representing glacier geometry changes in a semi-distributed hydrological model, Hydrol. Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018, 2018. a, b

Shafii, M. and Tolson, B. A.: Optimizing hydrological consistency by incorporating hydrological signatures intomodel calibration objectives, Water Resour. Res., 51, 3796–3814, 2015. a, b, c

Shea, J. M. and Immerzeel, W. W.: An assessment of basin-scale glaciological and hydrological sensitivities in the Hindu Kush-Himalaya, Ann. Glaciol., 57, 308–318, https://doi.org/10.3189/2016AoG71A073, 2016. a

Shea, J. M. and Moore, R. D.: Prediction of spatially distributed regional-scale fields of air temperature and vapor pressure over mountain glaciers, J. Geophys. Res.-Atmos., 115, D23107, https://doi.org/10.1029/2010JD014351, 2010. a

Singh, S., Kumar, R., Bhardwaj, A., Sam, L., Shekhar, M., Singh, A., Kumar, R., and Gupta, A.: Changing climate and glacio-hydrology in Indian Himalayan Region: A review, Wiley Interdisciplinary Reviews: Climate Change, 7, 393–410, https://doi.org/10.1002/wcc.393, 2016. a

Sorensen, J. P. R., Finch, J. W., Ireson, A. M., and Jackson, C. R.: Comparison of varied complexity models simulating recharge at the field scale, Hydrol. Process., 28, 2091–2102, https://doi.org/10.1002/hyp.9752, 2014. a

Stewart, I. T., Ficklin, D. L., Carrillo, C. A., and McIntosh, R.: 21st century increases in the likelihood of extreme hydrologic conditions for the mountainous basins of the Southwestern United States, J. Hydrol., 529, 340–353, https://doi.org/10.1016/j.jhydrol.2015.07.043, 2015. a, b

Stoffel, M., Wyżga, B., and Marston, R. A.: Floods in mountain environments: A synthesis, Geomorphology, 272, 1–9, https://doi.org/10.1016/j.geomorph.2016.07.008, 2016. a

Tabachnick, B. G. and Fidell, L. S.: Using Multivariate Statistics Sixth Edition, Pearson Education Limited, Essex, United Kingdom, sixth edn., 2014. a

Taylor, K. E., Stouffer, R. J., and Meehl, G. a.: An Overview of CMIP5 and the Experiment Design, B. Am. Meteorol. Soc., 93, 485–498, https://doi.org/10.1175/BAMS-D-11-00094.1, 2012. a

Teutschbein, C., Grabs, T., Karlsen, R. H., Laudon, H., and Bishop, K.: Hydrological response to changing climate conditions: Spatial streamflow variability in the boreal region, Water Resour. Res., 51, 9425–9446, https://doi.org/10.1002/2015WR017337, 2015. a

Thorsteinsson, T. and Björnsson, H.: Climate Change and Energy Systems: Impacts, Risks and Adaptation in the Nordic and Baltic countries, Tech. rep., Nordic Council of Ministers, Copenhagen, 2012. a

Vaughan, D., Comiso, J., Allison, I., Carrasco, J., Kaser, G., Kwok, R., Mote, P., Murray, T., Paul, F., Ren, J., Rignot, E., Solomina, O., Steffen, K., and Zhang, T.: Observations: Cryosphere, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T., Qin, D., Plattner, G.-K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 358–359, 2013. a

Vetter, T., Huang, S., Aich, V., Yang, T., Wang, X., Krysanova, V., and Hattermann, F.: Multi-model climate impact assessment and intercomparison for three large-scale river basins on three continents, Earth Syst. Dynam., 6, 17–43, https://doi.org/10.5194/esd-6-17-2015, 2015. a, b, c, d

Vetter, T., Reinhardt, J., Flörke, M., van Griensven, A., Hattermann, F., Huang, S., Koch, H., Pechlivanidis, I. G., Plötner, S., Seidou, O., Su, B., Vervoort, R. W., and Krysanova, V.: Evaluation of sources of uncertainty in projected hydrological changes under climate change in 12 large-scale river basins, Climatic Change, 141, 419–433, https://doi.org/10.1007/s10584-016-1794-y, 2017. a, b, c

Viviroli, D. and Weingartner, R.: The hydrological significance of mountains: from regional to global scale, Hydrol. Earth Syst. Sci., 8, 1017–1030, https://doi.org/10.5194/hess-8-1017-2004, 2004. a

Viviroli, D., Dürr, H. H., Messerli, B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43, 1–13, https://doi.org/10.1029/2006WR005653, 2007. a

von Storch, H. and Zwiers, F. W.: Statistical Analysis in Climate Research, Cambridge University Press, Cambridge, United Kingdom, 1999. a

Wijngaard, R. R., Lutz, A. F., Nepal, S., Khanal, S., Pradhananga, S., Shrestha, A. B., and Immerzeel, W. W.: Future changes in hydro-climatic extremes in the Upper Indus, Ganges, and Brahmaputra River basins, PLOS ONE, 12, e0190224, https://doi.org/10.1371/journal.pone.0190224, 2017. a, b

Willis, I.: 168: Hydrology of Glacierized Basins, in: Encyclopedia of Hydrological Sciences: Part 14. Snow and Glacier Hydrology, edited by: Anderson, M. G. and McDonnell, J. J., John Wiley & Sons, Ltd, Chichester, UK, 2005. a

Yadav, M., Wagener, T., and Gupta, H.: Regionalization of constraints on expected watershed response behavior for improved predictions in ungauged basins, Adv. Water Resour., 30, 1756–1774, https://doi.org/10.1016/j.advwatres.2007.01.005, 2007.  a, b

Yilmaz, K. K., Gupta, H. V., and Wagener, T.: A process-based diagnostic approach to model evaluation: Application to the NWS distributed hydrologic model, Water Resour. Res., 44, W09417, https://doi.org/10.1029/2007WR006716, 2008. a

Yuan, F., Zhao, C., Jiang, Y., Ren, L., Shan, H., Zhang, L., Zhu, Y., Chen, T., Jiang, S., Yang, X., and Shen, H.: Evaluation on uncertainty sources in projecting hydrological changes over the Xijiang River basin in South China, J. Hydrol., 554, 434–450, https://doi.org/10.1016/j.jhydrol.2017.08.034, 2017. a, b

Zappa, M. and Kan, C.: Extreme heat and runoff extremes in the Swiss Alps, Nat. Hazards Earth Syst. Sci., 7, 375–389, https://doi.org/10.5194/nhess-7-375-2007, 2007. a

Zemp, M., Frey, H., Gärtner-Roer, I., Nussbaumer, S. U., Hoelzle, M., Paul, F., Haeberli, W., Denzinger, F., Ahlstrøm, A. P., Anderson, B., Bajracharya, S., Baroni, C., Braun, L. N., Càceres, B. E., Casassa, G., Cobos, G., Dàvila, L. R., Delgado Granados, H., Demuth, M. N., Espizua, L., Fischer, A., Fujita, K., Gadek, B., Ghazanfar, A., Hagen, J. O., Holmlund, P., Karimi, N., Li, Z., Pelto, M., Pitte, P., Popovnin, V. V., Portocarrero, C. A., Prinz, R., Sangewar, C. V., Severskiy, I., Sigurdsson, O., Soruco, A., Usubaliev, R., and Vincent, C.: Historically unprecedented global glacier decline in the early 21st century, J. Glaciol., 61, 745–762, https://doi.org/10.3189/2015JoG15J017, 2015. a