A simple regional snow hydrological process-based snow depth model and its application in the Upper Yangtze River Basin

Snow depth plays a significant role in the regional water balance, for which snowfall is usually determined by a fixed temperature threshold in regional snow research. This study developed a regional hydrological process-based snow depth model in the Upper Yangtze River Basin by using spatially distributed critical temperature data derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data and station data. Based on meteorological station and remotely sensed data, daily snow hydrological components from 1 August 2003 to 31 July 2015 were simulated. Results show that the simulated snow depth patterns agreed with those of the observed snow depth. The multi-year average of the starting date and ending date of snow duration was 16 October and 6 June, respectively, and that of maximum annual snowfall was approximately 455 mm, of which canopy interception comprised 10%, with a maximum value of 50 mm. The proportion of snow sublimation was less than 20%, which was contributed by interception sublimation (40%), snow surface sublimation (40%) and sublimation underneath the canopy (20%). The maximum annual snow sublimation was 29 mm. Snow melting was the primary snow consumption pathway, and approximately 70% of snowfall melted. This research is significant for the assessment and management of water resources in this region. doi: 10.2166/nh.2019.079 om http://iwaponline.com/hr/article-pdf/50/2/672/549380/nh0500672.pdf 2020 Yanqun Ren Suxia Liu (corresponding author) Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China E-mail: liusx@igsnrr.ac.cn Yanqun Ren University of Chinese Academy of Sciences, Beijing 110049, China Suxia Liu College of Environment and Resources/SinoDanish Center, University of Chinese Academy of Sciences, Beijing 110049, China


INTRODUCTION
Global snow cover, which can influence climate change, energy balance and water cycles, reaches approximately 46 million km 2 in winter (Koivusalo & Heikinheimo ; Zhang et al. ). The amount of snow, which can be represented by snow depth, is an indicator of surface radiation and water budget, is an important input of hydrological models focused on runoff prediction, and affects the magnitude of the main source of water to many streams in spring and early summer (Clark et al. ).
Therefore, snow depth is crucial for human society, especially in mid-latitude regions, such as the Upper Yangtze River Basin (UYRB), which account for most of the snow in the world.
Studies of snow depth calculation, at point to regional scales and with various ground observation methods, mathematical statistics, and a focus on geodesy and physical mechanisms, have been conducted in the past few decades (Strack et al. ; Molotch & Margulis ; Revuelto et al. ). Since the early 1960s, remotely sensed data remote sensing is difficult to use for snow depth retrieval and is strongly influenced by clouds (Dai et al. ), but microwave (especially passive microwave) remote sensing has the advantage of strong penetrability and can monitor spatial and temporal variations of snow depth at both global and regional scales (Chang et al. ). However, it can be affected by snow characteristic parameters, such as snow grain and density, which can create uncertainty in the extracted snow depth. This is particularly prominent in plateau zones with thin and/or discontinuous snow cover distribution (Clark et al. ).
Snow hydrological processes are complex, and technology cannot clearly represent changes in these processes.  (Bowling et al. ). However, the CRHM has a more complete range of processes for this environment (e.g. blowing snow, snow interception, energy balance, and infiltration). The CRHM is adapted to cold regions where snow is characterised as dry, but the UYRB houses frigid, temperate and subtropical zones, and the snow processes in this area are more complex than those of other cold regions. Moreover, the type of snow is not homogeneous and varies with temperature.
Among these snow hydrological processes, the snowfall calculation is the first and most important step of the model operation. Precipitation separation is the core of the snowfall processes, as it determines the amount of solid precipitation and changes with time and terrain (Berghuijs et al. ). Rainfall and snowfall are two common precipitation patterns, which are different in the hydrological cycle (Harder & Pomeroy ). Rainfall can infiltrate soil or form runoff immediately, but snowfall dynamics change depending on the energy budget of the land surface. In most hydrological models, the air temperature threshold is a fixed empirical value, which does not agree with the real environment. A distributed value is more likely to reflect reality. Research on precipitation-phase partitioning has been conducted in the past few decades ( Most of the study area is covered by snow during the winter, and snowmelt in spring and early summer is one of the most important sources of discharge in the Yangtze River Basin. Therefore, a more accurate spatial distribution of snow depth significantly influences runoff calculation. Overall, this study aimed to build a distributed snow depth model to simulate components of the snow hydrological cycle during the snow season in the UYRB, China. The objectives of this study were to: (1) establish a more accurate simulation model of snow depth based on meteorological data (air temperature, land surface temperature (LST), precipitation, wind speed, relative humidity and sunshine duration), land type data (leaf area index (LAI) and normalised difference vegetation index (NDVI)) and earth parameters, which considers the canopy interception of snow and snow blowing; and (2) test the ability of the model built in this study to simulate components of the hydrological cycle and to evaluate the adaption of this model compared with different models or data resources.

MODEL DESCRIPTION
In this study, we established a distributed snow depth model for the UYRB based on hydrological processes adapted to this unique area by using some CRHM modules, including sublimation (Fang & Pomeroy ) and snow interception processes (Ellis et al. ). The model includes five modules, namely snowfall, canopy interception, snow sublimation, snowmelt, and snow accumulation modules.
Infiltration was not considered because low temperatures induce freezing of the soil. The main inputs included climate driving factors, land surface parameters, and earth parameters. The model structure and mass change of the hydrological processes are shown in Figure 1. Details of each module are described in the following sections.

Snowfall module
Research has shown that the single critical temperature method is simple, easy to use, and has had a high precision in many areas of China (Chen et al. ). To determine precipitation as snowfall or rainfall, a single air temperature threshold (T 0 ) was applied (Leavesley et al. ), where snowfall (P s ) is calculated as follows: where To weaken the influence of cloud pixels, a cloud removal algorithm with four steps of sensor combination, temporal combination, snowline determination, and seasonal filter (Dietz et al. ) was adopted in this study. By combining the sensors of MOD10A1 and MYD10A1, the cloud-cover percentage decreases to 55% (line MCD in Figure 2). After applying the seasonal filter, the cloud-cover percentage decreases significantly and finally reaches about 1% (line SF in Figure 2).
Using the improved snow cover dataset, the beginning date of snow duration (BSD) was estimated for each year.
The BSD refers to the day after which snow starts accumulating on the surface. To avoid including occasional snowfall, the BSD was calculated as the first day of the first continuous 3 days with snow cover. The BSD's temperature, which is the dividing point between the liquid and solid states of precipitation, was referred to as the 'critical temperature'. Furthermore, the spatial distribution of critical temperature for each year was computed and corrected by observed data. It is worth noting that only two precipitation patterns (snowfall and rainfall) were considered in this study.

Canopy interception module
The snowfall on forested land was comprised of two parts: the snowfall intercepted by the canopy and the snowfall falling directly on the ground across the canopy. The snowfall that was intercepted by the canopy was further divided into two parts: snowfall directly intercepted by the canopy, and snow that drops from the canopy. The mass balance equation of canopy interception is as follows: where M I is the amount of snow that remains on the canopy (mm), I is the amount of snow that the canopy intercepts from snowfall before unloading occurs (mm) and Ds is the amount of snowfall (mm) that unloads from the canopy to the ground. The hypotheses considered were that the interception efficiency decreases with canopy snow load and increases with canopy density and that snow unloading increases with time. I was calculated by the following equation (Hedstorm & Pomeroy ): where I max is the maximum intercepted snow load (mm) that can be retained by the forest canopy given the current canopy structure and temperature conditions, and I 0 is the intercepted snow load (mm) at the start of a snowfall event. Because the simulation period included several snow seasons, we set it to zero at the beginning of each snow season. Additionally, C c is the canopy cover, p is snowfall and c is the maximum ratio of snow to leaf contact area per unit area of ground (mm). In the canopy interception process, the maximum snowfall-holding capacity of the canopy leaves directly affects the amount of snow falling to the ground surface. The maximum canopy load I max can be calculated as a function of LAI, and the maximum snow load per unit branch area (S) (mm) (Schmidt & Gluns ).
The unit used for fresh snow density (ρ snow ) is mm.

Schmidt & Gluns () made extensive measurements
suggesting that the value of the maximum snow load per unit branch area (S p ) for spruce, which is the main forest type of the UYRB, is 5.9 mm. As ρ snow is not normally available in the meteorological record, an empirical relationship where T a is the air temperature (K) around the canopy.
The amount of snow that unloads from the canopy of canopy temperature and wind speed and was defined as follows:

Sublimation module
Snow sublimation is composed of pure snow surface sublimation and forest area sublimation. Pure snow surface sublimation (E s ) is often accompanied by blowing snow (MC blow ). The sublimation in the forest area (S f ) also consists of two parts: sublimation of canopy interception and sublimation underneath the canopy. To calculate the total sublimation accurately, all of these situations should be considered. Thus, snow sublimation (Ms) was calculated as follows: where the units of Ms, E s , MC blow , and S f is mm.

Pure snow surface sublimation
For sublimation of snow by heat absorption, the latent heat (Q LE ) was used to calculate the snow surface sublimation (E s ): where λ v is the latent heat of snow sublimation.

Blowing snow process modelling
Blowing snow consists of snow transport and sublimation.
Snow transport includes saltation transport and suspension transport. Sublimation rates are calculated for a column of saltating and suspended blowing snow extending to the top of the boundary layer. The calculation was performed as follows (Li et al. ): where MC blow is the blowing snow loss, Msalt (mm) is the snow transport in the saltation layer, Msusp is the snow transport in the suspension layer and EB (mm) is the blowing snow sublimation. The rate of mass transport (MT) (mm) by blowing snow across a unit width perpendicular to the wind includes suspension-layer and saltation-layer blowing snow: In place of the numerical integration used to evaluate Equation (8) where u 10 (m/s) is the wind speed at the height of 10 m, and Ta ( C) is the air temperature.
Analysis indicates that the variation of blowing snow flux with height is consistent with the Gamma distribution.
The blowing sublimation was calculated by the following as follows: where L S is the latent heat of sublimation, M is the molecular weight of water, R is the universal gas constant, λ T is the thermal conductivity of air, D is the diffusivity of water vapour in air, and ρ s is the water vapour saturation density.

Snow sublimation on forest land
The two components of snow sublimation on forest land are snow sublimation occurring on the canopy and snow sublimation occurring underneath the canopy, respectively.
The sublimation in the forest (S f ) (mm) is calculated as follows: where E I (mm) is the snow sublimation on the canopy and E F (mm) is the sublimation underneath the canopy. The formula of sublimation from canopy interception (Parviainen & Pomeroy ) (E I ) is as follows: where D is the diffusivity of water vapour in the air, ρ a is the air density near the canopy layer, sh is the Sherwood number, q 1 is the specific humidity of air near snow, q 2 is the specific humidity of snow surface and C is a scale coefficient (Parviainen & Pomeroy ): where T c is the canopy temperature, r is the radius of snow particles, u z is the windspeed near the canopy, v is the air viscosity factor, ρ i is the density of ice and k is the state factor of snow.
The snow sublimation underneath canopy can be calculated as follows (Ellis et al. ): where u avg is the mean daily wind speed and ea avg is the mean daily water vapour pressure.

Snowmelt module
The simulation of snowmelt is based on the energy balance principle, which is shown in Figure 3.
The energy balance method strictly follows the energy balance equation to calculate the snowmelt process. The energy budget components are net shortwave radiation (R ns ), net longwave radiation (L a À L t ), sensible heat (Q H ), latent heat (Q LE ), soil heat (Q G ) and energy input from precipitation (Q P ). The net energy (Q net ) and snowmelt (Sm) are calculated as follows: where λ v is the latent heat of snow sublimation.

Net shortwave radiation
The daily shortwave radiation is usually calculated via observation data of sunshine hours (Masahiko & Takanori ): where A is the albedo of the snow surface and Q s is the incident shortwave radiation.

Atmospheric longwave radiation
The function of atmospheric longwave radiation (L a ) is given as follows: where σ is the Stefan-Boltzmann constant, T a is the air temperature (K), and ε ac is the atmospheric emissivity (Campbell ).

Snow surface longwave radiation
The snow surface longwave radiation (L t ) is computed by the following function: where ε s is the snow emissivity and T s is the snow surface temperature (K).

Energy input from precipitation
The heat input of precipitation changes the energy status of the snow layer, and the energy input to the snow cover from precipitation (Q P ) can be expressed as follows (Walter et al.

):
Q P ¼ ρ w c w T a P r þ ρ w c s T a P s where ρ w is the water density, c w is the specific heat capacity of water and c s is the heat capacity of snow (these three parameters are all constant), P r is the rainfall, and P s is the snowfall.

Soil heat
In this study, the influence of soil heat (Q G ) in the period of snow accumulation was negligible, so this part of the energy budget was therefore not computed.

Sensible heat
Sensible heat flux (Q H ) refers to the energy exchange between the atmosphere and the underlying surface under the temperature change. The methods of calculation follow (Brutsaert ): where ρc p is the volumetric heat capacity of air and r as is the aerodynamic drag, which can be calculated as follows (Campbell ): where Z μ is the reference height of wind speed, x is the reference altitude of air temperature, Z m , Z h and k 1 are constants and μ 2 is the windspeed at the height of 2 m. The conversion method of wind speed (μ y ) at any height (y) is as follows (Allen et al. ): μ y ¼ μ 2 4:87 ln(67:8y À 5:42)

Latent heat
The latent heat (Q LE ) of the snow surface can be calculated as follows (Zhao et al. ): where λ v is the latent heat of sublimation, R d is the dry air constant, e(T a ) is the air vapour pressure, and e(T s ) is the vapour pressure of snow surface. We assumed that the snow surface vapour pressure is the saturated vapour pressure, e(T a ), presented as the product of the saturated vapour pressure under air temperature (e svp (T a )) and relative humidity (RH).

Snow accumulation module
By defining 'net snowfall' as the amount of snowfall after canopy interception sublimation, blowing snow, and melting snow are removed, snowmelt was calculated as follows: (1) If the net energy (Q net ) is negative, there is no snowmelt; (2) If Q net is positive and less than the melting energy for net snowfall, only part of the snowfall is melted; (3) If Q net is positive and greater than the melting energy for net snowfall, all snowfall is melted. Daily snow depth was calculated according to this principle. The snow season begins when snow cover persists for three consecutive days. The snow depth on day t (SD t ) was calculated as follows: where SD tÀ1 is the snow depth on day t-1; P st , M It , M St , and Sm t are the snowfall, canopy interception, snow sublimation, and snowmelt on day t, respectively.

Study area
The UYRB (Figure 4) with an area of about 40,000 km 2 (25-36 N, 90-105 E), is located in the southwest of the Tibetan Plateau. About 70% of the area is covered by seasonal snow, and the elevation of the area is 173-6,492 m. The climate of the UYRB ranges from frigid to temperate to subtropical zones and is impacted by warm oceanic airflow and the West Pacific subtropical high pressure. Summer precipitation is higher than at other areas of the same latitude.
In the eastern part of the basin, precipitation is abundant due to warm and humid currents from the Pacific and Indian Oceans, and the annual average precipitation ranges from 800 to 1,500 mm. The distribution of precipitation decreases from east to west and south to north. The type of precipitation in winter is snowfall, and the snowmelt runoff causes river discharge in spring.

Input data and process
Meteorological data

Parameters
The parameters used in the model are listed in Table 1.

RESULTS AND DISCUSSION
Critical temperature The multi-year average of BSD and critical temperature were calculated, as shown in Figure 5. Aside from the unusable observed records and station data, only a small number of stations had available snow depth records that could be used. Figure 6 illustrates that both the BSD and critical temperature showed significant spatial characteristics. The BSD decreased from north to south because of lack of snowfall.
The BSD of the south was close to 0, that of most of the north was less than 100, and that of the centre was between 100 and 140. These phenomena mainly occur because of the large elevation range and climate variation. Furthermore, the temperature decreased with higher elevation, which caused the early snowfall. The southern part of the UYRB mostly comprises a non-snow zone or has just a few days of snowfall, therefore, BSD was not present as in the snow covered areas. By analysing the applicability of these stations, only six stations observed snow depth data that could be used to validate the simulated results, as shown in Figure 6(a). Based on the multi-year average BSD, the critical temperature of the UYRB based on MODIS snow cover and station observed temperature data was obtained ( Figure 6(b)). The spatial pattern of critical temperature showed that critical temperature increased from north to south and that the value in the centre was larger than that in the north or south. The critical temperature ranged from 2.2 to 4.4 C. Under the same air temperature, the probability of snowfall increased with higher elevation, which agreed with the elevation change characteristic of decreasing from north to south. The characteristics of the  Figure 7 and Table 2.
The results indicate that the changing amplitude of the         capacity decreased as temperature decreased gradually.
Later on in the snow year, gradually increasing temperature caused snowmelt capacity to increase until the start of July, when it reached a maximum of approximately 100 mm.
The annual variation patterns of snow depth showed an inverted 'V' shape: beginning in October, reaching a peak value (36 mm) in February, and ending in June. In total, according to the simulation of all snow hydrology processes, the mass balance of snow in the UYRB was clearly described by our model.

Spatial patterns of snow components
The daily mean snow depth and its coefficient of variation are show in Figure 12.
The spatial distribution of snow depth, significantly affected by changes of climate and terrain, varied regionally. The snow depth in the centre of the study region was deeper than that in the north and south, with the south having the shallowest snow depth. The maximum mean daily snow depth was 46.5 mm. Due to the blocking of topography, the west cold trough lies on the west side of the plateau, which causes more snowfall there. The southern part of the plateau, where precipitation is abundant, is affected by the warm oceanic air from the Indian Ocean. Because the climate in this region is warm, it receives more rainfall.

CONCLUSIONS
In this study, a distributed snow depth simulation model based on snow hydrological processes was established to evaluate the snow hydrological components of the UYRB in China. This model consists of five modules, namely snowfall, canopy interception, snow sublimation, snowmelt, and snow accumulation modules, with climate driving factors, land surface parameters, and earth parameters as the main inputs. The validation results showed that the model simulated snow depth well, as it generally agreed with observed snow depth. However, the simulated value of snow depth was slightly larger than the observed snow depth. The change amplitude of the BSD varied from station to station in both simulated and observed data, but the differences between these results were negligible. Furthermore, the model built in the study can acquire the seasonal and multi-year evolution of snow depth and can accurately simulate peak values.
The simulated snow components showed that the sum of sublimation, including sublimation underneath the canopy (20%), interception sublimation (40%), snow surface sublimation (40%), was less than 20% of the total snowfall. The value of canopy interception was about 10%.
The main source of snow loss was snowmelt, as more than 70% of the snowfall in the study area was melted. The snowfall was mainly distributed in the north and centre of the study region, and the coefficient of variation was small in these regions. The snowmelt capacity was negatively correlated with altitude but positively correlated with temperature.
Although our model has improved snow depth simulation, some problems exist. One important influencing factor is limited data accuracy caused by a shortage and uneven distribution of meteorological stations in the study area. Moreover, the meteorological data obtained by interpolation using the GIDS method may have some errors due to the shortage of stations. This problem of data accuracy is inevitable in the current situation. With the simulated snow depth of the UYRB, the next step is to calculate the snowmelt runoff process. This research has significance in the assessment of water resources and water management in this region, such as the evaluation of the discharge of the Yangtze River Basin and its water resource management.