The control of climate sensitivity on variability and change of summer runoff from two glacierised Himalayan catchments

The response of catchment runoff to climate forcing is determined by its climate sensitivity. We investigate the sensitivity of summer runoff to precipitation and temperature changes in winter-snow dominated Chandra (western Himalaya), and summer-rain dominated upper Dudhkoshi (eastern Himalaya) catchments in order to understand the nature of climatechange impact on the mean summer runoff and its variability. The runoff over the period 1980–2018 is simulated with a semi-distribute hydrologic model, which is calibrated using available discharge and glacier mass loss data. An analysis of 5 the interannual variability of the simulated summer runoff reveals that the runoff from the glacierised parts of the catchments is sensitive to temperature changes, but is insensitive to precipitation changes. The behaviour of the summer runoff from the non-glacierised parts is exactly opposite. Such precipitation-independent runoff from the glacierised parts stabilises the catchment runoff against precipitation variability to some degree. With shrinking glacier cover over the coming decades, the summer runoff from the two catchments is expected become more sensitive to the precipitation forcing and less sensitive to the 10 temperature forcing. Because of these competing effects, the impact of the glacier loss on the interannual variability of summer runoff may not be significant. However, the characteristic ‘peak water’ in the long-term mean summer runoff, which is caused by the excess meltwater released by the shrinking ice reserve, may lead to a detectable signal over the background interannual variability of runoff in these two catchments.

. A summary of the characteristics of Chandra and upper Dudhkoshi catchments. The meteorological variables are from bias-corrected reanalysis data (Hersbach et al., 2020), the hydrological data are from model simulations (the present study). The glacier mass-balance and area-loss estimates are from the existing literature (supplementary Table S3).

Catchment
Chandra Annual precipitation (mm yr −1 ) 1610 1531 Summer precipitation / winter precipitation 0.5 6.8 Liquid precipitation/ solid precipitation 0.5 9.7 Glacier area loss (% decade −1 ) 1.1-5.5 1.2-4.2 Glacier mass balance (m w.e. yr −1 ) −0.13 ± 0.11 to −0.56±0.38 −0.26 ± 0.13 to −0.52±0.22 Annual runoff (m yr −1 ) 1.25 0.99 the western disturbances (Bookhagen and Burbank, 2010), and the influence of Indian summer monsoon is relatively weak. In upper Dushkoshi catchment, more than 80% of the precipitation happens during the summer months ( Fig. 1e) with a dominant influence of Indian summer monsoon. Here, and in the rest of this paper, summer refers to the high-discharge period from May to September (Azam et al., 2019). Due to the contrasting seasonality of precipitation, the ratio of solid to liquid precipitation in Chandra and upper Dudhkoshi catchments are 0.5 and 9.7, respectively. Apart from the above differences in the precipitation 60 regimes, the two catchments are quite similar in terms of the mean annual temperature, catchment hypsometry, elevation range, specific summer runoff, glacier fraction, and the recent rates of glacier loss (Figs. 1d-1f, and Table 1).

Data and methods
We simulated river runoff of the above two catchments over the period 1980-2018 using VIC model (Liang et al., 1996). The model was forced with bias-corrected reanalysis data (Hersbach et al., 2020), and was augmented with a temperature-index

Climate sensitivity of runoff
The climate sensitivity of summer runoff Q is defined as the change in runoff due to a unit perturbation in a meteorological forcing parameter (e.g. Zheng et al., 2009). Here, we concentrate of the sensitivity of the specific summer runoff (Q) from the 75 studied catchments to changes in annual precipitation (P ) and mean summer temperature (T ). Note that summer runoff and annual precipitation are in m yr −1 , and mean summer temperature is in • C. Hereinafter, for any given variable X, we use the symbols X 0 , δX, ∆X, and σ X to denote the long-term mean, the anomaly for a given year, the change in the long-term mean, and the standard deviation, respectively.

Estimating the climate sensitivities 80
The sensitivities of summer runoff can be used (e.g. Zheng et al., 2009) to relate the anomalies of summer runoff (δQ), annual precipitation (δP ), and summer air-temperature (δT ) as follows.
Here, precipitation sensitivity is denoted by s P . = ∂Q ∂P = ∂δQ ∂P , and temperature sensitivity is denoted by s T . = ∂Q ∂T = ∂δQ ∂T . In Eq. (1), a possible bilinear interaction term proportional to δT δP (Lang., 1986) was not considered. We confirmed this correc-85 tion term, when included in the regression for the catchment studied, was not significant.
In order to estimate the sensitivity coefficients s T and s P , we regressed simulated time series of δQ for the catchments during 1997-2018 with the corresponding time series of δT and δP . The standard error of the fits obtained the corresponding uncertainties. The sensitivities estimated from the simulated δQ over 1997-2018 were validated using the simulated variability during 1980-1996, with the help of metrics Nash-Sutcliffe efficiency (NSE) and root mean squared error (RMSE). 90 We also considered the runoff from glacierised part of the catchments Q (g) .
= Q (g) 0 +δQ (g) , and that from the non-glacierised part of the catchments Q (r) .
= Q (r) 0 + δQ (r) . The corresponding sensitivities were defined in a similar way, These climate sensitivity coefficients and the corresponding uncertainties were estimated using the simulated annual anomalies 95 δQ (g) and δQ (r) , and the corresponding P and T anomalies over the period 1997-2018. With glacier fraction denoted by x, the following relations connect the quantities defined for the glacierised and non-glacierised part of the catchments to onesdefined for the whole catchment.
Apart from the sensitivities of summer runoff, we also computed the precipitation and temperature sensitivities of glacier mass balance using the corresponding simulated interannual variability. For the precipitation sensitivity of glacier mass bal-ance, we defined it to be the mass-balance change due to a 10% change in precipitation for the ease of comparison with the corresponding values available in the literature (e.g. Mölg et al., 2012). The climate sensitivities defined above determine the interannual variability of summer runoff given those of P and T as, where σ Q , σ P , and σ T are standard deviations of Q, P and T , respectively. An implicit assumption here is that δP and δT are uncorrelated over the simulation period, which we verified to be true at p < 0.05 level. Equation (7), together with Eqs. (5)   110 and (6), can be used to explain the first order differences in σ Q between different catchments in terms of the differences in the corresponding climate sensitivities, glacier fraction (x) , climate variability (σ T , and σ P ). The first order changes in σ Q under a changing climate can also be estimated using the above equation using projected values of σ T , and σ P . For the studied catchments, we used this relation to investigate the impact of changing glacier cover on the future interannual variability of summer runoff σ Q .

Multidecadal changes in mean summer runoff
The climate sensitivity coefficients defined above can also be used to predict changes in summer runoff (∆Q) over the next few decades for any given changes annual precipitation (∆P ) and mean summer temperature (∆T ). For glacierised catchments, apart from changing P and T , the changes in glacier fraction x will also influence the future runoff. Defining time-varying glacier fraction x . = x 0 + ∆x, the following linear-response equation can be constructed, In this linear-response approximation, the terms that are higher order in ∆ are ignored. The implicit assumptions in this formulation are: 1) the climate sensitivities of runoff from glacierised and non-glacierised parts do not change appreciably over a few decades, and 2) the contribution of a changing x to changes in summer runoff is well represented by the recent difference between the mean runoff of the glacierised and non-glacierised parts. We note that a similar linear-response approache been 125 used to analyse glacier-compensation effect (Chen and Ohmura., 1990), but without explicitly referring to climate sensitivity.
For Periche sub-catchment of upper Dudhkoshi catchment, possible changes in total runoff by 2050 were studied previously using hydrological model simulation (Soncini et al., 2016). We used the corresponding projected changes in P , T and x to obtain the summer runoff changes using Eq. (8) described above. These climate-sensitivity based predictions were compared with the available estimates from hydrological model simulations (Soncini et al., 2016) assuming the recent ratio of winter to 130 summer runoff to remain unchanged.
We investigated the ability of Eq. (8) to predict the timing and the magnitude of the 'peak water' (Hock et al., 2005;Huss and Hock, 2018) in the studied catchments. For these computations, the projected future changes in glacier area for Ganga and Indus basin for three different climate scenarios (RCP 2.6, 4.5, and 8.5) (Huss and Hock, 2018) were applied to Chandra, and upper Dudhkoshi catchments, respectively. The corresponding future temperature changes were obtained from projections available for the western and eastern Himalaya, respectively (Kraaijenbrink et al., 2017). The precipitation changes were not significant within the uncertainties (Kraaijenbrink et al., 2017), and were ignored here. The above estimates of the magnitude and the timing of the peak water were compared with the corresponding gridded values available in the literature (Huss and Hock, 2018). For this comparison of estimated peak-water, we considered the future changes in runoff of the parts of the catchments that were glacierised at 2000, following the convention of Huss and Hock (2018 Randolph Glacier Inventory (RGI 6.0) (Arendt et al., 2017) was used for the glacier boundaries in both the catchments. The boundaries corresponded to the glacier extent in 2002. We considered 8 available geodetic mass-balance observations for each 155 of the catchments that spanned a decade or more (see supplementary Table S3 for details).

Reanalysis data and bias correction
We used hourly 2m air temperature, precipitation, and wind-speed from fifth-generation European Center for Medium-Range Weather Forecasts Atmospheric Reanalysis of the global climate (ERA5) from 1980 to 2018 (Hersbach et al., 2020), at a spatial resolution of 0.25 • ×0.25 • , as model inputs. Following the existing hydrological studies of various high Himalayan catchments 160 (Immerzeel et al., 2013;Khadka et al., 2014;Soncini et al., 2016;Shrestha et al., 2017;Azam and Srivastava, 2020), the precipitation and temperature data were bias corrected as discussed below.
The available observed air temperature data at the Himansh station (Chandra catchment), and at Phadking (Dudhkoshi Catchment) were used to compute the mean monthly temperature biases for both the catchments (supplementary Fig. S1). The computed biases were corrected for, assuming they remain the same for the whole catchment and over the whole simulation 165 period. While comparing mean monthly temperature of the stations and that of the corresponding grids of ERA5, mean monthly lapse rates (supplementary Fig. S2) were used to correct for the station elevation.
To correct for possible biases in the ERA5 precipitation data, a constant scale factor α P was used for each of the catchments following the existing studies from the region (e.g. Bhattacharya et al., 2019;Krakauer et al., 2019;Zaz et al., 2019;Kanda et al., 2020;Azam and Srivastava, 2020). The scale factor was calibrated using the observed runoff and glacier mass-balance 170 data employing a Bayesian procedure as described in Sect. 3.2.4. In some of the existing studies from the region, an elevationdependent precipitation scaling have also been employed (e.g. Immerzeel et al., 2013). As an elevation-dependent correction may potentially introduce additional uncertainties (e.g. Johnson and Rupper, 2020), we preferred to use a constant α p for each catchments, keeping the number of calibration parameters to a minimum.
We scaled the incoming shortwave radiation obtained from VIC model by a catchment specific constant to match the cor-

Details of the hydrological model
We used VIC model (version 4.2.d, accessible from https://vic.readthedocs.io/en/master/; Liang et al. (1996)), forced by biascorrected ERA5 precipitation and temperature reanalysis data, to solve the energy and water balance and simulate runoff for 180 the above two catchments. Within this model, the partitioning of rain and snow was done based on a threshold temperature of Table S2).
The area-elevation distributions for both the catchments were obtained from 30 m resolution Shuttle Radar Topography Mission digital elevation model (Farr et al., 2007). We used 10 elevation bands for each grid box. Depending on the elevation range within a grid box, the size of the bands varied in the range 100-300 m. Soil data were from Harmonized World Soil While the two studied catchments have more than 20% glacier cover, VIC model does not have the capability to compute glacier melt or the evolution of the glacier extent (Liang et al., 1996). Here we augmented VIC model with a temperature-index glacier mass-balance module (Hock, 2003). As the observed shrinkage of glacierised fraction was small (0.01-0.05) over the simulation period for both the catchments (Table 1), we assumed a static glacier cover. Biases due to such static-glacier 195 approximation was shown to be small in another glacierised Himalayan catchment over the same period (Azam and Srivastava, 2020). Note that a dynamic description of glaciers would, however, be essential for predicting the long-term future changes in runoff as large changes in glacier extent are expected over the coming decades (e.g. Kraaijenbrink et al., 2017).
The hydrological fluxes over the non-glacierised parts were simulated with the standard VIC model. On the glacierised parts, a separate VIC model simulation was used to compute snow melt and snow-covered fraction in each elevation band. For the 200 snow-free glacier surface in each elevation band, a temperature-index model (Hock, 2003) with a catchment-specific degreeday factor (DDF) obtained the ice melt. For melt calculations, temperature of the individual elevation bands were obtained from bias-corrected ERA5 described before using mean monthly lapse rates derived from the same data set (supplementary Fig. S2).
The DDF values for each of the catchments were calibrated using a Bayesian procedure as described Sect.3.2.4. Glacier runoff was defined as the sum of snow melt, ice melt, and rainfall on the glaciers (e.g. Radić and Hock., 2014).

205
For model spin-up, we extended the meteorological input data back by repeating the data from 1980 to 1984. Subsequently, simulations were run over the period 1980-2018. Streamflow routing was implemented following Lohmann et al. (1998) to route the total surface runoff, baseflow, and meltwater fluxes generated from glacierised and non-glacierised parts of each grid box at an hourly time step.
Glacier mass balance was computed by subtracting the annual ice and snow melt over the glaciers from the corresponding 210 total snowfall. Complex mass-balance processes like avalanches (Laha et al., 2017) and debris effects (Kraaijenbrink et al., 2017) were ignored for the sake of simplicity.

Model calibration
With limited observed runoff data, the calibration of a large number of tunable model parameters in complex hydrological models may lead to issues like equifinality (Beven and Freer., 2001;Beven., 2006) and over-fitting, compromising the ability 215 of the model to capture the corresponding climate response accurately. Therefore we avoided calibrating a large set of model parameters (Jost et al., 2012) in the present study. We calibrated for only two model parameters: 1) catchment-wide precipitation scale factor α P , and 2) DDF of the glacier melt model. These two key parameters determine the catchment wide water balance, and the glacier mass balance. For the rest of the VIC model parameters, we used the central values of the recommended range (supplementary Table S2). Note that calibrated VIC model parameters used in the existing studies from the region (Zhao et al.,220 2015, 2019) were comparable to the ones used here to within 10-20%, suggesting that the VIC model parameters used in the present study were realistic.
To calibrate for the parameters α P and DDF, we used the following Bayesian approach (e.g. Tarantola, 2005). For given a set of available observations d and a set of model parameters θ, the posterior probability of the model parameters given the observations was, Here p(θ) was the prior distribution of the model parameters α P and DDF. assumed to be, Here the superscript obs (mod) denoted the observed (modelled) values for a given θ. Q i was the total summer runoff for the i-th year, and the summation was over all the years with observed runoff data for a given catchment. The uncertainty σ Q in summer runoff was taken to be 10% of the mean summer runoff based on existing estimates for other Himalayan rivers which 235 varied from 5% to 10% (Haritashya at al., 2010;Srivastava et al., 2014;Kumar et al., 2016). b j was the observed regional geodetic glacier mass balance, with the index j denoting individual records. The corresponding uncertainties were in the range 0.05-0.32 m w.e yr −1 (supplementary Table S3). The empirical factor of 1 2 in the exponent associated with the mass-balance observations in Eq. (10) ensured that the two exponential weights had similar magnitude for the most-probable model (e.g. Tarasova et al., 2016;Meyer et al., 2019).

240
For each of the catchments, a total 11 × 29 (319) model runs scanned the two-dimensional parameter space with step sizes of 0.2 for α P , and 0.5 mm • C −1 day −1 for DDF. These runs mapped out p(θ|d) over the parameter space. The most-probable pair of model parameters, i.e., the values of α p and DDF for which p(d|θ) was the maximum, was used for simulating the runoff and glacier mass balance. The weighted ensemble of all the 319 models obtained the corresponding 2σ uncertainties.
Due to the limited observed runoff data, the full runoff data set were utilised for the above calibration procedure without any 245 validation period. As an additional check, for upper Dudhkoshi catchment, the calibration was performed using data from four consecutive years, with remaining three year's data utilised for validation. This experiment was repeated 4 times using runoff data from different sets of four consecutive years for calibration.
Following earlier studies (e.g. He and Pang., 2015;Isenstein et al., 2015), the parameter sensitivity of simulated results was estimated with the help of additional 22 simulations where one of the 11 VIC model parameters (supplementary Table S2 The Bayesian calibration procedure obtained best-fit DDF values of 5.0 and 7.5 mm day −1 • C −1 for Chandra and upper Dudhkoshi catchments, respectively. These best-fit DDF values were in the same ballpark range as previously used in studies in and around Chandra (Azam et al., 2014(Azam et al., , 2019Pratap et al., 2019) and Dudhkoshi catchments (Pokhrel et al., 2014;Khadka et al., 2014;Nepal, 2016). The best-fit α P was 1.4 for both the catchments which was within the range of values 0.8-2.2 used 260 in the existing studies in the Himalaya to correct various reanalysis products (Bhattacharya et al., 2019;Krakauer et al., 2019;Zaz et al., 2019;Azam and Srivastava, 2020;Kanda et al., 2020).  The Bayesian calibration procedure which simultaneously fitted the glacier mass balance and the summer runoff data, yielded unique best-fit models for the catchments studied . In contrast, the existing studies from the region where DDF was calibrated using only discharge data, usually encountered equifinality issues (e.g. Azam and Srivastava, 2020) such that 265 the best-fit model was not uniquely determined.
The calibrated models reproduced the observed summer runoff reasonably well (Fig. 3) with RMSEs of 11-12% of the mean summer runoff, and NSEs of 0.88-0.80 for the two catchments. The RMSEs and NSEs obtained here are comparable or smaller than those reported in the existing studies from the region (Nepal, 2016;Mimeau et al., 2018;Bhattacharya et al., 2019;Azam et al., 2019;Azam and Srivastava, 2020). Four additional calibration experiments for upper Dudhkoshi catchment, each one 270 using a different set of 4 consecutive years of data for calibration, obtained comparable best-fit values of DDF (7.2±1.1 mm day −1 • C −1 ) and α P (1.43 ± 0.03), with NSEs and RMSEs similar to those mentioned above. In the sensitivity analysis, the absolute changes in summer runoff were less than ∼ 1.5% for the perturbed values of all 285 but two parameters (Figs. 2b-2d). The slightly higher summer-runoff sensitivities of 1.8-2.5 % for the longer time scales in the routing model became less than 1% When the annual runoff was considered. In the additional 80 simulations where two parameters were perturbed simultaneously, the resultant runoff sensitivities were approximately equal to the sum of those obtained in the corresponding pair of simulations with a single perturbed parameter (supplementary Fig. S4). supplementary Table S3). The RMSE between modelled and observed mass balance of Chandra (upper Dudhkoshi) catchment was 0.10 m w.e. yr −1 (0.11 m w.e. yr −1 ). Our simulations indicated that glacier runoff contributed 39±9% and 36±11% of the total summer runoff in upper Dudhkoshi and Chandra catchments, with th glacier ice loss amounting to 9% and 4% of the 295 respective total summer runoff. The estimated glacier contributions to runoff were largely consistent with the existing studies from the region when differences in fractional glacier cover were taken into account (supplementary Table S4).

Simulated glacier mass balance and its climate sensitivity
The sensitivity of the modelled glacier mass balance to temperature was −475 ± 93 and −274 ± 46 mm yr −1 • C −1 for Chandra and upper Dudhkoshi catchments, respectively. The corresponding precipitation sensitivities for these catchments were 200 ± 42 and 49 ± 20 mm yr −1 for a 10% change in precipitation. These sensitivities were all significant at p < 0.01 300 level, and were consistent with the known estimates for Chhota Shigri (Azam et al., 2014) and Zhadang (Mölg et al., 2012) glaciers, and available regional-scale values (Shea and Immerzeel, 2016;Wang et al., 2019) (supplementary Table S5). Only the reported value for Dokriani Glacier from the central Himalaya (Azam and Srivastava, 2020) was exceptionally high in comparison. Dudhkoshi (R 2 =0.93) catchments. The best-fit form for the calibration period 1997-2018 reproduced the variability of summer runoff over the validation period 1980-1996 reasonably well (Fig. 5) for both the catchments with RMSE < 0.04 m yr −1 and NSE > 0.93. These results indicated that a linear-response framework (Eq. (1)) provided a good description of the interannual variability of summer runoff in these catchments. We also confirmed that the interannual variabilities of annual precipitation and summer temperature were uncorrelated in these catchments at p < 0.05 level.

Climate sensitivities of catchment runoff
The best-fit temperature sensitivities of summer runoff s T were 117±8 and 116±34 mm yr −1 • C −1 for Chandra and upper Dudhkoshi catchments, respectively. The best-fit s P were 0.39±0.03 and 0.47±0.06 mm yr −1 mm −1 for Chandra and upper Dudhkoshi catchments, respectively. These sensitivities were all significant at p < 0.01 level. Note that s T and s P for both the catchments were the same up to the uncertainty, despite the differences in their climate setting. Also, in both the catchment 320 s P was significantly smaller than 1 mm yr −1 mm −1 , possibly indicating that an interannual change of the snow/ice storage in response to precipitation changes dampened of the precipitation sensitivity of summer runoff.
To the best of our knowledge, only two other studies from the Himalaya attempted to quantify the temperature sensitivity of catchment runoff. In Dokriani glacier, the reported s T of 620 mm yr −1 • C −1 (Azam and Srivastava, 2020) was significantly higher than the estimates presented above. This was possibly due to the exceptionally high temperature sensitivity of 325 glacier mass balance estimated for the glacier as discussed before, and a relatively higher glacierised fraction of 0.5. Previously reported s T from Dudhkoshi catchment was 5.7 ± 0.3 mm yr −1 • C −1 (Pokhrel et al., 2014) -an order of magnitude smaller compared to that obtained here and by Azam and Srivastava (2020). Our estimates of s P from both Chandra and upper Dudhkoshi catchments were comparable with those reported earlier (Pokhrel et al. (2014); Azam and Srivastava (2020); supplementary Table S6).  Interestingly, summer runoff from both winter-accumulation type glaciers in Chandra catchment and summer-accumulation type glaciers in upper Dudhkoshi catchment was approximately independent of the corresponding precipitation variability. This precipitation-independent glacier runoff is a novel finding that may not has been discussed in the existing glacio-hydrological studies in the Himalaya. Since the patterns of the climate sensitivities of the glacier runoff discussed above were similar for the two catchments despite their contrasting climate regime, it may be worth exploring if this is a general feature of the climate 340 sensitivity of glacier runoff in the Himalaya or elsewhere.
According to our simulations, snowfall constituted 97% and 86% of the total precipitation over the glaciersied parts of Chandra and upper Dudhkoshi catchments, respectively. Consequently, a positive precipitation anomaly did not translate into a higher summer runoff of the glaciers. For the years with high precipitation (equivalently, high snowfall), snow and ice melt output from the glaciers reduced due to a higher albedo, which reduced the available melt energy. This led to a negative 345 correlation (r = −0.5, p < 0.05) between simulated summer runoff and precipitation in these two catchments. These effects made the runoff of the glaciersied parts relatively insensitive to precipitation variability (upper Dudhkoshi catchment) or have a small negative precipitation sensitivity (Chandra catchment). The latter was likely related to a higher snow-to-rain ratio in the winter snow dominated Chandra catchment where any change in precipitation leads to a opposite change in glacier melt. For example, in years with a higher snowfall the glacier melt was suppressed due to a strong albedo feedback, leading to a slight 350 reduction of the glacier runoff.
The negligible s (g) P discussed above implied a stabilisation of the total runoff of the glacierised catchments against precipitation variability, as the runoff contribution from the glacierised fraction x was essentially independent of precipitation. The magnitude of this stabilising effect scaled with the glacier fraction (x) in the catchment. The above stabilising effect is consistent with the reported buffering of catchment runoff by glaciers during the extreme drought years across High Mountain Asia 355 (Pritchard, 2019).

The relatively high s (g)
T in both the catchments as presented above can be understood as follows. A higher temperature implied a higher available energy leading to higher meltwater flux from the glaciers. The glaciers effectively acted as infinite reservoirs over an annual scale so that the meltwater volume was limited only by the available energy. These arguments were consistent with a high correlation (r > 0.9, p < 0.05) between the summer temperature and summer runoff of the glacierised 360 parts for both the catchments.

Climate sensitivity of runoff from the non-glacierised parts
In the non-glacierised parts of Chandra and upper Dudhkoshi catchments, s T within the corresponding uncertainties. However, compared to 365 the glacierised parts, the climate sensitivities of the runoff from the non-glacierised parts showed an exactly opposite trend with summer runoff relatively insensitive to temperature anomalies and sensitive to precipitation anomalies. Because of the presence seasonal snow cover over the non-glacierised parts, a temperature dependence of the summer runoff may be expected. However, the amount of snowmelt was limited by the supply of seasonal snow, and not by the available energy. This led to a weak response of the total summer runoff from the non-glacierised parts to temperature forcing. This argument was supported by 370 the fact that summer runoff from the non-glacierised parts were uncorrelated with summer temperature and strongly correlated with summer precipitation (r > 0.9, p < 0.05).

Implications for future changes in the mean and variability of summer runoff
The estimated climate sensitivities from glacierised and non-glacierised parts of both the catchments suggested s δQ ≈ xs These simplified equation suggested that the key parameters that determine the climate reponse of the catchment to gievn forcing were s (g) T and s (r) P . Below we first analyse the predictive power of Eq. (15), and then discuss the implications of the above simplified linear-response description in the context glacier-compensation effect (e.g. Chen and Ohmura., 1990) and peak-water effect (e.g. Huss and Hock, 2018) for the catchments.  (Soncini et al., 2016). The corresponding reported changes in summer runoff were A) 0.10, and B) 0.12 m yr −1 (Soncini et al., 2016) . For the above to scenarios, 390 Eq. (15) yielded similar summer runoff changes of A) 0.08±0.01, and B) 0.10±0.02 m yr −1 . This indicated the reliability of the climate sensitivities of summer runoff presented here, and also confirmed that accurate climate sensitivities can be used for an accurate prediction of catchment response to given forcing over multidecadal scales with minimal computational cost (e.g., Vano and Lettenmaier, 2014).  Then,Eq. (14) implies that σ Q is a non-monotonic function of x (Fig. 6a). It has a high value in the limit x → 0 due to a strong precipitation sensitivity of the off-glacier runoff, with σ Q ≈ (1 − x)s (r) P σ P . In the opposite limit of x → 1, σ Q is again high due to a high temperature sensitivity of glacier runoff, with σ Q ≈ s (g) T σ T . These to competing effects imply a minimum value of σ Q at an intermediate value of x (Fig. 6a). This is the well known glacier compensation effect (e.g., Chen and Ohmura., 1990). 400 Thus, Eq. (14) provides a theoretical explanation of the glacier-compensation effect. Note that Eq. (14) suggests a hyperbolic behaviour of σ Q (x), while some of the existing studies used an empirical parabolic curve to describe the depndence of the coefficient of variation of runoff from glacierised catchments on x (e.g., Chen and Ohmura., 1990). The results discussed above are also consistent with a reported strong correlation between runoff and precipitation (temperature) in the limit of small (extensive) glacier cover (Van Tiel et al., 2020).

405
It was proposed by Chen and Ohmura. (1990) that the parabolic dependence of the coefficient of variation of runoff on glacierised fraction can be utilised to estimate the change in σ Q for any given change in glacier cover. However, it has been reported that the corresponding best-fit parabola changes with time, making such predictions unreliable (Van Tiel et al., 2020). This effect can be explained based on Eq. (14) as follows. It is apparent from Eq. (14) that together with the glacier cover x, if σ P and σ T were to change under a changing climate, then the Eq. (14) describes a different hyperbolic curve altogether. This 410 can be verified from fig. 6a: the simulated σ Q for both the catchments for the period 1980-1996 were higher than those during the period 1997-2018, even as x remained constant in our simulations. These changes were due to differences in σ P and σ T in the corresponding period. Note that the simulated σ Q values during 1980-1996 could be predicted well using Eq. (14) together with σ P and σ T values over the same period (Fig. 6b). In fact, it is apparent from Fig. 6a that for the two studied catchments, the future changes in σ Q are likely to be dictated by those of climate variability (i.e., σ P and σ T ). The influence of a changing 415 glacier cover on the runoff variability may be relatively minor. Overall, the above exercise indicated that Eq. (14) may be used to quantify the changes in the variability of future runoff, and the departures of observed/modelled catchment response from the glacier-compensation-effect curve (Van Tiel et al., 2020).

Peak water due to shrinking glaciers
The future changes of glacier runoff from Chandra and upper Dudhkoshi catchments as obtained from Eq. (15) captured the 420 'peak water' successfully ( Fig. 7). For all the three climate scenarios (supplementary Fig. S6), an initial steady increase in the summer runoff from the glacier was followed by a steady decline. In Chandra catchment, the estimated timing of peak water were 2033, 2045, and 2055, respectively, for the RCP 2.6, 4.5, and 8.5 scenarios. The corresponding peak glacier runoff were 1.73, 1.85, and 1.97 m yr −1 , compared to the recent value of 1.54 m yr −1 (supplementary Table S7). In upper Dudhkoshi catchment, the estimated timing of peak water were 2022, 2031, and 2034, respectively, for the above scenarios.

425
The corresponding peak summer runoff from the glaciers were 1.74, 1.80, and 1.86 m yr −1 , compared to the recent value of 1.59 m yr −1 (supplementary Table S7). These predicted 'peak water' magnitude may not lead to a observable signal in the catchment runoff given that the present σ Q values of 0.08 (upper Dudhkoshi catchment) and 0.13 (Chandra catchment) m yr −1 (Fig. 6a) in these catchments. Except for the more ambitious RCP 2.6 scenario, in other scenarios the 'peak water' signal is likely to be observed in the catchment runoff over the background interannual variability of runoff (σ Q ) in these two 430 catchments.
It is interesting that the simple climate-sensitivity based estimates obtained timing and amplitude of the 'peak water' that were comparable with the previously reported values (supplementary Table S7) computed using state-of-the-art glaciohydrological simulations (Huss and Hock, 2018). In upper Dudhkoshi catchment, however Eq. (15) predicted somewhat quicker and smaller 'peak water' compared to the previously reported value (supplementary Table S7). This inconsistency may be re-435 lated to the use of different glacier inventories and meteorological forcing. For example, Huss and Hock (2018) obtained a specific glacier runoff in Ganga basin which was about twice as large as obtained here for upper Dudhkoshi catchment. Also, itation changes were ignored in the present study due to the large uncertainty. We note that the present formulation (Eq. (15)) ignored the feedback of the evolving glacier geometry on mass balance, and thus, possible future long-term changes in Q (g) 0 .
440 However, despite the limitations, the above exercise overall indicated that a climate-sensitivity based approach was able to capture the basic characteristic of the long-term response of glaciersied catchment, namely, the peak water effect, both qualitatively and quantitatively.

Model limitations
An obvious limitation of the present study is a lack of long-term hydro-meteorological field data. However, this is a common 445 problem that affects glacio-hydrological studies in the Himalaya. We have used runoff time series generated with VIC model simulations to circumvent the issue. A bias correction of the reanalysis product using the available meteorological field data, the use of total summer runoff and the decadal-scale glacier mass balance data for calibration, a Bayesian calibration procedure, calibrating for only two parameters to avoid over-fitting ensured robustness of the model results presented. The model sensitivity to the specific choice of parameters were generally low. Available estimates from the existing literature with the 450 corresponding results presented here indicated reasonable model performance as well.
We reemphasise that the climate sensitivity discussed are only applicable over the range of forcing used in the corresponding fits. Over shorter time scale of multiple decades, as long as the changes in forcing is within the range of the fits, the climate sensitivities presented here may be expected to characterise the catchment response accurately. Also, over longer time scales the dynamic glacier geometry which was ignored in our simulation may lead to systemtic changes in Q (g) 0 that cannot be 455 captured in a linear response approach.

Summary and conclusions
In this paper, we simulate the runoff of Chandra (western Himalaya) and upper Dudhkoshi (eastern Himalaya) catchments over 1980-2018, using VIC model augmented with a temperature-index glacier-melt module. With a minimal calibration strategy of calibrating only two parameters to fit the available data of summer runoff and decadal-scale geodetic glacier mass balance, 460 the model produces good match with both observations and the results available in the literature. The interannual variability of summer runoff simulated with the VIC model is then utilised to obtain the climate sensitivity of summer runoff to temperature and precipitation forcing in the catchments, which leads to the following conclusions.
-For both the catchments, summer runoff from the glacierised parts have a high temperature sensitivity, and that from the non-glacierised parts is essentially temperature independent. Temperature sensitivity of summer runoff in Chandra and -Precipitation sensitivity of summer runoff from the non-glacierised parts are high, but that of the runoff from the glacierised parts is negligible. Precipitation sensitivity of summer runoff in Chandra and Upper Dudhkoshi catchments are 0.39 ± 0.03 and 0.47 ± 0.06 mm yr −1 mm −1 , respectively.
-In the limit of low glacier cover, an increasing glacier cover stabilises catchment runoff against precipitation variability 470 due to a precipitation-independent runoff of the glaciers. In contrast, at the limit of high glacier cover, the variability of catchment runoff increases with glacier cover due to a strong temperature depndence of glacier runoff. This leads to the so-called 'glacier-compensation effect'. The present climate sensitivities and the future variability of precipitation and summer temperature, will determine together the future variability of summer runoff.
-A knowledge of the climate sensitivities of runoff of the glacierised and the non-glacierised parts allows an estimation of 475 future changes of catchment runoff and the amplitude and the timing of 'peak water', under any given climate scenario.
Data availability. All the observed hydro-meteorological data, basic simulation results, and codes developed will be made public.
Author contributions. AB conceived the study, and developed the theoretical framework. SL ran the simulations. SL and AB did the data analysis. AB wrote the manuscript with inputs from SL. AS, PS, and MT did the field measurements at Chandra catchment, and analysed the field data. All the authors participated in the discussions, and edited the manuscript.