Description and application of a distributed hydrological model based on soil–gravel structure in the Qinghai–Tibet Plateau

The Qinghai–Tibet Plateau, known as the “Asian Water Tower”, has a thin soil layer with a thick gravel layer 15 underneath. Its unique geological structure, combined with widespread snow and frozen soil in this area, profoundly affect the water circulation processes of the entire region. To thoroughly study the water cycle mechanism of the Qinghai–Tibet Plateau, this study considered the geological and climatic characteristics of this area and selected the Niyang River Basin as the study area. The Water and Energy transfer Processes in the Qinghai–Tibet Plateau (WEP-QTP) model was constructed based on the original Water and Energy transfer Processes in Cold Regions (WEP-COR) model. This model divides the single soil structure 20 into two types of media: the soil layer and gravel layer. In the non-freeze–thaw period, two infiltration models based on the dualistic soil–gravel structure were developed based on the Richards equation in non-heavy rain periods and the multi-layer Green–Ampt model in heavy rain periods. During the freeze–thaw period, a hydrothermal coupling model based on the continuum of the snow–soil–gravel layer was constructed. This distributed hydrological model can dynamically simulate the changes in frozen soil and flow processes in this area. The addition of the gravel layer corrected the original model’s 25 overestimation of the moisture content of the soil layer below the surface soil and reduced the moisture content relative error (RE) from 33.74 % to −12.11 %. The addition of the snow layer not only reduces the temperature fluctuation of the surface soil, but also works with the gravel layers to revise the original model’s overestimation of the freeze–thaw speed of the frozen soil. The temperature RE was reduced from −3.60 % to 0.08 %. In the non-freeze–thaw period, the dualistic soil–gravel structure improved the regulation effect of groundwater on flow, stabilizing the flow process. The maximum RE at the flow 30 peak and valley decreased by 88.2 % and 21.3 %, respectively. In the freeze–thaw period, by considering the effect of the snow–soil–gravel layer continuum, the change in the frozen soil depth of WEP-QTP lags behind that of WEP-COR by approximately one month. There was more time for the river groundwater recharge, which better shows the “tailing” process after October. The flow simulated by the WEP-QTP model was more accurate and closer to the actual measurements, with https://doi.org/10.5194/hess-2021-538 Preprint. Discussion started: 8 November 2021 c © Author(s) 2021. CC BY 4.0 License.

average precipitation is affected by the Indian Ocean tropical ocean monsoon. Under the effect of the Indian low pressure, the southwest monsoon pushes a large amount of warm and humid air from the Bay of Bengal along the Yarlung Zangbo River Valley to the Niyang River Basin, causing precipitation in the basin with heavy rainfall and large vertical changes. The average annual precipitation in the basin is 1416 mm, and the average annual temperature is approximately 8 °C. Obvious temperature 100 changes occur from east to west with elevation.
In this study, a monitoring experiment of the coupling processes of water and heat during the seasonal freezing and thawing period was carried out on the mountainside of the Sejila Mountain in the lower reaches of the Niyang River Basin. The longitude and latitude of the study site are 94°21′45″ E and 29°27′12″ N, respectively, and the altitude is 4607 m a.s.l. The

Data description
The data required for this study were mainly divided into two categories: the first category was the data required for model construction (mainly including meteorology, geology and landform, terrain, soil type, land-use type, vegetation index, and https://doi.org/10.5194/hess-2021-538 Preprint. Discussion started: 8 November 2021 c Author(s) 2021. CC BY 4.0 License. glacier data), and the second category was the data used to verify the model results (including historical and experimental monitoring data). 115 Temperature, relative humidity, sunshine hours, and wind speed data were collected by the Nyingchi Meteorological Station in the basin and the Jiali Meteorological Station outside the basin, from 1961-2018. The data were obtained from the China Meteorological Data Network (http://data.cma.cn). In addition to the two meteorological stations in Nyingchi and Jiali, the rainfall data sources also included six rainfall stations (2013)(2014)(2015) including Gongbujiangda, Gengzhang, Baheqiao, Niqu, Kelaqu, and Zengba in the watershed and the contour map of annual precipitation in the Tibet Water Resources Bulletin (2012-120 2017). The temperature, relative humidity, sunshine hours and wind speed in the basin were interpolated from the meteorological station data by the reversed distance squared method with elevation correction. As for the precipitation data, the rainfall stations in the study area were concentrated in the valley (Fig. 1). If only these data were used to interpolate precipitation, there would be large errors in high altitude areas, affecting the accuracy of runoff simulation. Therefore, in this study, the precipitation-elevation relationship was first determined based on the contour map of annual precipitation in the 125 Tibet Water Resources Bulletin and the station precipitation data. Then, the daily precipitation data in the basin were obtained through elevation interpolation (Wang et al., 2017).
Due to the combined effects of plate tectonics, weathering, and erosion, a unique geological structure was formed in the Qinghai-Tibet Plateau with a thin soil layer on the top and a thick gravel layer on the bottom (Fig. 2). According to the geological characteristics of the Qinghai-Tibet Plateau, the Niyang River Basin, a first-level tributary of the Yarlung Zangbo 130 River, was selected in this study to represent a typical area. From the source of the river to the estuary, 24 sampling points at different altitudes were selected to conduct field investigations on the soil texture ( Fig. 1). Among them, points 1 to 16 were along the river, and points 17 to 24 were from the foot of the mountain to the peak. The soil thicknesses and compositions at the 24 sampling points were measured and analyzed. The soil layer of the Niyang River Basin is mainly sandy loam, with average sand, silt, and clay contents of 55.89 %, 31.2 %, and 12.91 %, respectively. The gravel layer is mainly round gravel, containing pebbles; the gravel content is approximately 50 %-65 %, the clay content 140 is 5 %-10 %, and the pores are filled with medium and fine-grained sand. The thickness of the soil layer gradually decreases from the foot to the peak of the mountain. It is approximately 40 cm on the hillside with higher altitude and increases to more than 100 cm in the valley.
The elevation data (Digital Elevation Model, DEM) used in this study were from the SRTM90 (Shuttle Radar Topography Mission), which is jointly measured by the National Aeronautics and Space Administration (NASA) and the National Imagery 145 and Mapping Agency (NIMA) with an accuracy of 90 m.
Soil type data were obtained from the second national soil census and "Chinese Soil Records". Land-use data came from the Resource Environment Science and Database Center, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (http://www.resdc.cn), and the data resolution was 30 m.
Moderate-resolution Imaging Spectroradiometer (MODIS) data from 2000-2017 were selected as the data source. Among 150 them, the leaf area index accuracy was 500 m, and the normalized difference vegetation index accuracy was 250 m; these were mainly used to calculate evaporation and vegetation interception processes, respectively.
The glacier data included China's second glacier inventory data set (1:100 000) and Landsat TM/ETM+/OLI remote sensing images. The second glacier inventory data comes from the China Cold and Arid Regions Science Data Center (http://westdc.westgis.ac.cn/). Landsat data came from the data sharing platform of the United States Geological Survey 155 (USGS) (http://glovis.usgs.gov/). ENVI software was used to extract glaciers, and the boundaries of the glaciers were finally determined with reference to Google Earth imagery. According to China's second glacier cataloging rules, the glaciers in the basin were classified and the glacier area was calculated. The volume of the glacier was calculated by the area-volume empirical formula (Grinsted, 2013;Radić and Hock, 2010).
Model verification data included historical data and experimental monitoring data. Historical data included daily measured 160 flow data from the Gongbujiangda Hydrological Station (2013-2016, 2018), the Baheqiao Hydrological Station (2013, and the Duobu station (2013-2018). The experimental monitoring data included the soil temperature and volumetric water content of the experimental site from 2016-2017. At the test point, a time-domain reflectometry sensor for monitoring the water content of the liquid water, a PT100 sensor for measuring the temperature, and a TensionMark sensor for measuring the potential of the substrate were installed every 10 cm in the vertical depth of the experimental pit at a depth of 1.6 m. Water, 165 heat, and potential energy were automatically monitored during freezing and thawing.

Introduction of WEP-COR
The WEP-QTP model was developed based on the Water and Energy transfer Processes in Cold Regions (WEP-COR) model.
For the sake of understanding and comparison, the WEP-COR model is briefly introduced in this section. The vertical structure 170 of WEP-COR is divided into the vegetation canopy or building interception layer, the surface depression storage layer, the aeration layer, the transition zone layer, and the groundwater layer. To accurately simulate the changes in soil moisture and heat from the surface to the deep layers and to reflect the influence of soil depth on the evaporation of bare soil and the water absorption and transpiration of vegetation roots, the aerated zone soil is divided into 11 layers (Fig. 3a). Among them, R surface represents runoff from the surface, and R i represents lateral flow or soil flow in the i-th layer of soil. R i is related to slope and 175 soil moisture content. E 1 represents soil evaporation, E r represents vegetation transpiration, Q i is gravity drainage of the i-th layer, P is precipitation, T a is atmospheric temperature, T i is the temperature of the i-th layer, and G i is the heat flux caused by the temperature difference between the i-th layer of soil and the adjacent soil layers. The thickness of the first and second layers was set to 10 cm, and the thickness of layers 3-11 was set to 20 cm.  The WEP-COR model divides the infiltration process into two scenarios for simulation: the heavy rain period and non-heavy rain period. The division standard is whether the daily rainfall exceeds 20 mm. In the non-heavy rain period, the Richards equation is used for daily scale simulation (Jia et al., 2009). During the heavy rain period, the multi-layer unsteady rainfall Green-Ampt model proposed by Jia and Tamai is used (Jia and Tamai, 1998). The relationship between the water and heat 190 transport of frozen soil is mainly manifested in the dynamic balance of the moisture content of the unfrozen water and the negative temperature of the soil. According to the principle of energy balance, the energy change of each layer in the freezethaw system was used for the soil temperature change and water phase change in the system.
In terms of the horizontal structure, the WEP-COR uses the contour bands inside small sub-basins as the basic calculation unit ( Fig. 3b), and it can fully consider the vertical changes of vegetation, soil, air temperature, precipitation, and other factors in 195 the basin with the elevation. Each unit is divided into five types according to the land-use type: water body, soil-vegetation, irrigated farmland, non-irrigated farmland, and impervious area. The calculation result of the water and heat flux in each type was weighted by area to obtain the water and heat flux of the contour band. Evapotranspiration of water and soil was calculated using the Penman formula, and the vegetation canopy evaporation was calculated using the Penman-Monteith formula (Monteith, 1973). The subsurface runoff was calculated based on slope and soil hydraulic conductivity, and groundwater 200 movement was calculated by the Boussinesq equation (Zaradny, 1993). The confluence of overland flow and channel flow were calculated using the kinematic wave method (Jia et al., 2001). The "degree-day factor method" (Hock, 1999)

Model improvement 205
Based on the WEP-COR model, this study developed the improved WEP-QTP (WEP-Qinghai-Tibet Plateau) model. In contrast to the general cold areas where the WEP-COR model is applied, the widespread dualistic soil-gravel structure in the Qinghai-Tibet Plateau has a great impact on the water cycle processes in the basin. According to the geological characteristics of the Qinghai-Tibet Plateau, this study improved the simulation methods of the non-freeze-thaw period and the freeze-thaw period. 210 In the non-freeze-thaw period, the calculation object of water movement above the groundwater was defined as the dualistic soil-gravel structure (Fig. 4a). The upper layer is soil, and its 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).

220
In non-heavy rain periods, as in the WEP-COR model, the water movement process was described by the one-dimensional vertical Richards equation as follows: where θ l is the volumetric content of liquid water in the soil or gravel layer (cm 3 /cm 3 ); D(θ l ) and K(θ l ) are the unsaturated soil 225 hydraulic diffusivity (cm 2 /s) and hydraulic conductivity (cm/s); and t and z are the time and space coordinates (positive vertically downward).
The Van Genuchten function (Van Genuchten, 1980) was used to describe the upper soil water retention curves: (2) https://doi.org/10.5194/hess-2021-538 Preprint. where θ s is the saturated water content (cm 3 /cm 3 ); θ r is the residual water content (cm 3 /cm 3 ); h is the matric suction (cm); α is 230 an empirical parameter (cm -1 ); n and m are empirical parameters affecting the shape of the retention curve; and m = 1-1/n. However, because gravel can neither conduct nor store water, the gravel, which accounts for 50 %-65 % of the gravel layer, hinders the movement of water and affects the water retention curves (Cousin et al., 2003). Therefore, the revised formula for water retention properties of the soil-gravel mixture was used to describe the lower gravel layer water retention curves : where A is an empirical parameter, λ is the pore-size distribution parameter (λ < 1), and ω gravel is the volume ratio of the gravel in the gravel layer.
In the heavy rain period, the multi-layer Green-Ampt equation was used to calculate the infiltration process when the infiltration front (INF) was in the soil layer (Fig. 5), which is the same as in WEP-COR. When the INF reached the m-th layer 240 of soil, the soil infiltration capacity was calculated by the following formulas: where f is the infiltration capacity (mm/h); A i-1 is the total water capacity of the soil above the i layer (mm); B i-1 is the error 245 caused by the different soil moisture content of the soil above the i layer (mm); F is the cumulative infiltration (mm); k i is the hydraulic conductivity of the i-th soil layer (mm/h); L i is the soil thickness of the i-th layer (mm); SW m is the capillary suction pressure at the INF of the m-th layer (mm); and ∆ = − . The cumulative infiltration quantity F when the INF reaches the m-th layer is calculated based on whether there is water accumulation on the ground surface. If the ground surface has accumulated water when the INF reaches the m-1th layer, Eq.
(7) was used; otherwise, Eq. (8) was used (Jia and Tamai, 1998): 255 where t is the time; t m-1 is the time when the INF reaches the interface between the m−1 and m layers; t p is the start time of the water accumulation; F p is the cumulative infiltration quantity at t p ; and I p is the precipitation intensity at t p .
When the INF moves to the interface of the soil and gravel (Layer itf), the front movement slows down because the water suction of the gravel layer is less than that of the soil (Mao and Shang, 2010). 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. 5). Therefore, 265 a new multi-layer Green-Ampt model based on the soil-gravel structure was proposed as follows: where f gravel is the stable infiltration rate after breaking through Layer itf; k soil is the saturated hydraulic conductivity of the soil layer (mm/h); A itf is the total water capacity of the soil above the interface (mm); B itf is the error caused by the different soil moisture content of the soil above the interface (mm); and F itf is the cumulative infiltration when the front breaks through 270 Layer itf (mm).
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 (Fig. 4a). After the INF breaks through the interface, the infiltration water preferentially recharges the groundwater through the macropores. At this time, the accumulated infiltration quantity is as follows: 275 where Q gd is the quantity of groundwater recharge by infiltration, = ( − ), and t itf is the time when the INF breaks through the interface.
During the freeze-thaw period, in addition to considering the impact of gravel on hydrothermal transfer, the higher reflectivity of snow to shortwave solar radiation and its contribution to heat insulation were also taken into consideration. The 280 hydrothermal coupling simulation object was defined as the snow-soil-gravel layer continuum (Fig. 4b). A snow layer was added on top of the dualistic soil-gravel structure, the thickness of which was determined by the snow water equivalent and snow density. For the heat transfer process, assuming that the upper boundary of the system is the atmosphere, which controls the input and output of the system energy. When there is snow on the surface, the atmosphere first exchanges energy with the snow layer, 285 and then the snow layer exchanges energy with the soil. When there is no snow, the atmosphere directly exchanges energy with the soil, and the upper boundary energy can be calculated by meteorological elements. The lower boundary at the bottom is the transition layer or groundwater layer, assuming that it maintains a constant temperature.
The energy balance equation of the surface can be expressed by the following equation (Jia et al., 2001): where RN is the net radiation flux (MJ/m 2 /d); LE is the latent heat flux (MJ/m 2 /d), which was calculated from the melting and evaporation of the snow layer and the freezing, thawing, and evapotranspiration of the soil layer moisture; H is the sensible heat flux (MJ/m 2 /d), which was obtained from the remainder of the surface energy balance equation; and G is the heat flux (MJ/m 2 /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. 295 The surface snow or soil temperature (when there was no snow cover) was calculated by the forced recovery method (Douville et al., 1995;Pitman et al., 1991). For soil and gravel layers, the average temperature was represented by the temperature in the middle of the layer. The temperature difference between the atmosphere and the surface is the source of heat conduction; after the surface temperature was determined, the heat flux and temperature of each layer were calculated by the following formula (Shang et al., 1997;Wang et al., 2014): 300 where C v and λ are the volumetric heat capacity (J/[m 3 ·°C]) and thermal conductivity (W/[m·°C]) of the soil or gravel layer, respectively; L f is the latent heat of ice melting (3.35 × 10 6 J/kg); T is the temperature (°C) of the soil or gravel layer; ρ I is the ice density (kg/m 3 ); θ I is the volumetric content of ice in the soil or gravel layer (cm 3 /cm 3 ); and z is the layer thickness (m).
The redistribution of snow on the Qinghai-Tibet Plateau is affected more by the elevation difference than by wind. The daily 305 variation of snow water equivalent was calculated as follows: where S is the daily variation of snow water equivalent (mm/d); S p is the snow water equivalent from precipitation (mm/d), S p is equal to the daily precipitation when the average temperature of the day was T a < 2 °C, otherwise S p = 0; S d is the snow water equivalent variation due to snow sliding down (mm/d), when the difference in snow thickness between contour bands in the 310 same sub-basin exceeds the threshold, the snow slides downwards until the snow thickness is the same; S m is the quantity of snow melting equivalent (mm/d), which was calculated by the degree-day factor method (Hock, 1999) as follows: where d f is the degree-day factor (mm/[°C·day], i.e., 4 mm/[°C·day] in this study); T S is the critical temperature of snow melting (°C, i.e., −1 °C in this study), assuming snow melt starts when T a is > T S .
During the freeze-thaw period, only liquid water migrates. The soil vertical heat flux transfer can be written as follows (Shang et al., 1997;Wang et al., 2014): where ρ l is the water density (kg/m 3 ).
Temperature is the driving force of the water phase change. The relationship between the water and heat transport of frozen 320 soil is mainly manifested in the dynamic balance of the moisture content of the unfrozen water and the negative temperature of the soil: where θ m (T) is the maximum unfrozen water moisture content corresponding to a negative soil temperature.
For the snow layer added to improve the model, the main hydrothermal parameters include thermal conductivity, volumetric 325 heat capacity, and snow density. The calculation formulas of each parameter are as follows: The thermal conductivity and volumetric heat capacity of snow were calculated as follows (Goodrich, 1982;Ling and Zhang, 2006 where ρ s is the snow density (kg/m 3 ); T a is the atmospheric temperature (°C); λ s is the thermal conductivity of the snow (W/[m·°C]); and C Vs is the volumetric heat capacity of the snow (J/[m 3 ·°C]).
As opposed to a single soil medium, the presence of gravel has a great influence on the hydrothermal transfer parameters of 335 the gravel layer. The main hydrothermal parameters of the soil-gravel layer include volumetric heat capacity, thermal conductivity, and soil hydraulic conductivity. The calculation formulas for each parameter were as follows: Volumetric heat capacity were calculated by Eq. (23) (Chen et al., 2008): where θ s , θ l , and θ I are the saturated volumetric water content, volumetric liquid water content, and volumetric ice content of 340 the soil or gravel layer, respectively; C s , C l , and C I are the volumetric heat capacity (J/[m 3 ·°C]) of the soil or gravel layer, water, and ice, respectively; at 0 °C, the soil and gravel layers have values of 1. The thermal conductivity calculation referred to the IBIS model, as follows (Foley et al., 1996) where λ and λ st are the actual thermal conductivity of the soil or gravel layer and the thermal conductivity in the dry state (W/[m·°C]), respectively, and ω gravel , ω sand , ω silt , and ω clay are the volume ratios of the gravel, sand, silt, and clay, respectively.
The hydraulic conductivity was calculated as follows (Chen et al., 2008): where θ r is the residual water content of the soil or gravel layer; K (θ l ) is the hydraulic conductivity (cm/s) of the soil or gravel layer when the liquid water content is θ l ; K s is the saturated hydraulic conductivity of the soil temperature correction (cm/s); and n is Mualem's constant.
K s can be calculated as follows (Chen et al., 2008;Jansson, 2004): where K is the initial saturated hydraulic conductivity (cm/s) and K 0 is the minimum hydraulic conductivity (cm/s) under freezing conditions. Considering the difference in the hydrodynamic properties of the soil and gravel layer, the K 0 value for soil was considered to be 0 cm/s. For the gravel layer, due to the larger pores, K 0 has a value of > 0 cm/s; T is the temperature of the soil or gravel layer (°C); and T f is the critical temperature (°C) corresponding to the minimum hydraulic conductivity.

Model evaluation criteria 360
Data from January 2013 to December 2018 were used to evaluate the simulation results of daily flow rates at Gongbujiangda, Baheqiao, and Duobu stations. The performance of the model was first evaluated using a qualitative assessment via graphs and then assessed quantitatively using statistical metrics including the Nash-Sutcliffe efficiency (NSE) and relative error (RE).
The NSE and RE were calculated as follows: where N is the number of observations; O i is the observed value; � is the mean observed value; and S i is the simulated value.

Model calibration and validation
We calibrated and verified the daily flow process of Gongbujiangda, Baheqiao, and Duobu stations-located upstream, on the 370 largest tributary, and downstream, respectively-from 2013 to 2018. The data from Duobu station were split into two parts: data from 2013 to 2015 were used for calibration and those from 2016 to 2018 for validation. The discontinuous, measured flow data from Gongbujiangda and Baheqiao stations from 2013 to 2018 were used to verify the model.
The parameters of the model were mainly divided into four categories: underlying surface parameters, vegetation parameters, soil parameters and aquifer parameters. All parameters have physical meaning and can be estimated based on observational 375 experimental data or remote sensing data. The sensitivity of the above four types of parameters was analyzed (Jia et al., 2006), and the sensitivity of these parameters was divided into three levels: high, medium, and low. Highly sensitive parameters included soil thickness, soil saturated hydraulic conductivity, and riverbed material permeability coefficient. The model was  water storage capacity and permeability of the aquifer were considerably improved owing to the addition of the gravel structure.
The WEP-QTP simulation flow process was smoother, and a large flow peak was not easily formed. In general, the WEP-QTP model delivered an acceptable performance for the Niyang River Basin and achieved efficiency coefficients of NSE > 0.75 and RE < 10 % for the validation period. The simulated flow was able to be used for further analysis.

Soil-gravel temperature
The soil-gravel temperature simulation results of the WEP-QTP and WEP-COR models are shown in Fig. 8. The soil thickness at the test site was approximately 40 cm. For the 10 cm layer, because the simulation parameters of the WEP-QTP and WEP-COR models were the same, the temperature simulation results were consistent except during the snow cover period (i.e., the 415 time period when the snow thickness was > 5 cm, 27 December 2016-20 March 2017). Due to the heat preservation effect of the snow, the heat transfer and temperature fluctuations of the surface soil were reduced. The RE of WEP-QTP in the 10 cm soil layer was 11.1 % during the snow cover period, which was much less than the 46.0 % of WEP-COR. The simulation results of the two models in the 20 cm layer did not differ considerably. For the 40 cm and 60 cm layers, the moisture content of the soil in the WEP-COR model was greater than that of the gravel layer in the WEP-QTP model. The higher heat capacity 420 of water reduces the thermodynamic difference between the soil and the gravel, resulting in a small temperature difference between the WEP-QTP and WEP-COR models in the early freezing stage. However, as the temperature decreased, the moisture in the soil was converted into ice with a smaller heat capacity; thus, the difference in thermodynamic properties between the gravel and the soil gradually increased. During this period, the simulation difference between the two models reached a maximum of 1.41 °C (40 cm layer, 26 January 2017). For the temperature simulation below 60 cm, because the temperature 425 was higher compared with that in the 40-60 cm layer, the thermodynamic properties of the gravel in the WEP-QTP model and the soil in the WEP-COR model are not significantly different due to water phase change, and the difference in the non-snow cover period was not as great as in the 40-60 cm layer. During the snow cover period, the effect of snow on the temperature was also reduced due to the weakening of the upper soil-gravel layers. In general, the snow cover reduced the heat transfer and temperature fluctuations of the soil layer, which improved the simulation accuracy of the surface soil temperature. In 430 addition, by neglecting the special hydraulic and thermodynamic properties of the gravel layer, the WEP-COR underestimated the soil temperature. The snow and gravel layer collectively resulted in the temperature simulation difference. The average RE of WEP-COR was −3.60 %, and that of WEP-QTP was 0.08 %. The results of the WEP-QTP simulation were closer to the actual measurement; thus, it was able to accurately reflect the temperature changes of each layer during freeze-thaw.  Figure 9 shows the comparison between the simulated and measured values of the liquid water content in the freeze-thaw 440 period of the WEP-QTP and WEP-COR models at the experimental points in 2016-2017. During the freezing period (December-March), the upper layer liquid water content first dropped because of the temperature drop, and then the lower gravel layer liquid water content dropped, stabilizing in January-February. When the temperature increased in March, the upper layer initially began to melt, resulting in an increase in the liquid water content and then the subsequent melting of the lower layer. After the thawing period, the upper part of the soil-gravel layer had a higher water content than that in the lower 445 part due to the infiltration of snow melt, and it was also higher than before freezing. During the entire freeze-thaw period, the simulation of the 10 cm soil layer was less affected by the gravel layer because the water-holding capacity of the 20-40 cm soil layer was greater than that of the underlying gravel layer. The moisture in the WEP-QTP was more easily retained in the upper soil, and the moisture content was higher when the soil began to freeze and the snow melted, which is closer to the actual measurement. Below 60 cm, the WEP-COR did not consider the impact of the 455 gravel layer, and the water content was higher than the measured value throughout the freeze-thaw period. However, due to the large uncertainty of the compositions of the soil and gravel layer, the unstable water-holding capacity of the soil-gravel layer cannot be accurately reflected when the model is generalized, which also leads to a certain difference between the WEP-QTP simulation and the measured values. There may have been a soil interlayer at 160 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 460 WEP-QTP was smaller at −12.11 %. WEP-QTP was able to reflect the influence of gravel on the vertical migration of water.

Simulation and comparison of watershed flow process
To further analyze the improved WEP-QTP compared with the original WEP-COR, three sites with data for 2014 were selected, and the daily flow data of the three stations were compared with the simulated data of the two models (Fig. 10). It can be seen from Fig. 10 that the simulation difference between the two models was mainly from June to November. The simulation 465 performance of the WEP-QTP model was better than that of the WEP-COR model in three aspects: the peak value of the flood season was not too high; the valley value of the flow process was higher than that of the WEP-COR model; and there was a tailing process after October.  Figure 11 shows a comparison and analysis of the changes in hydrological cycle flux, and Fig. 12 shows the percentage of frozen soil area in the basin. It can be seen from Fig. 11 that the runoff from rainfall of WEP-QTP was smaller than that of and flow peaks are not easily formed. This is why the WEP-QTP model performed better in the peak simulation during the flood season. In addition, the larger saturated hydraulic conductivity of the gravel layer also increased the groundwater recharge and discharge of the WEP-QTP model, so that it performed better in simulating low values of the flow process. Moreover, the temporal and spatial changes of frozen soil also affect the quantity of groundwater recharge to the river. Due to the effect of snow and gravel, the change in frozen soil depth of WEP-QTP lagged behind that of WEP-COR; the area of WEP-485 COR frozen soil with a depth greater than 1 m reached its maximum in January and its minimum in August. For WEP-QTP, the value reached its maximum in March and its minimum in September (Fig. 12). This resulted in the groundwater discharge peak of WEP-QTP being one month later than that of WEP-COR (Fig. 12), and there was more time for the river groundwater recharge, which shows a better tailing process after October. In general, the snow-soil-gravel layer structure changed the water circulation flux process in the basin. The peak value simulated by the WEP-COR model was large and the valley value was small, and there may no significant difference in the 495 monthly runoff process. However, in the daily flow process, taking Duobu station in 2014 as an example, compared with the measured value, the maximum difference between the RE of WEP-COR and WEP-QTP was 21.3 % (−45.0 % and −23.7 %, respectively, on 9 August 2014) at the valley value and 88.2 % (130.5 % and 42.3 %, respectively, on 15 August 2014) at the peak value. Ignoring the snow-soil-gravel layer structure greatly impacts the hydrological forecast, reservoir regulation, and water resource utilization. 500

Conclusions
This study combined the geological characteristics of the thin soil layer on the thick gravel layer and the climate characteristics of the long snow cover period in the Qinghai-Tibet Plateau. With the Niyang River Basin as the research area, the WEP-QTP model was constructed based on the original WEP-COR model. This model divides the single soil structure into two types of media: soil and gravel layers. In the non-freeze-thaw period, two infiltration models based on the dualistic soil-gravel structure 505 were developed based on the Richards equation in non-heavy rain periods and the multi-layer Green-Ampt model in heavy rain periods. During the freeze-thaw period, a hydrothermal coupling model based on the continuum of the snow-soil-gravel layer was constructed. This model was used to simulate the water cycle processes of the Niyang River Basin, and the improvement effect of the model was analyzed by comparison with the WEP-COR model.
Compared with the simulation results prior to improvement, it was found that the addition of snow not only reduces the surface 510 soil temperature fluctuations, but also interacts with the gravel layer to reduce the soil freezing and thawing speed. The low estimation of temperature by WEP-COR was corrected, and the RE was reduced from -3.60 % to 0.08 %. At the same time, the WEP-QTP model can reflect the impact of the gravel layer under the soil on the vertical movement of water and accurately describe the dynamic changes in moisture in the soil and gravel layers; the RE of the moisture content was reduced from 33.74 % to −12.11 %. 515 According to the comparison of the WEP-QTP simulation and measured results of the main stations in the Niyang River Basin, the daily flow process simulated by the model is in line with the actual situation, and the flow simulation result is more accurate (Nash > 0.75 and |RE| < 10 %), which is a considerable improvement compared with the WEP-COR model. In the non-freezethaw period, the dualistic soil-gravel structure increased the recharge and discharge of groundwater and improved the regulation effect of groundwater on flow, stabilizing the water flow process. The maximum RE at the flow peak and valley 520 decreased by 88.2 % and 21.3 %, respectively. In the freeze-thaw period, by considering the effect of the snow-soil-gravel layer continuum on soil freezing and thawing processes, the change in frozen soil depth of WEP-QTP lagged behind that of WEP-COR by approximately one month. There was more time for the river groundwater recharge, which shows a better tailing process after October. In contrast to the general cold area, the special geological structure and climatic characteristics of the Qinghai-Tibet Plateau 525 change the water cycle processes in the basin. Ignoring the influence of the dualistic soil-gravel structure greatly impacts the hydrological forecast and water resource assessment.

Data and code availability
The datasets and model code relevant to the current study are available from the corresponding author on reasonable request.

Author contribution 530
PW performed the model programming and simulations. ZZ, KW and YL conducted the field experiments. ZZ, JL, YL, and JL contributed to the model programming. The writing was performed by PW and ZZ. CX, YJ and HW contributed to the writing of the paper.

Competing interests
The authors declare that they have no conflict of interest. 535