The paper is interesting to read, which presented the development of Water and Energy transfer Processes model in the Qinghai–Tibet Plateau (WEP-QTP) that modified based on the original Water and Energy transfer Processes model in Cold Regions (WEP-COR). In the presented model, the vadose zone processes considered three strategies under different conditions: (1) a dualistic soil-gravel structure using the Richards equation under non-heavy rain in the nonfreeze–thaw period; (2) a multi-layer Green-Ampt model in a heavy rain scenario in the nonfreeze–thaw period; and (3) a hydrothermal coupling model based on the continuum of the snow-soil-gravel layer during the freeze––thaw period. The modified model was then verified with measured river discharge in Niyang River Basin by comparing the simulated groundwater.

The study adopted a new conceptualization of the water and energy transfer in Qinghai–Tibet Plateau, which is considered as novel. However, significant improvement is needed before the consideration of publication to Hydrology and Earth System Sciences.

Dear Reviewer:

We appreciate the detailed and valuable comments, which have considerably improved the quality of our manuscript. Our responses to the comments are provided below.

Major comments:

1. Author stated in Line 52-62 that the existence of gravel in soil can significantly affect the soil water content and water transport. However, coupling of soil water and heat transport may be still not fully achieved in the modified version of WEP-QTP. When the dualistic soil-gravel structure was used in the nonfreeze–thaw period, the soil water transport may be decoupled with the thermal transport (see Line 296: “for soil and gravel layers, the average temperature was represented by the temperature of the middle layer”). It seems that the full coupling of water and heat transport can only be achieved for freeze––thaw period? Author should at least state whether the neglection of heat transport in nonfreeze–thaw period affect hydrological processes.

Reply: Thanks for the comment and suggestion. In the non-freeze–thaw period, all the water was in a liquid state, and the heat conduction had a minor effect on the water migration process. The model detected the non-freeze-thaw period based on ice content and temperature of each computing unit. In that period, only the moisture simulation was performed for simulation efficiency. We will make the supplementary modifications in Section 2.2 to address this point.

1. During the nonfreeze–thaw period, the soil hydrology was simulated with either a dualistic soil-gravel model or a Green-Ampt equation, and the selection of the two options depend on whether the rainfall intensity exceeded 20 mm/day (Line 677). Why was such threshold selected? Would the dualistic model more suit to the high intensity rainfall?

Reply: The runoff generation mechanism is different for the non-heavy and heavy rain scenarios: during non-heavy rain, there is saturation-excess, while during heavy rain, there is infiltration-excess. In heavy rain scenario, the Richards model is unstable for soil hydrology simulation while the Green-Ampt model is stable and has high computational efficiency. Therefore, in the WEP model, this threshold value was used to divide the simulations into two scenarios.

1. In the schematic figure shown in Fig. 3 (a), the author proposed a dualistic soil-gravel model, it is not clear whether the dualistic model is similar to the dual-porosity model proposed by Greke and van Genuchten (1993). Moreover, author should clearly state how to separate the water flow in such dualistic pore system.

Reply: Thanks for the insightful comment. In the dual-porosity model proposed by Greke and van Genuchten, the water transport medium is a mixture of soil and rock that is consistent from top to bottom; hence, the model generalized the medium into two systems: macropore and matrix pore. Our research object was the upper and lower layered medium with the thin soil layer and thick gravel layer in the Qinghai–Tibet Plateau. Based on this, we generalized the medium as an upper and lower dualistic soil–gravel structure to simulate the process of water and heat transport in different periods. In this model, the water flow is not separated like in the dual-porosity model. In future research, we will refer to the dual-porosity model to improve the mathematical description method of water transport in the gravel layer.

1. The soil water retention curve was described with van Genuchten model in Eq. B2 (Line 685), while the soil hydraulic conductivity function adopted a power function（B3）which is similar to that was used in Brooks-Corey model. Besides, the parameter n in Eq.B3 also adopted Mualem’s constant (Line 692). Such combination may be acceptable only if more cautions were taken for the parameterizations. Author should clarify why chosen to combine the selected soil water retention curve and soil hydraulic function, and how these soil hydraulic parameters were specified for distributed hydrological modeling.

Reply: We apologize for the unclearness on this part. The meanings of the two n values in Equation B2 and Equation B3 are different. In Equation B2, n and m are empirical parameters affecting the shape of the retention curve; m = 1–1/n. In Equation B3, n is Mualem’s constant. In the revised version, we will replace letters to clearly explain the meanings of different variables in the equation.

The combined application of the two models has been verified in the previous WEP COR model (Li 2019), which performs well in simulating water transport in frozen soil. These two models are mainly used to calculate the unsaturated soil hydraulic diffusivity D(θ_l) and the hydraulic conductivity K(θ_l) in Equation B1. Equation B1 was used to calculate the vertical movement of water in unsaturated soils.

1. In Page 20, Fig.8, why the simulated soil moisture differed between the two models in a freeze–thaw period (Line 414)? Modification in the proposed model may be solely focused on the nonfreeze–thaw period.

Reply: Differences in simulated soil moisture were caused by different model structures. The WEP-COR model did not consider the layered geological features of Qinghai-Tibet Plateau; the simulation object was homogeneous soil. Therefore, the simulated moisture of the WEP-COR model changed gradually in the vertical direction, and a large difference between simulated and measured values occurred below 40 cm (the soil layer thickness at the experimental site is 40 cm, with gravel layer below 40 cm.). The WEP-QTP model took this geological structure into account. Gravel layer has higher hydraulic conductivity and lower water-retention capacity, which is manifested in the simulated difference in the water content of the gravel layer. The simulated results from the improved WEP-QTP model were closer to the measured values.

The model improved in both the freeze-thaw period and non-freeze-thaw period. Like in the non-freeze-thaw period, the revised formula for water retention properties of the soil–gravel mixture was used to describe water retention curves for the lower gravel layer during the freeze-thaw period (Equation 1). The saturated hydraulic conductivity (K_s) of the soil or gravel layer was corrected by temperature (Equation 6). There are also improvements to the heat transport calculations; please see Section 2.2.2 for details.

1. All the figures have poor resolution. Please consider replacing all of them.

Reply: Thank you for your suggestion. In the revised version, we will replace all the figures with high-quality figures and improve figure layout.

Minor comments:

1. Line 132: “Temperature” should be “temperature”

Reply: Thanks, we will correct it and check for other potential errors.

1. Line 154: The “0” may be redundant.

Reply: Sorry for this mistake, we will correct it and check for other potential errors.

1. Line 164: The citation of MODIS data should be added.

Reply: The citation of MODIS data is: https://ladsweb.modaps.eosdis.nasa.gov/search/

1. Line 166: The citation China's second glacier inventory data set should be added.

Reply: Line 169 shows this citation (http://westdc.westgis.ac.cn/).

1. Line 167: add the citation of Water and Energy transfer Processes in Cold Regions (WEP-COR) model.

Reply: Thanks, we will add this citation in the revised version.

1. Line 267: The unit of saturated hydraulic conductivity Ks and snow water equivalent S should be consistent.

Reply: S is the daily variation of snow water equivalent (mm/d), which is consistent with the unit of daily precipitation. K_S is the saturated hydraulic conductivity, which was used for the calculation of water transport in Equation B1. Unit conversion was considered in the calculation process of the model.

1. Line 289: The unit of E，ρa, cp, and ra should be added.

Reply: We apologize for these omissions. We will supplement these contents as follows:  is the sum of the surface sublimation and evaporation rates (mm/day);  is the air density (kg/m³);  is the constant pressure specific heat of the air (MJ/kg/℃);  is the aerodynamic resistance (day/m).

1. Line 290: ra is aerodynamic resistance?

Reply: Yes, we will supplement its definition in the revised version.

1. Line 340: where can we found the calibrated parameters?

Reply: The calibrated parameters were in Line 344-354.

1. Line 454: Figure 10 In order to prove the conclusion in this paper that WEP-QTP can better simulate the measured runoff, it was suggested to plot the measured runoff data in the figure.

Reply: The simulated and measured runoff data are compared in Figures 5 and 9. Due to the limitations of the experimental site conditions, the hydrological cycle fluxes in Figure 10 have no measured values. Figure 10 was provided to compare the effect of model improvement on the runoff process. We will supplement this in the revised version.

1. Figure11 Legend and the scale is too small to read; It is recommended to mark the location of the three stations. What the source of the plotted data, measured snow thickness or the model simulation? If it is a map of measured snow, it is recommended to put it in the appendix. If it is a map of modeled results, suggest making a comparison with the actual measurement.

Reply: We apologize for the unclearness on this part. Figure 11 shows the snow thickness simulated by the model. The measured snow thickness was calibrated at the experimental site (Figure 6). We will modify this map as you suggested and supplement the site location.

References

Li J, Zhou Z, Wang H, et al. Development of WEP-COR model to simulate land surface water and energy budgets in a cold region [J]. Hydrology Research, 2019, 50(1): 99-116.

The manuscript “Application of a new distributed hydrological model based on soil-gravel structure in the Niyang River Basin, Qinghai-Tibet Plateau” applied a new model which considered the impact of gravel on water and heat transfer, as well as the snow cover in the study region. This work was necessary for the study region which has soils with large portion of gravel to a certain depth. However, when I read the manuscript, I felt that the work needs some major changes in order to make it clearer. The authors developed both the soil and snow processes, but they did not show how each process have improved the results. Besides, the description on the results are quite subjective and I did not see confident quantitative descriptions in multiple places, I will put the specific comments below. In general, I suggest a substantial revision to make this work more attractive and interesting to readers.

Dear Reviewer:

We appreciate the detailed and valuable comments, which have considerably improved the quality of our manuscript. Our responses to your comments are provided below.

Specific comments:

1. Line 63-70: I like this paragraph about the soil formation of QTP, but the position of this paragraph can be moved upward before introducing the gravel content impacts on soil heat transfer.

Reply: Thank you for this suggestion. We will revise this paragraph as you suggest.

1. Line 76-79: I think this sentence is repeated as it was already mentioned above.

Reply: Thanks for pointing this out, we'll revise this section by remove the repetition.

1. Line 79: how to adjust parameters? I think you mean calibration, but the calibration always depends on the function as the goal, i.e. soil temperature from surface or from lower layer, and/or soil water simulation accuracy.

Reply: Yes, it is parameter calibration. For the water and heat simulation of the underlying surface in the hydrological model, with the goal of improving the accuracy of soil water and soil temperature simulation, the hydrothermal parameters of the underlying surface were calibrated. In the revised version, we will clarify this part.

1. So I would suggest the authors to be more direct on demonstration of what you focus and why it is important. This would connect the whole storyline of this manuscript. Otherwise, I would not think you tell a good story on your work and the importance of it would be heavily lowered.

Reply: Thank you for your constructive feedback. We will modify the structure of the Introduction in the revised version, as you suggested, to help readers better understand the importance of our work.

1. Line 84: I did not see the impacts of rain intensity in your introduction.

Reply: We are sorry that we failed to make this clear enough in the previously submitted version. The runoff generation mechanism is different between the non-heavy and heavy rain scenarios. The non-heavy rain scenario indicates saturation-excess, and the heavy rain scenario indicates infiltration-excess. Therefore, in the WEP model, the runoff generation process was divided into two scenarios. In the revised version, we will clarify this part.

1. Line 143: to avoid

Reply: Thanks, we will correct it and check for other potential errors.

1. Line 154: what is “and0”?

Reply: The “0” was redundant. We will remove it and check for other potential errors.

1. Line 181-182: this is repeated.

Reply: Thanks for pointing this out, we'll remove the repetition.

1. Line 349-350: how were they determined?

Reply: These parameters were calibrated with actual measurement.

1. Line 420: Figure 8, again, interestingly, the model was good for frozen period because liquid water was very low, but during the thawing period, i.e. March and April, the model starts showing large discrepancy from observations, which I am curious about the possible causes.

Reply: The soil moisture simulation results of the WEP-QTP model in March and April were indeed not satisfactory, but there is still a clear advantage compared to the WEP-COR model. The lower accuracy in simulation can be attributed to soil inclusions within the gravel layers. Gravel layers are not homogeneous, but they were generalized in the model as a set of unified parameters; hence, there is a discrepancy in the model simulation results.

1. I just wonder the improvement is more about the hydraulic process than the thermal process? I can not agree that soil temperature was obviously improved from Figure 7 as the dynamics is still not well represented by the WEP-QTP model.

Reply: The improvement of the heat transfer process includes two aspects: one is to consider the thermal insulation effect of snow, and the other is to consider the thermodynamic parameters difference between gravel and soil. The addition of snow mainly reduces the temperature fluctuation in the surface layer, while the addition of gravel affects the heat conduction in the underlying layer. Impacted by the heterogeneity of the gravel layer, the improved model shows a discrepancy from the observations, but it is closer to the observations than the WEP-COR model.

1. Section 3.3: these are not quantitative and a little subjective. I would like to see how the new model is better and with some quantitative results.

Reply: Thanks for the comment and suggestion. In the revised version, we will add the discussion of the quantitative results in this section 3.3 to help readers understand the model improvements.

1. Line 457: I am wondering, how you compare the new model with the old model if it was changed in different processes, as you have developed the gravel related processes, and now you have a new snow model. I will be confused on what caused the differences in results.

Reply: By analyzing the sensitivity of snow and gravel layers to temperature simulation, we found that when only considering the influence of the thermodynamic parameters of gravel, the average RE was -4.38%, which is close to the average RE of WEP-COR (-3.6%). When only the insulation effect of snow was considered, the average RE was 0.70%, which is close to the average RE of WEP-QTP (0.08%). Therefore, snow is the main factor affecting the temporal and spatial variation of permafrost.

1. Conclusions: I am not so confident on the results for frozen period in the whole manuscript as I mentioned above, the snow and soil processes were both developed and I did not see the impact from one single process. I therefore additional sensitivity work should be done to quantitatively show how each process affects the results.

Reply: Thank you for your professional comment. As mentioned in the reply to question 13, we will add additional sensitivity analysis in the revised version to quantitatively show how each process impacts the simulation results. In addition, we will focus on revising the discussion section based on your comments above to enhance the persuasiveness of our conclusions.

General comments

To tackle the question of the cryo-hydrology of Tibetan catchments under climate changes, the authors use an already publish cryo-hydrological model and improve it in two ways: they use a stratigraphy of a soil lying on a soil with more gravels that they identify as widespread across the QTP and they account for the yearly cycles of growth and melt of a snowpack both thermally and hydrologically.

I think overall that the study is interesting because it contributes to give visibility of catchment scale cryo-hydrological modeling which is a key approach, currently under development, to understand how climate change will impact the water cycle in high mountain regions. I also think that the field observations on the stratigraphy give an important added value to the paper and they should be presented with more details to better assess the characterization of this stratigraphy and its spatial distribution. The study also conveys the interesting message and demonstrates that the stratigraphy is important both regarding the hydrology and the thermal behavior of the model setup. The importance of representing the snowpack is also interesting but probably more obvious, as you cannot really model realistically the hydrology and thermal regime of a catchment with a significant seasonal snowpack without accounting for it.

So for me, this is where the added value of the article is: interesting field observations that motivate an interesting sensitivity test on stratigraphy. Even though visible, the improvements on the model outputs are not stunning on the provided graphs but I believe that if the authors could provide observations or something we could consider “reality” on figures 10 and 12 this feeling could be improved. I realize that this study represents an important amount of work, the objective is clear and I think the global structure of the article is relevant.

Yet the manuscript has important flaws that needs major modifications. I detail that below.

Dear Reviewer:

We appreciate the detailed and valuable comments, which we addressed to substantially improve the quality of our manuscript. Our detailed responses to your comments are provided below.

1. In the first place, according to me, the idea that the whole QTP, a region of millions of km2 presents a more or less uniform stratigraphy with a layer of soil on top of a layer of gravels makes very little sense to me. I had a look at the studied catchment on Google Earth and I saw steep rock walls, colluvium, torrential streams, glacial valleys, all types of moraines, alluvial fans, braided fluvial systems… Meaning I saw the normal variety of landscape processes I could have expected and there is no reason they all produce this uniform stratigraphy at the scale of the catchment, not counting the whole QTP ! Not to mention that most of the sampling points are in the low lying parts of each valley. So the authors need to present better their field observations of the stratigraphy and discuss them in a much more cautious way.

Reply: There is indeed a variety of landscapes in the catchment, as you can see on Google Earth. The WEP-COR model generalizes this and divides the underlying surface into five classes: water body, soil-vegetation, irrigated farmland, non-irrigated farmland, and impervious area (Lines 655-Line 658). The water body class includes rivers, lakes and glaciers. The soil–vegetation class includes bare land, grassland and woodland. Impervious area consists of urban buildings and impervious surfaces. These five classes of underlying surfaces represent all landscapes of the catchment in the model. Simulation of soil water and heat transport is not applicable for water body and impervious area classes, while for other underlying surfaces, the improvement of the water and heat transport model was achieved in this study. Apologies for the lack of clarity regarding the surface classification in the previous version of the paper; we will describe this in the revised manuscript.

1. L76-77 “The geological features of the QTP are generally thin soil layers above the thick gravel layers with clear boundaries between them.” -> Any study supporting this at a large scale?

Reply: Affected by geological structure and freeze-thaw cycles, the phenomenon of the gravel layer under the thin soil layer is prevalent in the Qinghai-Tibet Plateau (Sun, 1996; Yang et al., 2009; Chen et al., 2015). Moreover, we confirmed this phenomenon in the present study with field observations in the study area, as shown in the figure below.

1. I also think mentioning “gravel” might be misleading. I think the authors should find a name to describe this type of formation (I see Wang et al. 2013 uses “soil-gravel mixtures” which sounds much more informative to me) because I think most readers reading “gravel” won’t think about this unsorted slope deposits but rather to well sorted alluvial gravel formation that are highly conductive and then it might sound very counter intuitive. Especially when the author says “since gravel can neither conduct nor store water” and that it “hinders the movement of water”.

Reply: Thanks for the comment and suggestion. In the revised version, we will replace "gravel" with "soil-gravel mixture (SGM) layer".

1. Another big problem lies in the presentation of the models. Many points are unclear or don’t make sense. I detail this below. A good example of this is the method used to calculate snow melt, which is first said to be based on a snow depth threshold between contour bands rather than on climatic variables (even though it is later the case, when the author mentions a PDD method). Not mentioning the so-called “snow sliding”. I suspect the problem doesn’t lie in the model itself since the output looks good, but rather in the description. I am happy to be shown I am wrong if it is the case, but from what I read and understand, there are major problems in the identifications of the processes and their representation as they are described (under the hypothesis that what is coded is different and correct).

Reply: Sorry for not being clear enough about this part in the previous version of the manuscript. There are some problems in our description, which we addressed in detail in answering questions 10 and 11.

1. Finally, I started correcting the English but I am not a native speaker and the task here is too important for a scientific reviewer. So I would recommend having the manuscript proofread by a native speaker because I think many formulations can be improved (see the examples I give below for the beginning).

Reply: Sorry for the grammatical and stylistic errors, before the revised version is submitted, we will invite native English speakers to polish the language.

Key specific comments

1. I had a look at Wang 2013 and, unless I am mistaken, one empirical parameter is missing in equation 1, where (1 - Wgravel) should be (1 – B x Wgravel) (Wang writes it B). Is this a mistake or is it me who missed a transformation of the equation ?

Reply: Thanks for professional questions. This is a mistake, and we will correct it in the revised version.

1. Since the model uses different types of equations depending on a threshold for rain, what happens at the transition between normal and heavy rain ? What if during a rain event the threshold is crossed, how smooth is the model in this regard ? Are there some data processing methods to smooth a potentially sharp transition ?

Reply: Thanks for the insightful comment. The flow in the model was simulated by day. When the daily precipitation exceeded the threshold, the Green-Ampt model was used. A sharp transition in the flow process has not been detected in the current simulations, but your question is very insightful, and we will conduct a series of studies on it in our future research.

1. L242

“The large portion of gravel in the gravel layer causes the formation of macropores, which are connected to form a fast channel for transporting water during heavy rains”

I am totally confused here. As I said earlier, it is indeed the general way gravel formations are treated (as a conductive layer). But before the authors wrote :

“However, since gravel can neither conduct nor store water, the gravel […] hinders the movement of water and affects the water retention curves”

So how does all of this work together? The authors have been, since the beginning, using studies on soils containing rock fragments to support a certain type of behavior from their bottom layer and now they argue towards another behavior because they have been calling these fragments “gravels”. And at this stage I am confused. Maybe I missed something here but if so, there is a lack of pedagogy/clarity in the way this model and stratigraphy works together.

Reply: Sorry for the unclearness on this part. There are macropores in the soil-gravel mixture (SGM) layer. In a saturated state, the macropores form a fast channel for transporting water. However, when the SGM layer is in an unsaturated state, the water mainly moves under the actions of the matrix potential and gravitational potential. Thus, in an unsaturated state, the macropores do not work, and the gravel will hinder the movement of water. In the revised version, we explain this in detail.

1. L255

“(the contour bands)”

The authors need to explain more clearly how the model works in the main text. I had to read the appendix to get a clearer idea of how this works. It is an unusual approach so it needs to be commented on. What decides the shape and extension of a band ? And it also needs some statistics: How many bands ? Average size of a band ? Average elevation range within a band ? I would also like to see a map with all these bands to see how all this looks like. Otherwise, what is done here remains very abstract.

Reply: We thank you for your valuable suggestion. Each sub-basin was divided into 1~10 contour bands as the basic calculation unit according to the elevation. The basin was divided into 871 contour bands. The average area of the contour bands is 20 km². The sub-basins of the watershed and the division of the contour bands are shown in the figure below.

1. L256-259

“When the snow thickness difference between two calculation units exceeded this threshold, snow meltdown occurred. The snow in the higher-altitude calculation unit slides into the next unit until the two units had the same snow thickness. The daily variation of snow water equivalent was calculated as follows”

Why is melt based on thickness difference ? Melt should be based on the climate input. All this makes very little sense, but I suspect these problems lie in the model description and not in the model itself.

Reply: Sorry for the unclear description. We will correct it in the revised version as follows:

"When the snow thickness difference between two adjacent contour bands exceeds this threshold, snow begins to slide between those contour bands. The snow in the higher-altitude contour band slides into the lower band until the two bands equalize in snow thickness."

1. L262-263

“when the difference in snow thickness between contour bands in the same sub-basin exceeds the threshold, the snow slides downwards until the snow thickness”

What is snow sliding ? I never heard of that and found nothing relevant on the net. The 2 important redistribution mechanisms I can think of are wind drift and avalanche. Snow creep also exists but is marginal in comparison. So what is the process here ?

Reply: In our study, "snow sliding" means avalanche, i.e. snow collapse driven by elevation difference. In the model, we generalized the avalanche as the redistribution of snow between two adjacent contour bands. In the revised version, we will replace the terms with "avalanche".

1. L290-292

“G is the heat flux (MJ/m2/d) conducted into the snow or soil, which was determined by the temperature difference between the soil or snow and the atmosphere near the surface. The above equation was combined with the ground heat conduction and energy balance equations”

I think this is wrong. G is the energy input in the ground that is used to drive heat conduction after the surface energy balance equation has been applied. So G is not derived from the atmosphere temperature near the surface, but H is. G is what you get when you sum the energy fluxes from the radiations and turbulent fluxes. Another problem is that the end of the sentence talks about heat conduction and energy balance equations. Conduction has not been introduced but energy balance is actually equation 7.

Reply: Sorry for the inaccurate description of this process. This part introduced the energy balance equation used in the model and the calculation method of each energy flux. G was calculated from the temperature difference between the underlying surface and the atmosphere near the surface and was the heat flux conducted into the underlying surface. The heat conduction in the underlying surface was only related to G. After it was determined, the heat flux and temperature of each layer were calculated via Equation 10. In the revised version, we will simplify this part, remove extraneous equations, and provide the calculation method for G.

1. L294

“We […] simplified the calculation by solving the H according to the energy balance equation after calculating the LE”

Well that is not what the authors say before. Equation 9 is clearly a way to calculate H from temperature inputs. It is impossible to deduce H from the energy balance equation because you deduce G from this equation knowing all the other terms.

Reply: Sorry for not being clear in the previous version of the manuscript on this part pointed out by you. As we replied in question 12, we will supplement the calculation method of G as follows:

G=C_Vu*d_u*(T_a-T_u)

where,C_Vu is the volumetric heat capacity of the underlying surface (MJ/m³/°C); d_u is the depth of the underlying surface affected by heat conduction (m); T_a is the air temperature on the day of simulation (°C); and T_u is the surface temperature of the underlying surface on the day before simulation (°C).

1. L297

“The temperature difference between the atmosphere and the surface is the source of heat conduction”

Why say this after calculating the surface energy balance ? The surface energy balance enables to calculate the energy change of the top cell, to work with temperature, the authors can then do ΔT = ΔE/Cp. Saying what I quote here after detailing an SEB module is more than confusing.

I don’t have this level of problematic issues with the rest of the paper. Yet I think that in general, the text of the result and discussion section could be lighter and more concise.

Reply: Thank you for the insightful comments. For clarification of the energy calculation part, you can refer to our replies to questions 12 and 13 here. Sorry for being not clear enough in the previous version of the manuscript on the aspects pointed out by you, which will be clarified in the revised version.

Specific comments

1. I don’t know where to put it so I write it here: to be able to understand what the new stratigraphy brings we need to have access to the WEP-COR stratigraphy, on Figure 3 for example.

Reply: The WEP-COR stratigraphy was shown in the previous version. See Figure B1 in the appendix for details.

1. L16

“The Qinghai-Tibet Plateau has a thin soil layer on top of a thick gravel layer”

I have 2 problems with this abstract opening:

Problem 1: See my previous comments, this cannot be true at the scale of a region as large as the QTP where one can find mountain peaks, peatlands, moraines, alluvial fans, blocky terrain… I suggest writing something like “For hydrological purposes, simplifying the representation of the QTP subsurface conditions to a thin soil layer on top of a thick gravel layer…” but this needs to be either demonstrated in a previous paper or in the present paper.

Reply: Thank you for the insightful comments. For the proof of the geological structure in the Qinghai-Tibet Plateau and the generalization method of the different underlying surfaces in the model, you can refer to the answers to questions 1 and 2.

1. Problem 2: I guess this is just a matter of personal preference, but I would recommend to start the abstract with a bit of context on what big question this study works with. Hydrology in mountainous cold regions and climate change…

Reply: Thank you for your professional suggestion. In the revised version, we will revise the abstract according to your recommendation.

1. L41-42

“plays an important role in ensuring the security of water resources in China and Southeast Asia”

Needs to be supported by a reference.

Reply: We added more references as follows:

[1] Liu X, Yang W, Zhao H, et al. Effects of the freeze-thaw cycle on potential evapotranspiration in the permafrost regions of the Qinghai-Tibet Plateau, China[J]. Science of the Total Environment, 2019, 687: 257-266.

[2] Yu L, Feng C Y. Recent progress in climate change over Tibetan Plateau[J]. Plateau and Mountain Meteorology Research, 2012, 32(3): 84-88.

1. L43-44

“cannot be ignored”

Needs also a reference. The sentence is also surprising. The authors could start the sentence by “The extensive glacier…” list the items and end the sentence with “have a major impact on the water cycle…”

Reply: Thanks for your professional suggestion. We will modify this sentence as you suggested:

“The extensive glacier, snow cover, and permanent and seasonal frozen soil have a major impact on the water cycle.”

Additional references are as follows:

[1] Yongjian D, Shiqiang Z, Jinkui W U, et al. Recent progress on studies on cryospheric hydrological processes changes in China[J]. Advances in Water Science, 2020, 31(5): 690-702.

[2] Zhiwei L I, Guo'an Y U, Mengzhen X U, et al. Progress in studies on river morphodynamics in Qinghai-Tibet Plateau[J]. Advances in Water Science, 2016, 27(4): 617-628.

1. L44-45

“On the surface of seasonal frozen soil and permafrost, seasonal thaw layers alternately freeze and thaw as seasons change.”

This is a convoluted way to say that both in permafrost and permafrost free areas, the ground undergoes seasonal freezing.

Reply: Thank you for pointing out this deficiency. We will modify it according to your suggestion.

1. L63-67

My expertise on the topic is limited but this section on the links between tectonics, sedimentology and granulometry of the Quaternary sediments could be better phrased and states obvious things that don’t show particular relevance for the study. I don’t understand the message the authors want to convey that is important for the paper.

Reply: Thank you for the comments. This paragraph introduces geological factors affecting the formation of the special underlying surface structure of QTP. In the revised version, we will simplify the description of this part.

1. L105

“8 °C”

Add the average elevation associated to this mean temperature

Reply: The average elevation of the catchment is 4688.6 m and this will be added to the revised version.

1. L111-112

“Permafrost accounts for approximately 23.65%, mainly distributed in the upper reaches ofthe basin and the high-altitude areas on both sides of the mainstream.”

Reference for this value ? Also I doubt one can reach such a precision in the significant numbers of the percentage.

Reply: This value was calculated from ground temperature in ArcGIS according to the definition of permafrost: ground that remains at or below 0 °C for at least two consecutive years (Biskaborn et al., 2019; Dobinski W, 2011). We will reduce the significant digits to the whole number, 24% of the area under permafrost.

1. L113-114

“The annual average temperature of the experimental site is 5.28 ℃, which is a seasonally frozen soil area.”

It is the first time the authors mention this site, maybe introduce it first.

Reply: Thanks for the comments. In the revised version, we will incorporate this sentence into the introduction of experimental site in (Lines 114-118).

1. Caption of Fig 2

Indicate where this is located on the map of Fig. 1

Reply: Figure 2 was taken near the experimental site shown in Figure 1, and we will make notes under this figure in the revised version.

1. L172-173

“The volume of the glacier was calculated by the area-volume empirical formula (Grinsted, 2013; Radić and Hock, 2010).”

Does this give access to volume changes along time ?

Reply: Yes, glacier volume changed over time. The area of glaciers was obtained from four remote sensing images from 1994, 2003, 2009 and 2015. We linearly interpolated the glacier volume calculated from the area, making its temporal change as the model input.

1. L188-189

“According to the geological characteristics of the QTP, this study improved the hydrothermal simulation methods of the non-freeze–thaw period and the freeze–thaw period.”

This is a conclusion, it should not be part of the methods.

Reply: Thank you for this comment. We will delete this part from the Methods section in the revised version.

1. L193-196

“In the non-freeze–thaw period, the calculation object of water movement was defined as the dualistic soil–gravel structure (Fig. 3a). The upper layer is soil whose thickness and number of layers are determined by the location of the calculation unit; the thickness of the soil layer gradually decreases from the foot to the peak of the mountain. The lower layer is the gravel layer (mixed layer of soil and gravel).”

This is really hard to read/understand, rephrase, with examples and tangible elements.

Reply: We rewrote this sentence and illustrated it with Figure 3c:

“In the non-freeze–thaw period, the calculation object of water movement was defined as the dualistic soil–gravel structure. The upper layer was soil, whose thickness was determined by the location of the calculation unit and which gradually decreased from the foot to the peak of the mountain (Fig. 3c). The lower layer was the gravel layer (mixed layer of soil and gravel).”

Figure 3c: Snow–soil–gravel layered structure.

1. L233-234

“Until the water has the same potential energy in the soil and the gravel, the INF breaks through the critical surface, and then the infiltration rate stabilizes (Fig. 4).”

I don’t understand this part. First I am unsure that potential energy is the good terminology (i.e. potential energy of the water in a dam), I assume it is the pressure head. And if it is so, the Green-Ampt model does not calculate the pressure head, it calculates the volume of infiltrated water or the depth of the infiltration front. So I don’t understand this sentence. Maybe I did not understand the situation correctly but then please clarify this point.

Reply: Potential energy here refers to soil water potential, including solute potential (not considered in this study), matric, gravity, and pressure potentials. The Green-Ampt model was derived by combining the Darcy’s law with the continuity principle (Green and Ampt, 1911). The volume of infiltrated water, depth of the infiltration front, and capillary suction pressure were used to calculate the potential gradient in Darcy's law. The specific derivation process can be found in the references.

1. Figure 4: Cumulative infiltration process of the WEP-QTP model

I don’t understand this figure. Please give an explanation in the caption. Re-explain the letters. Why are there 2 dashed lines, are they different scenarios ? I see now that this part of the figure is modified from Jia et al. (2001). I think that it should be cited as a source element of the figure. Also, now that I found this image from Jia, I understand that what is represented are the successive wetting fronts. Yet what I don’t understand is why we see these dashed curves. In the Green-Ampt model, the wetting front is horizontal.

Reply: Sorry for the obscurity on this part. A dashed line represents the wetting front at a moment, and the dashed lines in the figure represent the wetting fronts at different times: t₁, t_m, and t_itf. The dashed line is the actual wetting front, but the Green-Ampt model equates it to a straight line separating the saturated soil above from the soil below. In the revised version we will redraw this figure as you suggested.

1. Eq2 and L233-241

Here again it is hard to understand what the authors are doing. Where does this equation come from? Classically, the infiltration rate tends towards K because F(t) is at the place of Fitf here. But Fitf is hard to understand as it is a finite quantity and not a variable (i.e. the cumulative infiltration when the front breaks through). Also why ksoil when working with the “gravel” layers ? And What are “error caused by the different soil moisture content of the soil above the interface”. I think this paragraph needs more pedagogy to avoid giving the feeling that the authors are doing their own cooking with some well-established equations.

Reply: Equation 2 was modified from Equation B4, which calculates the stable infiltration rate after the infiltration front penetrates the soil-gravel interface. Because the saturated hydraulic conductivity of the upper soil was the upper limit of the infiltration rate of the lower soil-gravel mixtures layer, K_soil was used. The detailed calculation process of the infiltration rate and the calculation method of each parameter can be found in Equations B4-B6 in the appendix. In the revised version, we will improve the presentation of this section to help readers better understand our work.

1. L253

“determined by the snow water equivalent and snow density”

What is the approach ? Constant values ? Values derived from climate forcing data?

Reply: The thickness of the snow layer (cm) was calculated as follows:

h_snow=S_snow*ρ_l/ρ_s/10

where S_snow is cumulative snow water equivalent (mm); ρ_l is the densities of water (10³ kg/m³); ρ_s is the densities of snow (kg/m³), which was derived from climate forcing data as shown in Equations 11.

1. Equation 7

“RN = LE + H + G”

A more comprehensive way to write it is G=… because it shows how what the authors derive from the climate forcing data is used as an input from the model. Also what about this equation when there is snow ? Is it also applied ?

Reply: Thank you for your professional suggestion. In the reply to questions 12 and 13, we simplified the introduction of the energy calculation and provided the calculation method for G.

1. Equation 10

How is the ice content linked to the temperature ? I assume the authors also use a soil freezing curve.

Reply: The relationship between the water and heat transport of the frozen soil is mainly manifested in the dynamic balance of the moisture content of the unfrozen water and the negative temperature of the soil, which is shown in Equation B14 in the Appendix.

1. L350

“the riverbed conductivity was approximately 5.184 m/d”

Where is this value used ? I feel some part of the model description is missing. From what I understood there is the soil and the gravel layers, now it seems there is a riverbed layer. I just looked for the word riverbed, it does not appear before part 3. I have the hardest time understanding how this model works and I am trying hard.

Reply: The model improvement in this study does not involve the water exchange process of the river channel, but the riverbed conductivity is a sensitive parameter of the model; hence, the parameter calibration results were presented here.

1. L350-351

“The thickness of the soil layer at the mountaintop, mountainside, and foot of the mountain was 0.4 m, 0.6 m, and 1.0 m, respectively.”

Was there an attempt to characterize the stratigraphy based on the topography/morphology ? If so it is not mentioned in the method. And it is more than necessary for the message this study wants to convey. So if this effort has been made, please explain it in the methods. Also as I said before, the stratigraphic observations should be presented in detail somewhere. They really contribute to the added-value of the paper.

Reply: Yes, as mentioned in the reply to question 28, we have redrawn the generalized structure of the model, and in the revised version, we will supplement this part accordingly.

1. Figure 11

The legend is so small, with the resolution I got for the figure (which is low) I cannot read it.

Reply: Thank you for your suggestion. In the revised version, we will replace all the figures with high-quality figures and improve figure layout.

1. L429-424

“There might have been a soil interlayer at cm, and the measured water content was between the simulated values of the WEP-QTP and WEP-COR. The average RE of WEP-COR was 33.74%, and that of WEP-QTP was smaller at -12.11%. WEP-QTP could reflect the influence of gravel on the vertical migration of water.”

This discussion is really hard to follow. How speculative is the existence of this interlayer ? What is the physical process that makes it a relevant hypothesis ? I think if this suggestion is important it needs to be explained in more detail.

Reply: Thank you for the insightful comment. The soil-gravel mixtures (SGM) layer under the topsoil is not homogeneous, but in the model, we assume homogeneity of this layer and use a uniform set of parameters to describe its water and heat properties. In the SGM layers, the observed value of water content in the 160-cm layer was smaller than that in other layers and between the simulated values of the WEP-QTP and WEP-COR; hence, we speculate that there may be a soil interlayer.

1. Figure 10 and L453-454

“was inconsistent with the measured value”

If possible the authors should add observations on figure 10. Since the main message of the study is to show the benefits of the new stratigraphy, these benefits are visible only when compared to observations.

Reply: Sorry for the misunderstanding of the description here. The simulated and measured runoff values are compared in Figures 5 and 9. Due to the limitations of the experimental site, the hydrological cycle fluxes in Figure 10 have no measured values. Figure 10 presents the effect of model improvement on the runoff process. We will revise this in the revised version.

1. L455

“In addition, snow cover significantly contributed to the inconsistencies between the temporal and spatial changes of the frozen soil in the two models, which in turn caused variations in the groundwater recharge and discharge process.”

What supports these results ? Can the author provide a number or a graph that supports this ?

Reply: By analyzing the sensitivity of snow and gravel layers to the temperature simulation, we found that only considering the influence of thermodynamic parameters of gravel, the average RE was -4.38%, which is close to the average RE of WEP-COR (-3.6%). When only the insulation effect of snow was considered, the average RE was 0.70%, which is close to the average RE of WEP-QTP (0.08%). Therefore, snow is the main factor affecting the temporal and spatial variation of permafrost.

1. L656

“water body, soil–vegetation, irrigated farmland, non-irrigated farmland, and impervious area”

What about mountain bare lands above the tree line ? Which type was used?

Reply: As we replied to question 1, "Soil-vegetation includes bare land, grassland, and woodland". Therefore, mountain bare lands above the tree line belong to soil–vegetation.

Technical corrections

1. L16

“while”

I don’t see the point of this “while”, it is not putting 2 ideas in parallel. It is possible to end the first sentence before the while and restart with “This unique…”

Reply: Thank you for your suggestion; we will change the text accordingly.

1. L17

“To investigate the mechanism of the underlying surface structure on the…”

The sentence reads weird. The structure does not have mechanisms. And “to investigate” misses something like “the effect”, “the consequences”… I would rephrase it to something like “To understand the impacts of this subsurface structure on the water cycle of QTP catchments…”

Reply: Thank you for your professional suggestion, we will revise it in the revised version.

1. L17

“hydrothermal migration”

Is this term correct ? I googled it and found papers about the motion of magma. Which makes more sense because heat can trigger density gradients and thus motion. I do not really see it with fresh water on the continent.

Reply: Thank you for pointing this out. In the revised version, we will uniformly change this term to "water and heat transport".

1. L22

“the single soil”

I think “single” is not useful here.

Reply: Thank you, we will change “single soil” to “uniform soil profile” in the revised version.

1. L23-24

“In the non-freeze–thaw period”

When no phase change occurs in the ground

Reply: Thank you for your suggestion; we will apply it in the revised version.

1. L35

“the "tailing" process after October”

the observed "tailing" process after October (it if is indeed what the authors mean)

Reply: Thank you for your suggestion, we will revise it in the revised version.

1. L40

“typical”

This world carries little meaning. I would write “major”

Reply: Thank you for your suggestion; we will change it to “major” in the revised version.

1. L45-46

“Almost all ecological, hydrological, soil, and biological activities in the soil in cold regions occur here, hence it has been the focus of hydrological research in cold regions”

This sentence could be improved because of the “Almost all”, “soil in the soil” etc. I would rephrase “This region has received a lot of attention regarding hydrological research because of the great variety of biological and physical processes occurring at the surface and subsurface…”.

Reply: Thank you for your suggestion; we will revise it according to your suggestion in the revised version.

1. L84-85

“during the non-freeze–thaw period”

Under fully thawed conditions

Reply: Thank you for your suggestion; “under fully thawed conditions” sounds clearer. We will make relevant changes.

1. L86

“develop a hydrothermal coupling method”

This is unclear according to me, I would say “develop a modeling framework representing coupled heat and water transfers in the ground”

Reply: Thank you for your suggestion, we will incorporate your suggestions in the revised version.

1. L102

“litharenitewith”

Reply: Sorry for this mistake; it should be "litharenite with" here. We will make relevant changes.

1. L144-143

“stations. Avoid”

Reply: Sorry for this mistake, we will revise it in the revised version.

1. L153-154

“from the foot to the peak of the mountain”

decrease with elevation. I guess it gets to 0 even before the peaks.

Reply: Yes, this situation exists. However, for the sake of simulation efficiency, the soil thickness in the model was only divided by the mountaintop, the mountainside, and the foot of the mountain.

1. L154

“and0”

Reply: Sorry for this mistake; the 0 here was mistyped and has been deleted.

1. L279

“the upper boundary energy”

“energy fluxes”

Reply: Thank you for your suggestion; we will revise it according to your suggestion in the revised version.

1. L279

“calculated by meteorological elements”

“calculated based on the climate forcing”

Reply: Thank you for your suggestion; we will make relevant changes.

1. L289

“specific heat”

“specific heat capacity”

Reply: We will modify the phrases according to your suggestions.

1. L293

“the iterative method”

There are a few of them, be more specific.

Reply: Thanks for the comment and suggestion, we will rewrite this section in the revised version.

1. L389-390

“the heat preservation effect of the snow”

I guess the authors refer to the insulation effect of snow

Reply: Yes, we will revise the section to use professional terminology.

1. L400

“the thermodynamic properties”

Are the author talking about the temperature ? I am confused.

Reply: Thermodynamic properties here refer to heat capacity and thermal conductivity, which change the layer temperature by affecting heat transfer.

1. L429-430

“the unstable water-holding capacity”

I don’t understand why it is unstable.

Reply: The soil-gravel mixture (SGM) layers under the soil are not homogeneous, but in the model, we assume it is homogeneous and use a uniform set of parameters to describe its hydrothermal properties. Compared with the uniform parameters of the SGM layers, the water-holding capacity of each layer varies.

References

[1] Biskaborn B K, Smith S L, Noetzli J, et al. Permafrost is warming at a global scale[J]. Nature Communications, 2019, 10(1): 1-11.

[2] Chen H, Nan Z, Zhao L, et al. Noah modelling of the permafrost distribution and characteristics in the West Kunlun area, Qinghai‐Tibet Plateau, China[J]. Permafrost and Periglacial Processes, 2015, 26(2): 160-174.

[3] Dobinski W. Permafrost[J]. Earth-Science Reviews, 2011, 108(3-4): 158-169.

[4] Green W H, Ampt G A. The flow of air and water through soils[J]. The Journal of Agricultural Science, 1911, 4: 1-24.

[5] Hedstrom N R, Pomeroy J W. Measurements and modelling of snow interception in the boreal forest[J]. Hydrological Processes, 1998, 12(10‐11): 1611-1625.

[6] Sun H L, The formation and evolution of the Qinghai-Tibet Plateau, China: Shanghai Science and Technology Press, 1996, ISBN: 9787532340231.

[7] Yang K, Chen Y Y, Qin J. Some practical notes on the land surface modeling in the Tibetan Plateau[J]. Hydrology and Earth System Sciences, 2009, 13(5): 687-701.