The compensating effect of glaciers: Characterizing the relation between interannual streamflow variability and glacier cover

Meltwater from glaciers is not only a stable source of water but also affects downstream streamflow dynamics. One of these dynamics is the interannual variability of streamflow. Glaciers can moderate streamflow variability because the runoff in the glacierized part, driven by temperature, correlates negatively with the runoff in the non‐glacierized part of a catchment, driven by precipitation, thereby counterbalancing each other. This is also called the glacier compensation effect (GCE), and the effect is assumed to depend on relative glacier cover. Previous studies found a convex relationship between streamflow variability and glacier cover of different glacierized catchments, with lowest streamflow variability at a certain optimum glacier cover. In this study, we aim to revisit these previously found curves to find out if a universal relationship between interannual streamflow variability and glacier cover exists, which could potentially be used in a space‐for‐time substitution analysis. Moreover, we test the hypothesis that the dominant climate drivers (here precipitation and temperature) switch around the suggested optimum of the curve. First, a set of virtual nested catchments, with the same absolute glacier area but varying non‐glacierized area, were modelled to isolate the effect of glacier cover on streamflow variability. The modelled relationship was then compared with a multicatchment data set of gauged glacierized catchments in the European Alps. In the third step, changes of the GCE curve over time were analysed. Model results showed a convex relationship and the optimum in the simulated curve aligned with a switch in the dominant climate driver. However, the multicatchment data and the time change analyses did not suggest the existence of a universal convex relationship. Overall, we conclude that GCE is complex due to entangled controls and changes over time in glacierized catchments. Therefore, care should be taken to use a GCE curve for estimating and/or predicting interannual streamflow variability in glacierized catchments.

In mountain regions, glaciers can decrease the streamflow variability in partly glacierized catchments because the glacierized part of the catchment has an opposite runoff regime compared with the non-glacierized part of the catchment, thereby counterbalancing each other (Braithwaite & Olesen, 1988;Meier & Tangborn, 1961;Rothlisberger & Lang, 1987). The non-glacierized part is characterized by a rainfall-runoff regime and the glacierized part by a melt-dominated regime. This means that, in a catchment where both of these runoff regimes are present, during warm and dry periods, the melt from the glacierized part could compensate for the lack of precipitation in the non-glacierized part (e.g., Koboltschnig, Schöner, Holzmann, & Zappa, 2009;Zappa & Kan, 2007). On the other hand, during cold and wet periods, runoff from the non-glacierized part could compensate for the lack of melt from the glacierized part (e.g., Hopkinson & Young, 1998). This process is called the glacier compensation effect (GCE, e.g., Fountain & Tangborn, 1985;Rothlisberger & Lang, 1987).
How much the two different hydrological regimes in a glacierized catchment can counterbalance each other is considered to depend on the relative area of the glacierized part and the non-glacierized part (e.g., Fountain & Tangborn, 1985), that is, the catchment's glacier cover fraction ( g). If the non-glacierized part covers a large part of the catchment, precipitation variability will mainly control streamflow variability. If on the other hand the glacier covers a large part of the catchment, streamflow variability is dominated by temperature variability (e.g., Casassa, López, Pouyaud, & Escobar, 2009;Collins, 2006a;Rothlisberger & Lang, 1987). Several empirical studies looked at the relation between g and streamflow variability, expressed as the coefficient of variation (CV Q ) for different samples of catchments (Braithwaite & Olesen, 1988;Chen & Ohmura, 1990;Collins, 2006b;Fleming & Clarke, 2005;Fountain & Tangborn, 1985;Krimmel & Tangborn, 1974;Meier & Tangborn, 1961;Moore, 1992), and some of them found a convex nonlinear relationship with higher CV Q at low and high g (Chen & Ohmura, 1990;Fountain & Tangborn, 1985). This convex relationship (see conceptual curve in Figure 1) was found to have an optimum (minimum CV Q , maximum compensation, GCE opt ) at 36% g for a U.S. data set (Fountain & Tangborn, 1985) and at 39-44% g for a data set in the European Alps (Chen & Ohmura, 1990). Most of the other studies only looked at a small sample of glacierized catchments and were therefore only able to confirm parts of the relationship. Moreover, Moore (1992) and Collins (2006b) discuss that finding a relation between g and CV Q can be complicated because other factors than just g, for example, elevation range, glacier hypsometry, climate characteristics, and snow on the glacier, influence the interplay between runoff from the glacierized part (Q gl ) and the non-glacierized part (Q nongl ).
Despite these few empirical studies that suggest a nonlinear convex relationship between g and CV Q , later research has not built on these studies to establish a universal relationship for the influence of glaciers on streamflow variability. However, some studies used the relationships found to estimate and/or compare CV Q (e.g., Hopkinson & Young, 1998). Also, research into the processes' relative importance underlying the hypothesized curve has been limited. This is in stark contrast to the many studies that mention the compensation effect of glaciers as general motivation for why it is important to study the role of glaciers in hydrology. If a universal relationship could be found, it would help to better understand spatial differences in flow reliability. This would have a number of advantages for water resources planning such as decisions on where to install small hydropower plants or intakes for water supply (Gaudard et al., 2014;Schaefli et al., 2019).
Also, if a universal relationship between g and CV Q existed, it could be used to estimate how CV Q of a specific catchment would change over time due to deglaciation (following the GCE curve from right to left), which is crucial information for water management and planning. The intraregional comparison of catchments' CV Q and g in previous studies might thus be interpreted as a space-for-time substitution (Singh, Wagener, Werkhoven, Mann, & Crane, 2011). This F I G U R E 1 Conceptual understanding of GCE curve, after Chen and Ohmura (1990). CV Q is the coefficient of variation and represents interannual streamflow variability, g is the relative glacier cover, and GCE opt indicates the optimum g where CV Q is lowest. r P-Q and r T-Q are the correlations of streamflow with precipitation and temperature, respectively translation, going from a relationship based on different catchments to a relationship for one catchment over time, has recently been applied for the Budyko curve (Carmona, Sivapalan, Yaeger, & Poveda, 2014;Sivapalan, Yaeger, Harman, Xu, & Troch, 2011). Hock, Jansson, and Braun (2005) already interpreted the GCE curve as a space-fortime substitution. They hypothesized that under global warming, the interannual variability of streamflow will first decrease (or, depending on initial g, increase) and increase later on. However, we argue that before any time interpretations of the curve are made, a more systematic analysis of the compensation effect and the relationship between g and CV Q is needed.
The GCE curve, as sketched in Figure 1, could be interpreted in three ways: 1 several nested catchments with the same absolute glacier area but different g due to larger/smaller non-glacierized area (with outlet further downstream/upstream); 2 different catchments within a specific region (e.g., Chen & Ohmura, 1990); and 3 one catchment with changing g over time due to deglaciation (e.g., hypothesized by Hock et al., 2005).
A data set of the first type has not been used to define or corroborate a GCE curve yet, most likely because there is no appropriate streamflow data set available from several nested undisturbed catchments with the same absolute glacier area covering a g range large enough to fit or test a relationship. Neither has an examination of the climate drivers of streamflow variability in relation to the GCE curve been carried out. We hypothesize that precipitation is the dominant driver of streamflow variability below GCE opt , and temperature is the dominant driver above GCE opt . So, catchments to the right of the optimum would have a higher correlation of streamflow with temperature (r T-Q ) than with precipitation (r P-Q ) and vice versa for catchments to the left of the optimum ( Figure 1). If data can confirm this hypothesis, then the correlation of streamflow with temperature and precipitation might be an indicator of the catchments' location relative to the optimum.
In this study, we test the different interpretations of the GCE curve and its climate controls with the aim of challenging the existence of a universal (convex) relationship between streamflow variability and glacier cover. This will give a better insight in whether the relationship can be used for planning and predictions of future water reliability in glacierized catchments under climate change. By a universal relationship we mean that, for a certain region, a convex relationship can be fitted through some data points and that other catchments, as well as a particular catchment during deglaciation, in this region follow(s) this curve.
The study is organized as follows: First, we perform a model experiment with virtual nested catchments (Interpretation 1 as mentioned above) to isolate the effect of g on CV Q , to look at the streamflow components Q gl and Q nongl , and to investigate the climate drivers. Next, the modelled relationship is compared with observation data of catchments in the Swiss and Austrian Alps (Interpretation 2) to examine whether such a GCE relationship can be found for the region of the Alps and to confirm or reject, with more recent data, the previously found relationship for the Alps (from Chen & Ohmura, 1990). In addition to the g-CV Q relationship, correlations with climate drivers are analysed. Last, we explore the time aspect of the relationship (Interpretation 3). We evaluate the g-CV Q relationship in different periods and analyse changes of CV Q over time for a selection of long streamflow records.

| Model experiment
As a first step in analysing the relationship between CV Q and g, a model experiment was carried out (e.g., Rees & Collins, 2006). A virtual nested catchment model framework enabled the analysis of the compensating effect in a controlled way. Nested catchments form a natural laboratory to analyse the influence of differences in g while most other factors are similar. Specifically, we aim to test if realizing this nested catchment approach by hydrological modelling results in a convex relationship as hypothesized.
The setup of our virtual nested model catchments is based on a real catchment, using observed glacier geometries and area elevation distributions ( Figure 2). The catchment taken as starting point is the Schwarze Lütschine, located on the north side of the Swiss Alps. The glacier outline for the year 1973 and a digital elevation model were used to define the glacierized area and elevation zones with a 100-m interval. g in the original catchment, comprising subcatchments 1-4, is 20.1% (36.1 km 2 ). To design model catchments with 0% < g < 100%, we decreased (increased) the non-glacierized area of the original catchment to increase (decrease) g. In this way, we vary g, but keep glacier geometry constant and climate input comparable between the different nested model catchments.
Increasing g was done by removing the subcatchments (Table 1, Figure 2), one after the other, thereby reducing the original catchment area. Adding additional non-glacierized area to the original catchment required a virtual approach because the outlet in the real setting cannot be moved further downstream (another glacierized tributary joining and the river ending in Lake Brienz). Therefore, a clone of subcatchment 1 was added as many times as was needed to create a model catchment with g of 15%, 10%, 5%, 3%, and 1% (Table 1). For the lowest g (0%) model catchment, we took the virtual catchment g 1 (Table 1), and changed the glacier area in this catchment into nonglacierized area. For g 100 , the catchment area corresponded to the glacier area.
The HBV-light software (Seibert & Vis, 2012;Seibert, Vis, Kohn, Weiler, & Stahl, 2018) was used to model the streamflow of these nested glacierized catchments. HBV-light is a semidistributed model based on hydrological response units defined by elevation zones, aspect classes, and glacierized/non-glacierized areas. The model has different routines: snow and glacier, soil, response, and routing. The melt in the snow and glacier routine is calculated with a degree-day method and a higher degree-day factor is used for glacier ice compared with snow. The model simulates discharge in daily timesteps (Q) and also calculates the streamflow originating from the glacierized (Q gl ) and the non-glacierized part (Q nongl ) of the catchments, so that Q = Q gl + Q nongl whereby Q gl includes streamflow from rain, snowmelt, and ice melt, and Q nongl includes streamflow from rain and snow melt. For a detailed model description, we refer to Seibert and Vis (2012)

| Multi-catchment data set
The relationship between CV Q and g from the modelled nested catchments was compared with the relationship obtained from streamflow observations of 39 gauged catchments in the Swiss and Austrian Alps (Table A1, Figure 3). All catchments are situated at approximately the same latitude, but they are distributed from west to east over the Alps. Their g range from 0% to 73.5%, based on outlines for the year 2010 and 2012 for Switzerland and Austria, respectively (Table 2). For reference, the data set also contains four Swiss non-glacierized catchments.
Most catchment areas are between 10 and 100 km 2 with some larger ones between 200 and 517 km 2 . The mean elevation is around 2,000-2,500 m a.s.l. for the glacierized and around 1,000-1,500 m a.s.l. for the non-glacierized catchments. The highest catchment (3,122 m a.s.l. mean elevation) is the Vernagtbach catchment in Austria, which also has the highest glacier cover and is one of the smallest catchments. In a few catchments, the elevation range is more than 2,500 m (see Table A1). Glacier area ranges from smaller than 1 km 2 to larger than 100 km 2 .
Mean annual precipitation sums in the catchments range from 940 mm to more than 2,300 mm. Precipitation is highest in summer for all catchments, with a second smaller peak around November.
Monthly mean temperatures are in general above 0 C between May and October for the glacierized catchments and only below 0 C in winter (December-February) for the non-glacierized catchments.
Mean annual temperatures are higher for the non-glacierized catchments because of their lower elevation. Daily P and T data are taken from interpolated observations: the SPARTACUS gridded product (1 × 1 km, Hiebl & Frei, 2016 for Austria and the RhiresD (~1×1 km, precipitation) and TabsD   . These years were excluded from the analysis for these catchments. All streamflow records are assumed to be undisturbed.
However, two catchments are flagged by the FOEN as influenced (e.g., by subdaily hydropeaking; see Table A1). As no typical hydropeaking signal or other disturbance signal was evident in the daily streamflow time series, we assumed that the influence is negligible with the monthly and annual time resolutions used in this study.

| Analysis of streamflow variability
The coefficient of variation (CV) was used as metric for interannual streamflow variability (Chen & Ohmura, 1990;Fountain & Tangborn, 1985). CV Q was calculated for yearly values of annual streamflow (Q ann ) and August streamflow (Q aug ). August was selected because it can be assumed to be the month with highest ice melt because seasonal snow has melted from the glacier surface (Jost, Moore, Menounos, & Wheate, 2012;Lang, 1973;Stahl & Moore, 2006), although some studies also found September to be the month with highest glacier melt contribution (e.g., Frans, Istanbulluoglu, Lettenmaier, Fountain, & Riedel, 2018).
CV Q was plotted against g to test for the presence of a GCE curve ( Figure 1). Leap days were removed before aggregating daily streamflow series. Only streamflow sums from years or months without any data gaps were considered. CV Q was calculated for the period 1976-2015 both for the model experiment and the multicatchment data set.
To analyse the potential controls on the GCE curve (see Introduction) and GCE opt , we examined how streamflow variability relates to precipitation and temperature variability, by calculating the Spearman correlation coefficient between streamflow and the two climate variables (r P-Q , r T-Q ). Correlation coefficients were calculated for all catchments from the model experiment and the multi-catchment data set for the same period as CV Q (1976-2015).

| Time stability of GCE curve
For some of the Swiss catchments in the multi-catchment data set, long time series are available (>40 years). For a selection of those, CV Q was calculated for different subperiods (see Table A1). The definition of the subperiods was based on the reference years of available g data. We used the glacier outline data sets of 1935, 1973, and 2003 (Table 2). The corresponding 20-year subperiods to derive CV Q were 1932-1951, 1965-1984, and 1996-2015. Plotting, for each catchment, the CV Q of the subperiods against the g of the subperiods shows whether changes in g over time lead to expected changes in CV Q (following the curve, decreasing, or increasing CV Q ).
Furthermore, we also tested how stable the relationship is in time by calculating CV Q for two 20-year periods (1965-1984 and 1996-2015 Chen and Ohmura (1990), and compared for the two periods. Also here, the three to five gap years in three of the catchments in the different analyses periods were excluded from the analyses.

| Model experiment
The CV Q values of the virtual nested catchments show a convex relationship with g (Figure 4a,b). The GCE curve optimum lies around g = 15% for Q ann and around g = 10% for Q aug . CV Q is generally higher for Q aug than for Q ann , especially at high and low g. The CV Q of the g 100 catchment shows some deviation from a smooth curve. The correlations of streamflow with precipitation and temperature cross at g 15-20% for Q ann and at g 10-15% for Q aug , which corresponds with the optimum in the GCE curve (Figure 4c,d). As expected, streamflow of catchments with high g correlates more strongly with temperature than with precipitation, for the low g catchments it is opposite. With  In catchments with g lower than the optimum ( g 10 and g 3 , Figure 5, left column), the Q nongl anomalies are larger than the Q gl anomalies. A negative anomaly in Q nongl , for example, due to a precipitation deficit, cannot be offset by a positive anomaly in Q gl because Q gl is not large enough. On the other hand, for catchments with g above the optimum ( g 37 and g 20 , Figure 5, right column), the Q gl anomalies are much larger than those of Q nongl and, therefore, they can also not counterbalance each other. At g near the GCE curve optimum ( g 15 and g 10 , Figure 5, middle column), the Q gl and Q nongl anomalies are similar. Besides being of comparable magnitude, the direction of the anomalies should also be opposite to counterbalance each other. In most years, this is the case (e.g., 1977, 1978), but sometimes the anomalies have the same sign (e.g., 1996 for Q ann ). Total streamflow ( Figure 5) shows the effect of the interplay of the two components. For example, at low g, Q ann and Q aug are mostly determined by anomalies in Q nongl , but in years when Q gl is more substantial and opposite to Q nongl , total streamflow is closer to the mean streamflow. Near GCE opt , total streamflow shows large deviations when the anomalies have the same sign and close to average values when the anomalies are opposite. Figure 6 shows the same type of analysis as in Figure 4a-d, but now for the observed streamflow records from the multicatchment data set instead of the modelled nested catchments. In general, the CV Q of observed streamflow is higher than the modelled CV Q (dashed line). The CV Q of the multicatchment data set shows a scatter, and it is difficult to find a clear curve with an optimum.

| Multi-catchment empirical analysis
For Q ann , CV Q decreases with increasing g between 0% and 15%.
This part of the graph is, however, dominated by the four nonglacierized catchments showing high CV Q (Figure 6a,b). There are not many catchments with a high g. The Vernagtbach catchment (with highest g, 73.5%) is an important data point to indicate a potential increase in CV Q going from moderate to high g, especially for Q ann . However, its data are uncertain for Q ann , due to the infilling of mean winter streamflow values. For both Q ann and Q aug , F I G U R E 4 GCE curve, its drivers, and the streamflow components obtained from the simulations of the nested model catchments. Left column is Q ann , right column Q aug . (a) and (b) show the relation between g and CV Q , (c) and (d) the correlations r T-Q (red) and r P-Q (blue), (e) and (f) the ratio of the standard deviation of the streamflow components Q gl and Q nongl , and (g) and (h) show the ratio of total Q gl and Q nongl in boxplots (data of all years). The dashed lines in a, b, c, and d connect the results from the different catchments the Vernagtbach catchment shows a much higher CV Q compared with the model experiment. For Q aug , the decrease and subsequent increase in CV Q from non-to high-glacierized catchments is more evident than for Q ann , but the scatter near the modelled optimum is large. The correlations of streamflow with precipitation and temperature, r P-Q and r T-Q , are more comparable to the model results. have approximately the same value rather than one crossover point. This range is narrower for Q aug than for Q ann . The few catchments with g < 1% show a larger range of correlation values than that were modelled for the two nested catchments with lowest g, g 0 , and g 1 .

| Changes of streamflow variability over time
Another  Figure 7) shows a reversed pattern of the GCE curve because instead of decreasing and then increasing CV Q when passing the optimum g, the data shows increasing and then decreasing CV Q . Most of the catchments show a different trend in CV Q for Q ann compared with Q aug .
Using a larger subsample of the multicatchment data-set and fitting a second-order polynomial through the data points, as done by Chen and Ohmura (1990), result in very different curves for different time periods (Figure 8, left). The curve for the later period, 1996-2015, agrees most with the curve from Chen and Ohmura (1990), although in their study data from years until 1985 were used. Moreover, especially the earlier period, 1965-1984, shows that a second-order polynomial might not be the best function to describe the relation between CV Q and g. Also in the model experiment, the fitted curves differ between the two periods ( Figure 8b). The later period, 1996-2015, shows lower CV Q for the moderate-to high-glacierized catchments compared with the earlier period, 1976-2015.

| Is there a universal relationship between glacier cover and interannual streamflow variability?
The model experiment confirmed the hypothesis of a convex relationship between g and CV Q (Figure 4). In a controlled nested catchment setting, glacier cover thus influences streamflow variability. The left part of the curve (where rainfall-runoff regimes are dominant) was steeper than the right part (where melt regimes are dominant), more so for August streamflow. This suggests that precipitation and temperature may influence the streamflow variability differently and is broadly consistent with the asymmetric behaviour found for correlations of streamflow with precipitation and temperature. Results of this model experiment, however, do not support the simple quadratic shape of the GCE curve suggested by Chen and Ohmura (1990) for the Alps. The multicatchment data set had overall higher CV Q than the modelled virtual catchments; the modelled GCE curve appears to represent a lower envelope of the multicatchment CV Q . This is possibly due to the general tendency of underestimation of streamflow variability, including extremes, by a model (e.g., Whitfield, Wang, & Cannon, 2003) and due to our idealized nested catchment model framework in which, for example, parameters do not change with F I G U R E 7 CV Q of Q ann (left) and Q aug (right) for three 20-year time windows around glacier cover inventory dates (different symbols) for eight Swiss catchments (different colours) with long observed streamflow time series. The dashed grey curve in the background is from the model experiment as shown in Figure 4a and b catchment size, glaciers are modelled statically, and other simplifications. The scatter was high and the multicatchment data set did only roughly show a GCE curve ( Figure 6). Decreasing streamflow variability with increasing g in the low-glacierized part was evident for Q ann and Q aug . At the GCE opt of the model experiment for Q aug , the multicatchment data set showed a large range of CV Q. An increase in CV Q to the right of an assumed optimum was less evident and depended strongly on one data point for the catchment with the highest glacier cover (73.5% Vernagtbach), which has infilled streamflow data in winter due to discontinuous streamflow monitoring.
In the multicatchment data set and the model experiment, the correlations with climate drivers showed a crossover. From low-to highglacierized catchments, the correlation of streamflow with precipitation decreased, whereas the correlation with temperature increased.
However, only in the model experiments, due to the clear GCE curve, the crossover of correlations could be clearly related to the optimum in the GCE curve. In the multicatchment data set, a similar correlation in r P-Q might have resulted in different CV Q because of variations in, for example, mean streamflow or magnitude of streamflow anomalies due to dissimilarity in storage characteristics. Also, glacier characteristics can vary with respect to, for example, their water retention capacity (e.g., Jansson et al., 2003;Hock et al., 2005) or aspect (e.g., Hock, 2003Hock, , 2005, which influences the balance of Q gl and Q nongl on shorter or longer timescales and could cause the scatter in the CV Q and r P-Q and r T-Q for catchments with similar g. The correlations might therefore only be a very rough indicator of a catchment's location relative to the GCE optimum. We showed that a universal relationship between CV Q and g in the Alps that could be used for quantitative planning does not exist. A space-for-time interpretation of the curves postulated by Chen and Ohmura (1990) and Fountain and Tangborn (1985) is not advised based on our analyses for two reasons: (a) In our intraregional analysis of streamflow observations, we could not find a distinct relation between CV Q and g following a theoretical GCE curve, and (b) changes over time of CV Q do not follow one common curve and also the relation itself is not stable over time. Possible reasons for the differences between our study and Chen and Ohmura (1990) could be the different time periods analysed. We chose a common period for all catchments , whereas Chen and Ohmura (1990) took different lengths of time series as available for the different catchments. If a different time period had caused the different results, this would mean that the relationship is not stable over time. However, redoing the exact same analyses as Chen and Ohmura (1990) with the same catchments for the same period, we found that some catchments were apparently not included in the fit of the curve given by Chen and Ohmura (1990) although being mentioned in their paper ( Figure 9). We do not know the reason for leaving those out, but it might have influenced the results.
Relative glacier cover might not be the best indicator of the balance between the two runoff regimes in a partly glacierized F I G U R E 8 Relation between CV Q and g for different 20-year time periods for Q ann . Left shows the multicatchment data set, right the nested model catchments. A second-order polynomial was fitted through the points. Only catchments that have streamflow data for both periods were used. The curve from Chen and Ohmura (1990) is plotted for comparison in the left panel. In the right panel, the GCE curve from the model experiment (40-year data) is shown in grey F I G U R E 9 Reconstruction of the GCE curve by Chen and Ohmura (1990). Overlapping catchments in our study with Chen and Ohmura (1990) were plotted for the same time period. CV Q was calculated for annual streamflow, Hydrological years (HY), summer months June, July and August (JJA), and summer half year (MJJAS). CV Q is expressed as fraction, not percentage. Grey lines indicate catchments that were listed in the data section of Chen and Ohmura (1990) but not included in the figure of the quadratic fit catchment. When glaciers melt due to warming, they may thin, the snow and firn reservoirs on the glaciers may reduce, and the accumulation ablation ratio can change (e.g., Dyurgerov, 2003;Paul, Kääb, & Haeberli, 2007), all before notable decreases in g occur. These changes, however, could influence the amount of Q gl . These processes might play a role when comparing different catchments with similar g and when analysing changes over time for a particular catchment. Thus, also the balance of the glacier with climate can influence the GCE. A more balanced glacier has a different Q gl component than a glacier that is out of balance (e.g., Pritchard, 2019). European glaciers experienced close to balanced conditions around 1960-1980(e.g., Vincent et al., 2017. The period analysed in Chen and Ohmura (1990) mostly covers this rather stable period, whereas the more recent period used in this study  is characterized by negative mass balances and retreating glaciers (Huss, Bauder, Funk, & Hock, 2008).

| Limitations and other controls on GCE
The multicatchment analyses as well as the model experiment revealed a number of limitations for the analysis of GCE. For example, undisturbed streamflow observations from high-glacierized catchments are rare. Such data are, however, important for analysing the right part of a potential GCE curve, particularly the steepness of this part of the curve. The Vernagtbach catchment (g = 73.5%) proved its importance in our analysis, as well as in the study of Chen and Ohmura (1990), but it might be a unique catchment due to its small elevation range and a glacier that has retreated to a half-circle shape (with different aspects) and that has known surging behaviour (Hoinkes, 1969). Furthermore, the lack of observations of glacier cover at different historical times limits the combined analyses with long streamflow records.
In the study area, observed decreases in g over the time periods covered by streamflow and glacier area observations are rather small (in the order of 5-10%). That makes it difficult to attribute changes in observed streamflow variability to changes in glacier cover. The few catchments with long time series only represent a small part of the theoretical curve. Glaciers adjust their size and elevation distribution to be in balance with climate, that is, in a warming climate, they retreat to higher elevations. These adjustments affect Q gl and Q nongl , and the GCE may depend on glacier mass balance and on the phase of peak glacier melt water (e.g., Huss & Hock, 2018;Immerzeel, Pellicciotti, & Bierkens, 2013). Another modelling experiment would be needed to analyse how the streamflow variability of (different) high-glacierized catchments develops over time when the glacier is retreating and even disappears (and the climate is changing)(e.g., Hagg, Braun, Kuhn, & Nesgaard, 2006;Nolin, Phillippe, Jefferson, & Lewis, 2010 Collins (2006aCollins ( , 2006b and Moore (1992) discuss the importance of the interannual variability of snow on the GCE. In warm summers following dry winters, the snow line rises sooner and higher up-glacier, allowing more ice to melt over a larger area compared with years with cool (or also warm) summers following snowy winters. In addition, in long time series of streamflow, not only relative glacier cover can change but also the area of bare ice when the snow line rises transiently (Collins, 2006b;Dahlke, Lyon, Jansson, Karlin, & Rosqvist, 2014;Hock et al., 2005). Thus, the role of snow in the GCE is complicated. Further studies addressing particularly the role of snow within the GCE context could give more insights in the interannual streamflow variability that is not explained by glacier cover alone.
We explored potential relations with some other catchment characteristics that may explain the variability in the empirical relations of CV Q values with g ( Figures S1 and S2). Elevation characteristics, temperature, and catchment area showed some effect on CV Q . Mean elevation and temperature are however also related to g. Further model experiments and sensitivity analyses could give insights in the importance of other characteristics on streamflow variability in glacierized catchments, if models are able to represent hydrological (e.g., evapotranspiration and groundwater storage) and glaciological processes (glacier dynamics) in glacierized and non-glacierized catchments well.
Besides snow and catchment characteristics, antecedent conditions (for monthly flows) and changes and/or differences in climate characteristics might also be relevant to take into account. Annual mean temperature has been increasing over our analysed period, and precipitation amounts have been rather stable or slightly increasing.
However, the variability of T and P was fluctuating over time. If these fluctuations are large, glacier cover cannot compensate for them and changes in CV Q are then not only glacier cover related but also climate related (Collins, 2006a;Fleming et al., 2006). Another aspect of climate that could be important is the type of summer weather (Casassa et al., 2009;Moore, 1992). Moore (1992) suggests that in catchments with dry summers, in contrast to catchments with wet summers, the negative correlation between Q gl and Q nongl is more pronounced. In the Alps, the catchments have wetter summers compared with the summers in coastal North America, where most of the other studies on streamflow variability were conducted (e.g., Fleming & Clarke, 2005;Fountain & Tangborn, 1985;Moore, 1992). Combining different climatic influences within one regional sample may obscure a clear GCE curve. In our analysis, precipitation seasonality is similar for all catchments, but the relatively wet summers and varying precipitation amounts might complicate the search for a GCE curve. Related to that is the question of how the GCE appears in unstudied climates with distinctly differing amounts of precipitation and glacier melt and their interaction, possibly resulting in a different shape of a GCE curve.
Moreover, we might need to think of a relationship for the compensating effect of glaciers that also includes climate characteristics on the axes (like the Budyko curve) to be robust for changes in climate that are not only influencing g, and to distinguish different climates.
Such a relationship would then better allow for space-for-time substitution interpretations.

| CONCLUSION
This study analysed the potential to identify a universal relationship between glacier cover and interannual streamflow variability using complementary analyses of empirical data and a model experiment.
The nested catchment model experiment resulted in a characteristic convex GCE curve confirming the theory that streamflow variability is dependent on relative glacier cover. In this idealized case, the relative glacier cover determines if precipitation or temperature dominates streamflow variability as shown by the correlation of streamflow with precipitation and temperature. The optimum glacier cover corresponding to the lowest interannual streamflow variability derived from the model experiment was around 10% or 15% glacier cover, for annual and August streamflow, respectively. In contrast, the empirical analysis based on real observations from a multicatchment data set showed a considerable scatter, especially the hypothesized increased streamflow variability with increasing glacier cover was not evident.
Comparing the GCE curve from the multicatchment data with that from the model experiment shows that relative glacier cover does not completely overrule other factors that could influence streamflow variability, such as other catchment or glacier characteristics. This means that a relationship fitted to data from some (modelled) catchments cannot be transferred to catchments with other characteristics. Our results also showed that the relationship between glacier cover and interannual streamflow variability can change over time and that individual catchments do not clearly follow one curve over time. Differences between our results and previously published relationships shed doubt on the existence of an ideal and simple universal relationship between streamflow variability of different catchments with different glacier cover and therefore on the validity of the use of such a GCE curve for estimating streamflow variability in glacierized catchments.
Consequently, one needs to be careful about generalizing the role of glaciers as ideal buffers and assuming a characteristic GCE curve that can directly be used for space-for-time substitution applications.