We would like to thank the reviewer #1 Jonathan D Mackay for taking the time to provide comments on the manuscript. Please find our replies below.

Reply to reviewer #1  

A well written and clearly presented manuscript that introduces an updated version of the DECIPHeR hydrological model which includes a simple energy balance model to simulate snow and ice ablation and accumulation. Overall, only a few comments from me which should be easy for the authors to address (see attached PDF with comments). Of these comments, two are slightly more "major" points which I hope the authors can address in their response and manuscript revisions. The first is on the question of novelty. It's not entirely clear from the introduction what the novelty of the work is exactly. In the Conclusions section it states that:

"The motivation for this work was to develop hydrological model that can be used to simulate discharge in very large glaciated and snow-fed catchments, at a high spatial resolutions, whilst maintaining the ability to explore model uncertainty."

So I'm assuming that the novelty is the model itself, but my understanding is that there are already models out there that can be used to do this (I've mentioned some in the attached). Could the authors please spell out what the novelty of the work is in the introduction. If the novelty is the model then I think the manuscript would really benefit from a more explicit explanation of the limitations of current models available and what exactly this model offers to address these. A good starting point might be the review of Van Tiel et al. (https://wires.onlinelibrary.wiley.com/doi/full/10.1002/wat2.1483).

To clarify the novelty of the work is in using the DECIPHeR model.

We added the following text (L69) of the introduction to highlight the reasons we used DECIPHeR, instead of one of the other glacio-hydrological models in the literature.

“Many glacio-hydrological models already exist in the literature (van Tiel et al., 2020; Horton et al., 2022), however, we integrate a snow and glacier melt model into DECIPHeR for the following three reasons.  Firstly, DECIPHeR uses hydrological response units (HRUs) to model water flow in hydrologically similar parts of the catchment which allows the model to be run as a fully distributed (HRU for every single grid point), semi-distributed (multiple HRUs) or as a lumped model (1HRU). Depending on user requirements and the corresponding degree of complexity, topographic, land use, geology, soils, anthropogenic and climate attributes as well as points of interest (any gauged or ungauged point on the river network), can be supplied to define the spatially connected topology and thus differences in model inputs, structure and parameterization (Coxon et al. 2018). Other HRU based glacio-hydrological models exist, for example, SWAT (Omani et al., 2017), PREVAH (Koboltschnig et al., 2008) and HBV (Finger et al., 2015) but they don’t offer this level of flexibility within a single modelling framework.

Secondly, DECIPHeR is computationally efficient, which makes it suitable for modelling very large catchments. Many of the glacio-hydrological models in the literature are distributed (grid point based) for example, TOPKAPI (Pellicciotti et al.,2012), DHSVM (Frans et al., 2018), VIC (Schaner et al., 2012), GERM (Farinotti et al., 2012). The computational expense of modelling processes with adjacent grid points makes distributed models more suited to studying small catchments. Furthermore, computationally efficiency makes it possible to quantify uncertainties and run large ensembles which is important for understanding the uncertainties in future predictions.

Thirdly, the DECIPHeR code is open source which allows opportunities for further community development. In contrast, the glacier enhanced version of SWAT (Omani et al., 2017; Luo et al., 2013) is not open source.”

The other point regards the application of GLUE and the use of the top 0.5% of model simulations to represent a population of behavioural models.  The 0.5% seems arbitrary and, therefore, it's not clear to me what the merit of including these in the analysis is.  What do the uncertainty bounds of an arbitrary population of models mean? Could the authors please justify the use of using the top 0.5% of simulations rather than, say, defining a more objective set of "good behaviour" criteria e.g. based on the different metrics of model performance used in the study.

We added this to the discussion section “

In this study we set the behavioural models to the best 0.5% simulations in the ensemble as it was  important in our analysis to rank models according to their ability to capture seasonal discharge, particularly from spring snow melt and summer glacier melt. Often behavioural models are selected using threshold values for guidelines metrics. These metrics are calculated over the complete discharge timeseries, rather than for individual seasons. For example, metrics from Moriasi et al. (2007) are commonly used in the literature to categorise ‘acceptable’, ‘good’, or ‘very good’ simulations based on threshold values for NSE, PBIAS and RSR. Metrics calculated over the complete discharge time series are not a strong test of the model’s ability to predict seasonal discharge. To our knowledge, there are no standardised guideline thresholds in the literature for seasonal metrics, therefore we selected the best 0.5% of the ensemble. If we decided to define the behavioural models using a threshold for the seasonal RSR, then this would also be based on an arbitrary choice of value. A high threshold for seasonal RSR would be required to categorise the behavioural models because the summer values are high (See RSR_JJA in Table 4).

We explored the impact of selecting alternative threshold values (1%, 2.5%, 5% and 10%) on the calibrated NSE values (Fig. S24). To obtain NSE values > 0.7 at all the gauging stations requires a threshold smaller than 1%. This is notable at the Alfatun station where the NSE value at the 0.5% threshold is 0.74 but reduces to 0.63 at the 1% threshold (95^th percentile limit values). Fig S24 also shows how the uncertainty in the NSE values increases for higher threshold values. At the 10% threshold the uncertainties in the NSE values are much greater than at the 0.5% threshold.”

Please find attached an updated Supplementary material and replies to the minor comments in the zip file.     

Additional references

Aizen, V. B., Aizen, E. M., and Melack, J. M. (1995). Climate, Snow Cover, Glaciers, and Runoff in the Tien Shan, central Asia. JAWRA J. Am. Water Resources. Assoc. 31. doi:10.1111/j.1752-1688.1995.tb03426.x

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, HYDROLOGICAL PROCESSES, 26, 1909–1924, https://doi.org/10.1002/hyp.8276, 2012

Finger, D., Vis, M., Huss, M., and Seibert, J.: The value of multiple data set calibration versus model complexity for improving the performance of hydrological models in mountain catchments, Water Resources Research, 51, 1939–1958, https://doi.org/https://doi.org/10.1002/2014WR015712, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2014WR015712, 2015

Frans, C., Istanbulluoglu, E., Lettenmaier, D. P., Fountain, A. G., and Riedel, J.: Glacier Recession and the Response of Summer Streamflow in the Pacific Northwest United States, 1960-2099, WATER RESOURCES RESEARCH, 54, 6202–6225,720 https://doi.org/10.1029/2017WR021764, 2018

Horton, P., Schaefli, B., and Kauzlaric, M.: Why do we have so many different hydrological models? A review based on the case of Switzerland, WIREs Water, 9, e1574, https://doi.org/https://doi.org/10.1002/wat2.1574, https://wires.onlinelibrary.wiley.com/doi/abs/10.7601002/wat2.1574, 2022

Koboltschnig, G. R., Schöner, W., Zappa, M., Kroisleitner, C., and Holzmann, H.: Runoff modelling of the glacierized Alpine Upper Salzach basin (Austria): multi-criteria result validation, Hydrological Processes, 22, 3950–3964, https://doi.org/https://doi.org/10.1002/hyp.7112, https://onlinelibrary.wiley.com/doi/abs/10.1002/hyp.7112, 2008

Luo, Y., Arnold, J., Liu, S. Y., Wang, X. Y., and Chen, X.: Inclusion of glacier processes for distributed hydrological modeling at basin scale with application to a watershed in Tianshan Mountains, northwest China, Journal of Hydrology, 477, 72–85, https://doi.org/10.1016/j.jhydrol.2012.11.005, <GotoISI>://WOS:000313935200006, 2013

Moriasi, D. N., Arnold, J. G., Van Liew, M. W., Bingner, R. L., Harmel, R. D., and Veith, T. L.: Model evaluation guidelines for systematic825 quantification of accuracy in watershed simulations, Transactions of the Asabe, 50, 885–900, <GotoISI>://WOS:000248036800021, 2007

Omani, N., Srinivasan, R., Karthikeyan, R., and Smith, P. K.: Hydrological Modeling of Highly Glacierized Basins (Andes, Alps, and Central Asia), Water, 9, https://doi.org/10.3390/w9020111, <GotoISI>://WOS:000395435800041, 2017

Pellicciotti, F., Buergi, C., Immerzeel, W. W., Konz, M., and Shrestha, A. B.: Challenges and Uncertainties in Hydrological Modeling of Remote Hindu Kush-Karakoram-Himalayan (HKH) Basins: Suggestions for Calibration Strategies, MOUNTAIN RESEARCH AND DEVELOPMENT, 32, 39–50, https://doi.org/10.1659/MRD-JOURNAL-D-11-00092.1, 2012

Saks, T., Pohl, E., Machguth, H., Dehecq, A., Barandun, M., Kenzhebaev, R., Kalashnikova, O., and Hoelzle, M.: Glacier Runoff Variation Since 1981 in the Upper Naryn River Catchments, Central Tien Shan, Frontiers in Environmental Science, 9, https://doi.org/10.3389/fenvs.2021.780466, https://www.frontiersin.org/articles/10.3389/fenvs.2021.780466, 2022

Schaner, N., Voisin, N., Nijssen, B., and Lettenmaier, D. P.: The contribution of glacier melt to streamflow, ENVIRONMENTAL RESEARCH LETTERS, 7, https://doi.org/10.1088/1748-9326/7/3/034029, 2012

van Tiel, M., Stahl, K., Freudiger, D., and Seibert, J.: Glacio-hydrological model calibration and evaluation, WILEY INTERDISCIPLINARY REVIEWS-WATER, https://doi.org/10.1002/wat2.1483, 2022

We would like to thank the reviewer #2 for providing comments on the manuscript. Please find our replies below.

The paper is dedicated to an acute problem of development of glacio-hydrological models for the prediction of future changes in river runoff due to deglaciation. The presented study aims to develop a computationally efficient hydrological model that can be applied to large glaciated and snow-fed catchments. The paper is overall well-written and provides interesting results. However, there are few major and several minor recommendations to the authors, stated bellow:

The description of the DECIPHeR model needs to be extended: what hydrological processes are taken into account, how the water is routed, number of conceptual storages etc.

We have extended the description of DECIPHeR in Section 2.2 as follows:

“DECIPHeR simulates water storage, hydrologic partitioning, and surface/subsurface flow for steeper shallow soils and/or groundwater-dominated watersheds. The model structure (as implemented in Coxon et al, 2019) consists of three stores defining the soil profile (root zone, unsaturated and saturated storage), which are implemented as lumped stores for each HRU. Moisture is added to the soil root zone by rainfall input and removed only by evapotranspiration. Any excess precipitation is added to the unsaturated zone where it is either routed directly as overland flow or added to the saturated zone. Changes to storage deficits in the saturated zone are dependent on this recharge from the unsaturated zone, fluxes from upslope HRUs and downslope flow out of each HRU. Subsurface flows for each HRU are distributed according to a flux distribution matrix based on accumulated area and slope. Channel flow routing is modelled using a set of time delay histograms. For more detailed discussion of the original DECIPHeR model structure please see Coxon et al, (2019).

While DECIPHeR has been applied catchments in the UK (Coxon et al., 2019; Lane et al.,2021), it has not been used in a glaciated or snow-fed regions because these processes have not yet been included in the model. ….”

A clearer parameters calibration scheme should be added to the methods section. What is the initial and resulting range of the parameters? It is mentioned that degree day factor varies daily in the introduction – it gives the first impression that the values are calibrated for each day separately.

We have added an extra table to the Supplementary Material (Table S3) which lists the initial and calibrated parameter ranges for each sub-catchment. The table also lists the parameter values for the overall best simulations. We reference Table S3 in Section 3.2.2. The updated Supplementary Material is in the zip file attached to the replies to reviewer #1. 

We removed reference to the degree day factor varying daily in the introduction (L130), so as not to mislead the reader. The calculation of the degree day as a function of day number is described in Equation 3.

The 3.1 section provides information on the evaluation and validation period. It seems that for the evaluation the same period as for the calibration was used? It is not quite common. Authors should comment on that.

We replaced instances of “evaluation” with “calibration” to make clear that the period to which we refer to, is the calibration period (1951-1970).  

It is mentioned on P16 L 337 that the Nash–Sutcliffe efficiency (NSE) is used to evaluate high flows and the timing of peak discharge. Just below that a formula for mean monthly discharges evaluation using NSE metric is given

To clarify the point, we are using the NSE calculated from all monthly values as a metric to evaluate high flows and the timing of peak discharge. High NSE values will be found if the discharge predicts the monthly peak flows and timing of the monthly discharge well. We are not calculating NSE using the high flow values only (i.e. a subset of the monthly data).

Analysis of model performance using the MSC method compared to ISC method for other sub-catchments should be included as well in 3.2.2

We added table S2 to the supplementary material which lists the ranges of values for precipitation lapse rate and sublimation correction factors for the 10 best performing experiments for the individual sub-catchments.  

The text in section 3.2.2 is edited.

“Simulations that perform well in the sub-catchments (Figs. S10-S15) favour higher values for the precipitation lapse rates, in contrast to the global catchment parameters which range from 1%100m−1 − 10%100m⁻¹ (Fig. S9). This is visible in Table S2 which summarises the range of precipitation lapse rates for the 10 best performing simulations for each sub-catchment. The upper values for the precipitation lapse varies between 16 and 24%100m⁻¹ depending on the sub-catchment, which is higher than the global catchment upper bound of 10%100m⁻¹. Simulations also perform better in the sub-catchments when higher values for the sublimation factor E_sub are used, in contrast to the global values (0.005 - 0.2). The 10 best E_sub parameter ranges are also listed in Table S2. The upper bound values for E_sub vary between 0.6 - 1.0, depending on the sub-catchment, which is higher than 0.2 predicted by the global catchment values. E_sub controls the reduction in PET over snow and ice surfaces. “

Compare the range in glaciated area prediction with the observed glaciated area

We added this at L445

 “The model produces a large range of estimates for the glaciated area (680km² − 1, 196km²) (5^th -95^th percentile limits) at the end of the simulation period. This range is larger than the observed uncertainty range of 903 - 948 km².  The uncertainty range in the model is 516 km² (in 2007) which is more than 10 times greater than the uncertainty in the observed glaciated area (46 km²)”

The positive trend in snow melt and negative trend in rainfall component seems to be consistent over the territory that could be better emphasized in the text.

We added this “Small positive trends in the snow melt fraction are consistent across the catchment and are likely driven by warming temperatures. There are small negative trends in the rainfall fraction of less than 1% per decade which is associated with a small decrease in the APHRODITE precipitation. “

Discussion should be extended covering following aspects: 1) the 95th percentile simulations in all cases show an asymmetrically larger contribution of the rainfall compared to 5th and 50th percentile,

We add the following text at L501

“The discharge components are calculated using the 0.5% best calibration parameters for the Uch-Kurgan station located the outlet of the catchment.  “

Added at L507

“ Figure 10 shows that the rainfall component is larger at the 95^th percentile simulations than at the 5^th and 50^th percentile simulations. This is because the lapse rate at the 95th percentile simulations is higher (22%100m^-1) than at the 5th (1%100m^-1) and the 50^th (6%100m^-1) percentile simulations. “

2)analysis of the importance of including new calibration parameters in the DECIPHeR model. As the model performance seems to be not very sensitive to most of the calibration parameters values (FigS15),

We have added this into the discussion on the section on improving the evaluation to include additional observations (L605) as follows:

“This highlights the importance of including ancillary observations, such as glacier mass balance, snow depth or snow extent, in the evaluation to help constrain the predictions and parameter values. Currently, the model performance is not sensitive to many of the calibration parameter values (Fig S15). It is possible that some parameter combinations compensate each other. For example, a high snowfall correction factor may be compensated for by a lower precipitation lapse rate. More analysis needs to be conducted on the sensitivity of the new snow and ice parameters added to DECIPHeR as part of this study, both in time and space, and the types of data that may help to constrain these parameters. Remote sensing snow products have been used to evaluate models and studies indicated that the integration of data such as MODIS snow cover into hydrological models can improve the simulated snow cover while maintaining model performance with respect to runoff (Parajka and Bloschl, 2008)."

3) comparison of derived contributions of snow melt, glacier melt, rainfall with previous studies

We added this to the discussion at L529

“We used the model to calculate the relative contributions of snow, rain and glacier melt to the annual runoff. We found spatial variability in the relative contributions of each of the components. For the entire catchment (gauging station at Uch-Kurgan) the 50^th percentile contributions are snow (89%), rain (9%) and glacier melting (2%). These estimates are broadly consistent with Armstrong et al. (2019) who used MODIS imagery and degree day melt modelling to partition the runoff components in the Syr Darya river. Armstrong et al. (2019) found the runoff comprised of snow (74%), rain (23%) and glacier melting (2%). Our estimates are slightly higher for the snow melt contribution; however, our study focuses on the upper reaches of the Syr Darya river where the snow melt is more likely to dominate the runoff.

Snow melting is the dominate component of the runoff at the six gauging stations. Throughout the Tien Shan long-term hydrological records of the former USSR show that snow melt is the dominant source of runoff (Aizen et al., 1995). Further upstream in the Naryn sub-catchment the glacial melt contribution to the annual runoff is higher (4% - 15%) than at Uch-Kurgan. Our upper estimate (15%) is slightly lower than a study by Saks et al. (2022) who calculated that 23% of the runoff originates from glacier melting in upper Naryn river. A possible explanation for why our estimate is lower, is that our simulation period starts 30 years earlier (1951) than the study by Saks et al. (2022) which started in 1981.”

Other minor suggestions and technical corrections:

P 1 L12 The model reproduces the spatial extent in seasonal snow cover well, capturing 86% of the snow extent on average (2001-2007) for the median ensemble member of the best 0.5% evaluation simulations, when evaluated against MODIS snow extent. Better divide the sentence in 2-3 sentences to make the message clearer.

Rewritten as

“The model reproduces the spatial extent in seasonal snow cover well, when evaluated against MODIS snow extent. 86% of the snow extent is captured (mean 2001-2007) for the median ensemble member of the best 0.5% evaluation simulations.” 

L18 At all stations snow melting is the largest component, followed by the rainfall and the glacier melt component. Please provide estimation of shares

Added values “Snow melting is the largest component of the annual discharge (89%), followed by the rainfall (9%) and the glacier melt component (2%), where the values refer to the 50^th percentile estimates at the catchment outlet gauging station Uch-Kurgan. “

P2 L30 Sry Darya – Syr Darya

Typo corrected

 L32 semi-arid lows lands – semi-arid low lands

Typo corrected

 P3 L78-79 Section 3 describes the evaluation and validation of discharge. The evaluation and validation of modelled runoff?

We are calibrating and validating against discharge observations. 

Section 4 describes the validation of snow extent against MODIS observations of modelled snow extent?

Changed section heading from “Validation of snow extent against MODIS observations”

to “Validation of modelled snow extent against MODIS observations”

P4 L96-98 A high resolution irrigation map of the catchment derived from normalized difference vegetation index (NDVI) (Meier et al., 2018) shows that the irrigated area is low, in contrast to the Ferghana valley downstream. (Fig. S1). It would be better to add the numerical estimation to the comparison

We calculated that 3% of the catchment is irrigated and added this to the text.

“A high-resolution irrigation map of the catchment derived from normalized difference vegetation index (NDVI) (Meier et al., 2018) shows that the irrigated area is low (3% area is irrigated), in contrast to the Ferghana valley downstream. (Fig. S1).”

 P8 L163 “d” symbol doesn’t seem the best choice for the day of the year. As it is hard to distinguish from the first site in a formula abundant with letters d

Changed “d” to “j”

P26 L428 The dotty plots show… –> The dotty plots (see Fig.S9) show..

Edited as suggested.

Figure 3. Please add the transcription of the used indexes either in the caption or in the text.

We added a table containing the definitions of the symbols in Figure 3. The table is included in an appendix (inserted at the end of the paper)

Table 7. The addition of p-values would probably contribute to the informativeness of the table

This would perhaps make the table a bit messy. The p-Values < 0.05 are indicated by the bold font.

S1 The colors need to be explained

We edited the figure to include a legend.

S2 The color ramp is evidently different for the left and right half of the picture

The line visible at approximately 77.75^oE longitude is because the HRUs share a common climate grid, rather than an error in the colormap. The pattern is visible because the HRUs share the same climate input grid (APHRODITE and ERA5 data).  

 S3-S4 Glacier thickness seems to differ a lot between the pictures, though the corrections only for two glaciers are mentioned.

We have replotted these figures with a consistent colour ramp. The thickness is the same except for the corrected regions at the two glacier snouts.

Additional references

Aizen, V. B., Aizen, E. M., and Melack, J. M. (1995). Climate, Snow Cover, Glaciers, and Runoff in the Tien Shan, central Asia. JAWRA J. Am. Water Resources. Assoc. 31. doi:10.1111/j.1752-1688.1995.tb03426.x

Saks, T., Pohl, E., Machguth, H., Dehecq, A., Barandun, M., Kenzhebaev, R., Kalashnikova, O., and Hoelzle, M.: Glacier Runoff Variation Since 1981 in the Upper Naryn River Catchments, Central Tien Shan, Frontiers in Environmental Science, 9, https://doi.org/10.3389/fenvs.2021.780466, https://www.frontiersin.org/articles/10.3389/fenvs.2021.780466, 2022