Abstract

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.

INTRODUCTION

Global snow cover, which can influence climate change, energy balance and water cycles, reaches approximately 46 million km2 in winter (Koivusalo & Heikinheimo 1999; Zhang et al. 2014). 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. 2006). 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. 2004; Molotch & Margulis 2008; Revuelto et al. 2015). Since the early 1960s, remotely sensed data have been used in snow research, providing a new method for resolving spatial problems, and great progress has been made in this field (Hall et al. 1995; Riggs et al. 2006). Optical remote sensing is difficult to use for snow depth retrieval and is strongly influenced by clouds (Dai et al. 2012), 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. 1982). 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. 2006).

Snow hydrological processes are complex, and technology cannot clearly represent changes in these processes. Hydrological modelling, on the other hand, is an efficient way of simulating changes in mass and energy balances in the snow cover (Musselman et al. 2008; Fang & Pomeroy 2009; Musselman et al. 2015). In recent decades, with the development of snow hydrological models from one- to three-dimensional, the ability to simulate snow processes has progressed significantly (Anderson 1973; Wigmosta et al. 1994; Martinec et al. 1998; Bartelt & Lehning 2002; Pomeroy et al. 2007). The majority of these models that simulate seasonal snow evolution do not include blowing snow redistribution processes. The Cold Regions Hydrological Model (CRHM) is a modular model that incorporates appropriate hydrological processes for basins, selected from a library of process modules (Musselman et al. 2015; Weber et al. 2016). Cold regions' hydrological processes are well represented in models such as the Arctic Hydrological and Thermal Process Model (Zhang et al. 2000) and the Variable Infiltration Capacity Model (Bowling et al. 2004). 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. 2014). Rainfall and snowfall are two common precipitation patterns, which are different in the hydrological cycle (Harder & Pomeroy 2014). 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 (Maurer & Mass 2006; Ye et al. 2008; Thériault et al. 2010; Harder & Pomeroy 2013; Kulie et al. 2016), and data utilisation has varied from observation to radar data and methods from empirical to energy balance algorithms. Empirical algorithms are driven by available near-surface meteorological air temperature, which include single and double air temperature threshold methods (Leavesley et al. 1983) and were used in this study to simulate the variation of the snow hydrological cycle. Due to a lack of observations, a method for calculating the critical temperature was developed based on Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover and station interpolated temperature data. Compared with other critical temperature values, these data no longer comprise single values but vary over time and space. Based on these critical temperature data, a distributed snow depth model was built.

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 2009) and snow interception processes (Ellis et al. 2010). 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.

Figure 1

Mass balance of the snow hydrological cycle (left) and structure of the model built in this study (right).

Figure 1

Mass balance of the snow hydrological cycle (left) and structure of the model built in this study (right).

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. 2014). To determine precipitation as snowfall or rainfall, a single air temperature threshold (T0) was applied (Leavesley et al. 1983), where snowfall () is calculated as follows:  
formula
(1)
where is air temperature. At the basin scale, the critical temperature is not a constant parameter but varies with topography, elevation, and other factors (Marks et al. 2013). The parameterisation method based on critical temperature has been used to consider the spatial pattern of critical temperature (Chen et al. 2008), and some researchers have relied on meteorological station data to separate solid or liquid precipitation (Han et al. 2010). In this study, remotely sensed data of daily MODIS snow cover and station interpolated temperature data were used to obtain the critical temperature to separate precipitation types. The biggest problem in identifying snow cover area based on MODIS snow cover data is the cloud cover effect. The annual mean value of cloud-cover percentage in the UYRB reaches approximately 70% for MYD10A1 and 60% for MOD10A1. 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. 2013) 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).
Figure 2

Annual mean cloud-cover percentage of different snow cover data. MOD10A1 and MYD10A1 mean the cloud-cover percentage of MOD10A1 and MYD10A1 data, respectively; MCD means the cloud-cover percentage of data via ‘sensor combination’; SF means the cloud-cover percentage of data through all of the cloud removal steps.

Figure 2

Annual mean cloud-cover percentage of different snow cover data. MOD10A1 and MYD10A1 mean the cloud-cover percentage of MOD10A1 and MYD10A1 data, respectively; MCD means the cloud-cover percentage of data via ‘sensor combination’; SF means the cloud-cover percentage of data through all of the cloud removal steps.

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:  
formula
(2)
where 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 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. was calculated by the following equation (Hedstorm & Pomeroy 1998):  
formula
(3)
where is the maximum intercepted snow load (mm) that can be retained by the forest canopy given the current canopy structure and temperature conditions, and 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, 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 can be calculated as a function of LAI, and the maximum snow load per unit branch area (S) (mm) (Schmidt & Gluns 1991).  
formula
 
formula
The unit used for fresh snow density () is mm. Schmidt & Gluns (1991) made extensive measurements suggesting that the value of the maximum snow load per unit branch area for spruce, which is the main forest type of the UYRB, is 5.9 mm. As is not normally available in the meteorological record, an empirical relationship was used to relate to air temperature. The relationship follows (Hedstorm & Pomeroy 1998):  
formula
where is the air temperature (K) around the canopy.
The amount of snow that unloads from the canopy (Ds) (Niu & Yang 2004; Xia et al. 2009) is a function of canopy temperature and wind speed and was defined as follows:  
formula
(4)
where if . When v = 5 m/s, .

Sublimation module

Snow sublimation is composed of pure snow surface sublimation and forest area sublimation. Pure snow surface sublimation () is often accompanied by blowing snow (). The sublimation in the forest area () 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 () was calculated as follows:  
formula
(5)
where the units of is mm.

Pure snow surface sublimation

For sublimation of snow by heat absorption, the latent heat was used to calculate the snow surface sublimation (:  
formula
(6)
where 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. 2012):  
formula
(7)
where is the blowing snow loss, (mm) is the snow transport in the saltation layer, is the snow transport in the suspension layer and (mm) is the blowing snow sublimation. The rate of mass transport () (mm) by blowing snow across a unit width perpendicular to the wind includes suspension-layer and saltation-layer blowing snow:  
formula
(8)
In place of the numerical integration used to evaluate Equation (8), was evaluated using an approximation of the following form (Essery et al. 1999):  
formula
(9)
where u10 (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 equation (Essery et al. 1999):  
formula
(10)
depends on temperature (Thotpe & Mason 1966) as follows:  
formula
where is the latent heat of sublimation, M is the molecular weight of water, R is the universal gas constant, is the thermal conductivity of air, D is the diffusivity of water vapour in air, and 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 () (mm) is calculated as follows:  
formula
(11)
where (mm) is the snow sublimation on the canopy and (mm) is the sublimation underneath the canopy. The formula of sublimation from canopy interception (Parviainen & Pomeroy 2000) () is as follows:  
formula
(12)
where D is the diffusivity of water vapour in the air, is the air density near the canopy layer, is the Sherwood number, is the specific humidity of air near snow, is the specific humidity of snow surface and C is a scale coefficient (Parviainen & Pomeroy 2000):  
formula
 
formula
 
formula
 
formula
where is the canopy temperature, r is the radius of snow particles, is the windspeed near the canopy, v is the air viscosity factor, 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. 2010):  
formula
(13)
where is the mean daily wind speed and 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.

Figure 3

Energy balance structure of snow layer.

Figure 3

Energy balance structure of snow layer.

The energy balance method strictly follows the energy balance equation to calculate the snowmelt process. The energy budget components are net shortwave radiation (), net longwave radiation , sensible heat , latent heat , soil heat and energy input from precipitation . The net energy () and snowmelt () are calculated as follows:  
formula
 
formula
(14)
where 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 1995):  
formula
(15)
where A is the albedo of the snow surface and is the incident shortwave radiation.

Atmospheric longwave radiation

The function of atmospheric longwave radiation () is given as follows:  
formula
(16)
where is the Stefan–Boltzmann constant, is the air temperature (K), and is the atmospheric emissivity (Campbell 1998).

Snow surface longwave radiation

The snow surface longwave radiation () is computed by the following function:  
formula
(17)
where is the snow emissivity and 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 () can be expressed as follows (Walter et al. 2005):  
formula
(18)
where is the water density, is the specific heat capacity of water and is the heat capacity of snow (these three parameters are all constant), is the rainfall, and is the snowfall.

Soil heat

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

Sensible heat

Sensible heat flux () refers to the energy exchange between the atmosphere and the underlying surface under the temperature change. The methods of calculation follow (Brutsaert 1982):  
formula
(19)
where is the volumetric heat capacity of air and is the aerodynamic drag, which can be calculated as follows (Campbell 1977):  
formula
where is the reference height of wind speed, x is the reference altitude of air temperature, , and are constants and is the windspeed at the height of 2 m. The conversion method of wind speed at any height () is as follows (Allen et al. 1998):  
formula

Latent heat

The latent heat () of the snow surface can be calculated as follows (Zhao et al. 2007):  
formula
(20)
where is the latent heat of sublimation, is the dry air constant, is the air vapour pressure, and is the vapour pressure of snow surface. We assumed that the snow surface vapour pressure is the saturated vapour pressure, , presented as the product of the saturated vapour pressure under air temperature () and relative humidity ().

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 () is negative, there is no snowmelt;

  • (2)

    If is positive and less than the melting energy for net snowfall, only part of the snowfall is melted;

  • (3)
    If 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 was calculated as follows:  
    formula
    (21)
    where is the snow depth on day t–1; , , , and are the snowfall, canopy interception, snow sublimation, and snowmelt on day t, respectively.

STUDY AREA, INPUTS AND PARAMETERS

Study area

The UYRB (Figure 4) with an area of about 40,000 km2 (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.

Figure 4

Location and land use of the study area in the Upper Yangtze River Basin.

Figure 4

Location and land use of the study area in the Upper Yangtze River Basin.

Input data and process

Meteorological data

Meteorological data from 35 stations of the Chinese Meteorological Data Sharing Network (http://data.cma.cn/), including air temperature, precipitation, wind speed, sunshine duration, vapour pressure and relative humidity, were used. To obtain the spatial distribution of meteorological data, the station data were interpolated. The Gradient plus Inverse Distance Squared (GIDS) method combined with multiple linear regressions of climatic gradients and elevation were used to interpolate the station data (Nalder & Wein 1998). The meteorological data were interpolated to a spatially and temporally continuous form. Because all of the stations were not at the same height, we assumed a logarithmic wind profile and vertically interpolated the wind speed data to a common 10 m level using a roughness length of 0.001 m under the assumption of a completely covered snow surface and used the GIDS method to interpolate.

Remotely sensed data

Six types of data were used. (1) MODIS daily snow cover data included MOD10A1 and MYD10A1 datasets with a spatial resolution of 500 m and a temporal resolution of 1 day. (2) Land surface temperature is an important parameter in energy balance computation, and the MODIS level 3 product MOD11A1 LST dataset was used. The study area was covered by four scenes of images, with a spatial resolution of 1,000 m and a temporal resolution of 1 day. (3) NDVI data from the MOD13A3 dataset with a spatial resolution of 1,000 m and temporal resolution of 1 month was used. The effective value range of the data was 2,000–10,000, and the scale conversion coefficient was 0.0001. (4) LAI data from the MOD15A2 dataset with a temporal resolution of 8 days were used. The effective value ranged from 0 to 100, and the scale conversion coefficient was 0.1. All of the MODIS data products used in this study were derived from the Earth Observing System Data and Information System (EOSDIS) (https://earthdata.nasa.gov), which is a key core capability in NASA's Earth Science Data System (ESDS) Program. (5) The passive microwave radiometer is widely used in the field of cryospheric studies (Che et al. 2008). Therefore, the daily SSM/I-based snow depth data with a spatial resolution of 25 km were also used to validate our results and compared with the validation of the station observed values. The SSM/I-based snow depth data that we adopted are widely used in China (http://westdc.westgis.ac.cn/data/df40346a-0202-4ed2-bb07-b65dfcda9368). (6) Snow station observation daily data were obtained from 35 stations distributed in this area in 2003–2010. The pre-processing of the records included the removal of the records of less than 0 cm and missing records. Furthermore, the records with temperatures larger than 0 °C were removed to reduce the influence of wet snow.

Parameters

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

Table 1

Input parameters used in the model

Symbol Definition model dimension Adopted current value Reference 
 Maximum snow load per unit branch area 5.9 mm Schmidt & Gluns (1991)  
 Latent heat of snow sublimation 2.834 × 106 J kg–1 Essery et al. (1999)  
 Molecular weight of water 2.99 × 10−26 kg Essery et al. (1999)  
 Universal gas constant 8.314 J mol–1 K–1 Essery et al. (1999)  
 Thermal conductivity of air 0.0244 W K–1m–1 Essery et al. (1999)  
 Degree of water vapour saturation 1.29 kg m–3 Essery et al. (1999)  
 Radius of snow particular 0.3 mm Essery et al. (1999)  
 Air viscosity factor 17.9 Parviainen & Pomeroy (2000)  
 Density of ice 900 kg m–3 Parviainen & Pomeroy (2000)  
 State factor of snow 0.0114 Parviainen & Pomeroy (2000)  
 Diffused shortwave radiation constant 0.25 Masahiko & Takanori (1995)  
 Received shortwave radiation constant 0.5 Masahiko & Takanori (1995)  
 Stefan–Boltzmann constant 5.67 × 10–8 W m–2 K–4 Constant 
 Density of water 1.0 × 103 kg m–3 Walter et al. (2005)  
 Specific heat capacity of water 4.18 × 103 J kg–1°C–1 Walter et al. (2005)  
 Specific heat capacity of snow 2.09 × 103 J kg–1°C–1 Walter et al. (2005)  
 Volumetric heat capacity of air 1,205 J m–3 k–1 Brutsaert (1982)  
 Reference height of wind speed 2 m Campbell (1977)  
 Reference altitude of air temperature 1.2 m Campbell (1977)  
 Dynamic impedance coefficient 0.001 m Campbell (1977)  
 heat and water vapour impedance coefficient 0.0002 m Campbell (1977)  
1 Von Karman constant 0.4 Campbell (1977)  
 Dry air constant 287 J kg–1 K–1 Constant 
 Snow emissivity 0.97 Constant 
Symbol Definition model dimension Adopted current value Reference 
 Maximum snow load per unit branch area 5.9 mm Schmidt & Gluns (1991)  
 Latent heat of snow sublimation 2.834 × 106 J kg–1 Essery et al. (1999)  
 Molecular weight of water 2.99 × 10−26 kg Essery et al. (1999)  
 Universal gas constant 8.314 J mol–1 K–1 Essery et al. (1999)  
 Thermal conductivity of air 0.0244 W K–1m–1 Essery et al. (1999)  
 Degree of water vapour saturation 1.29 kg m–3 Essery et al. (1999)  
 Radius of snow particular 0.3 mm Essery et al. (1999)  
 Air viscosity factor 17.9 Parviainen & Pomeroy (2000)  
 Density of ice 900 kg m–3 Parviainen & Pomeroy (2000)  
 State factor of snow 0.0114 Parviainen & Pomeroy (2000)  
 Diffused shortwave radiation constant 0.25 Masahiko & Takanori (1995)  
 Received shortwave radiation constant 0.5 Masahiko & Takanori (1995)  
 Stefan–Boltzmann constant 5.67 × 10–8 W m–2 K–4 Constant 
 Density of water 1.0 × 103 kg m–3 Walter et al. (2005)  
 Specific heat capacity of water 4.18 × 103 J kg–1°C–1 Walter et al. (2005)  
 Specific heat capacity of snow 2.09 × 103 J kg–1°C–1 Walter et al. (2005)  
 Volumetric heat capacity of air 1,205 J m–3 k–1 Brutsaert (1982)  
 Reference height of wind speed 2 m Campbell (1977)  
 Reference altitude of air temperature 1.2 m Campbell (1977)  
 Dynamic impedance coefficient 0.001 m Campbell (1977)  
 heat and water vapour impedance coefficient 0.0002 m Campbell (1977)  
1 Von Karman constant 0.4 Campbell (1977)  
 Dry air constant 287 J kg–1 K–1 Constant 
 Snow emissivity 0.97 Constant 

RESULTS AND DISCUSSION

Critical temperature

For the precipitation separation temperature, the spatial distribution value of China was calculated based on the observed precipitation patterns of 643 stations during 1961–1979 (Han et al. 2010). However, the lack of observed data in recent decades has led to difficulties in obtaining the critical temperature. By using the MODIS data, the method for estimating the critical temperature was developed during the study period. The BSD of the study area varied both spatially and temporally (Figure 5).

Figure 5

Spatial distribution of the beginning date of snow duration (from 1 September to 31 August).

Figure 5

Spatial distribution of the beginning date of snow duration (from 1 September to 31 August).

The results show that the BSD in most of the snow area was distributed in the range of 40–60 days after 1 September, while the most widely distributed value was around 60 days, especially in the snow year of 2013–2014. In general, the BSD increased from north to south. There was a significant difference in the northern part of the study area, because of thin snowpack. The results agreed with the real situation based on the station records.

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 BSD and a comparison between observed and simulated BSD values are shown in Figure 7 and Table 2.

Table 2

Comparison of multi-year average values between simulated and observed beginning date of snow duration

Snow station no. Simulated (MM-DD) Observed (MM-DD) Difference (day) 
52908 10-09 10-05 
56004 11-24 11-05 19 
56034 10-26 10-21 
56038 10-13 10-10 
56146 10-18 11-6 –18 
56357 11-18 11-14 
Snow station no. Simulated (MM-DD) Observed (MM-DD) Difference (day) 
52908 10-09 10-05 
56004 11-24 11-05 19 
56034 10-26 10-21 
56038 10-13 10-10 
56146 10-18 11-6 –18 
56357 11-18 11-14 
Figure 6

Variation of multi-year average of beginning date of snow duration (a) and critical temperature (b) in the study area.

Figure 6

Variation of multi-year average of beginning date of snow duration (a) and critical temperature (b) in the study area.

Figure 7

Variation of station-observed beginning date of snow duration (BSD). The lines in the boxes represent the median value of the data series. The BSD starts on 1 September.

Figure 7

Variation of station-observed beginning date of snow duration (BSD). The lines in the boxes represent the median value of the data series. The BSD starts on 1 September.

The results indicate that the changing amplitude of the BSD varied from station to station during 2003–2010. This occurred mainly because of the difference in latitude and elevation. The BSD advanced from high to low latitude and elevation. About half of the recorded BSDs of all stations ranged from the 30th to the 70th day, i.e. October–November, except for stations 56034 and 56146. The differences between simulated and observed BSDs of most stations were small, except for stations 56004 and 56146, with an anomaly at about 20 days. The average difference between simulated and observed BSDs was approximately 9 days in the validation period of 2003–2010. The observed BSD was influenced by terrain, climate factors, and underlying types. Because of the wind and unloading of the canopy, intercepted snow can redistribute the snow mass. The stations in the open field with no trees and wind yielded relatively more accurate results.

Model validation

Remotely sensed snow depth data retrieved from passive microwave brightness temperature data have been used widely (Gao et al. 2012; Wang et al. 2016). Here, observed station data and SSM/I-based snow depth were used to validate the simulated values. According to the pre-processing of the station snow depth records, we obtained daily snow depths from the six stations during 2003–2010, and about 2000 observed snow depth data records were used to validate the model. Figure 8 shows the relationship of simulated and observed snow depth data and simulated and SSM/I-based data. Additionally, the coefficient of determination (R2) and Nash–Sutcliffe efficiency coefficient (NS) were also calculated.

Figure 8

Relationship between observed snow depth and simulated snow depth (a) and between simulated and SSM/I-based snow depth (c) in the Upper Yangtze River Basin from 2003 to 2010. Variation of monthly average snow depth between simulated and observed values (b) and between simulated and SSM/I-based snow depth (d). The dashed blue line is the linear fit trend and the red line is the 1:1 slope. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.079.

Figure 8

Relationship between observed snow depth and simulated snow depth (a) and between simulated and SSM/I-based snow depth (c) in the Upper Yangtze River Basin from 2003 to 2010. Variation of monthly average snow depth between simulated and observed values (b) and between simulated and SSM/I-based snow depth (d). The dashed blue line is the linear fit trend and the red line is the 1:1 slope. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2019.079.

The NS and R2 between simulated and observed snow depth were 0.81 and 0.87, respectively, and these values changed to 0.47 and 0.53 between simulated and SSM/I derived snow depth, respectively. Although both the simulated and SSM/I-based snow depth can explain the seasonal evolution of snow depth in the UYRB to a large extent, there was a comparatively better match between the simulated and observed values. The slope of the linear regression in Figure 8(a) was less than the 1:1 tendency line, which means that the simulated snow depth value was larger than the observed value. The slope of the linear regression in Figure 8(c) was also less than the 1:1 tendency line, indicating that the SSM/I-based snow depth value was larger than simulated value. The difference between simulated and observed snow depth was divided into the same interval to determine the distribution characteristics of the simulated error. The results showed that the difference between simulated and observed snow depth was similar to the normal distribution, with approximately 90% of the difference ranging from −2 to 2 mm (Figure 9). In summary, the distributed snow depth model built in this study indicated the key signals of changes in the real snow depth.

Figure 9

Distribution of the difference between simulated and observed snow depth; a positive value means that the simulated value is larger than the observed value, and a negative value has the opposite meaning.

Figure 9

Distribution of the difference between simulated and observed snow depth; a positive value means that the simulated value is larger than the observed value, and a negative value has the opposite meaning.

Temporal variations in snow components

Based on the critical temperature shown in Figure 6, the snow hydrological processes were simulated for the period of 2003–2015 (Figure 10).

Figure 10

Variation of different snow hydrological processes from 1 August 2003 to 31 July 2015, precipitation type including snowfall, rainfall and the ratio of snowfall and rainfall (a), sublimation type including sublimation underneath the canopy, interception sublimation and snow surface sublimation (b), daily interception (c) and daily snowmelt capacity (d).

Figure 10

Variation of different snow hydrological processes from 1 August 2003 to 31 July 2015, precipitation type including snowfall, rainfall and the ratio of snowfall and rainfall (a), sublimation type including sublimation underneath the canopy, interception sublimation and snow surface sublimation (b), daily interception (c) and daily snowmelt capacity (d).

During the study period, solid precipitation fell mainly in winter, with a ratio of snowfall/precipitation close to 1. Because the daily wind speed was relatively small, the blowing snow losses were negligible. Intercepted and snow surface sublimation comprised approximately 80% of the total sublimation. Furthermore, the interception of the forest canopy and snowmelt capacity both showed significant seasonal change characteristics.

Figure 11 further shows that the variation of daily snowfall within a year had a double-peak pattern. The first peak was on 19 October, with a maximum daily snowfall of 0.45 mm. The second peak was on 10 April with a maximum daily snowfall of 0.8 mm.

Figure 11

Changes of daily average values during the period of 1 August 2003 to 31 July 2015; snowfall (a), canopy interception (b), sublimation (c), snowmelt capacity (d), real snowmelt (e) and snow depth (f).

Figure 11

Changes of daily average values during the period of 1 August 2003 to 31 July 2015; snowfall (a), canopy interception (b), sublimation (c), snowmelt capacity (d), real snowmelt (e) and snow depth (f).

Snowfall intercepted by the canopy is an important part of snow loss, especially in mountainous regions with high tree cover and complex terrain. Commonly, canopy interception has a positive relationship with NDVI and LAI, i.e. it increases with increasing LAI and NDVI. For example, more than 60% of snowfall in the northern forest lands of Canada is intercepted by the canopy (Strasser et al. 2011), and the value in northern Idaho in the USA is 30% (Satterlund & Haupt 1970). Many models in China have not considered this part of the snow loss. In our study, canopy interception was calculated rigorously, and the results indicate that the annual variation of canopy interception was consistent with the snowfall variation, and that less than 10% of snowfall was intercepted by the canopy. The maximum canopy interception was 0.02 mm.

Evaporation, which is called snow sublimation in the snow hydrological processes in cold regions, is an important component of the hydrological cycle. Snow sublimation occurs mainly in regions with open terrain and high wind speeds. Some research has shown that sublimation can occupy more than 80% of the total snow loss in some cold and windy regions, such as the Rocky Mountains in Canada (Macdonald et al. 2010). An important source of sublimation is from canopy interception; in some areas approximately 30% of interception is sublimated (Satterlund & Haupt 1970). Mean wind speed in the UYRB is relatively small, therefore sublimation loss was slight. The amount of snow sublimation in the snow season was 0.03–0.04 mm. The maximum daily value reached 0.06 mm and was less than 20% of all snowfall.

Snowmelt calculated from the energy balance was the main source of snow loss, with more than 70% of total snowfall in this area lost through this process. Snowmelt capacity means the maximum amount of snow that can be melted without considering the limitations of the actual snow amount on the earth's surface. In terms of annual variation, our results show that snowmelt capacity was high in both the early and late periods of the snow year. Early in the snow year (before 20 October), snowmelt 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.

Figure 12

Spatial distribution of the multi-year daily mean snow depth (a) and coefficient of variation (b).

Figure 12

Spatial distribution of the multi-year daily mean snow depth (a) and coefficient of variation (b).

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.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the support of National Key Research and Development Program of China project (Grant No. 2017YFA0603702) and National Basic Research Program of China (973 Program) (2012CB957802). Thanks to all data providers.

REFERENCES

REFERENCES
Allen
R. G.
,
Pereira
L. S.
&
Raes
D.
1998
Crop Evapotranspiration – Guidelines for Computing Crop Water Requirements
.
FAO Irrigation and Drainage Paper 56
.
Food and Agriculture Organization of the United Nations
,
Rome
.
Anderson
E. A.
1973
National Weather Service River Forecast System – Snow Accumulation and Ablation Model
.
NOAA Technical Memorandum NWS HYDRO-17
, p.
217
(copy available on request from: Hydrology Laboratory, Office of Hydrologic Development, NOAA)
.
Bartelt
P.
&
Lehning
M.
2002
A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model
.
Cold Reg. Sci. Technol.
35
(
3
),
123
145
.
Berghuijs
W. R.
,
Woods
R. A.
&
Hrachowitz
M.
2014
A precipitation shift from snow towards rain leads to a decrease in streamflow
.
Nat. Clim. Chang.
4
(
7
),
583
586
.
Bowling
L. C.
,
Pomeroy
J. W.
&
Lettenmaier
D. P.
2004
Parameterization of blowing-snow sublimation in a macroscale hydrology model
.
J. Hydrometeorol.
5
(
5
),
745
762
.
Brutsaert
W.
1982
Evaporation Into the Atmosphere
.
Springer
,
Dordrecht
,
The Netherlands
, p.
299
.
Campbell
G. S.
1977
An Introduction to Environmental Biophysics
.
Springer
,
New York
, p.
159
.
Campbell
G. S.
1998
An Introduction to Environmental Biophysics
,
2nd edn
.
Springer
,
New York
, p.
286
.
Chang
A.
,
Foster
J.
,
Hall
D.
,
Rango
A.
&
Hartline
B.
1982
Snow water equivalent estimation by microwave radiometry
.
Cold Reg. Sci. Technol.
5
(
3
),
259
267
.
Che
T.
,
Li
X.
,
Rui
J.
,
Richard
A.
&
Zhang
T. G.
2008
Snow depth derived from passive mirocwave remote-sensing data in China
.
Ann. Glaciol.
49
(
1
),
145
154
.
Chen
R.-S.
,
Lu
S.-H.
,
Kang
E.-S.
,
Ji
X.-B.
,
Zhang
Z.
,
Yang
Y.
&
Qing
W.
2008
A distributed water-heat coupled model for mountainous watershed of an inland river basin of Northwest China (I) model structure and equations
.
Environ. Geol.
53
(
6
),
1299
1309
.
Chen
R.-S.
,
Liu
J.-F.
&
Song
Y.-X.
2014
Precipitation type estimation and validation in China
.
J. Mt. Sci.
11
(
4
),
917
925
.
Clark
M. P.
,
Slater
A. G.
,
Barrett
A. P.
,
Hay
L. E.
,
Mccabe
G. J.
,
Rajagopalan
B.
&
Leavesley
G. H.
2006
Assimilation of snow covered area information into hydrologic and land-surface models
.
Adv. Water Resour.
29
(
8
),
1209
1221
.
Ellis
C. R.
,
Pomeroy
J. W.
,
Brown
T.
&
Macdonald
J.
2010
Simulation of snow accumulation and melt in needleleaf forest environments
.
Hydrol. Earth Syst. Sci.
14
(
6
),
925
940
.
Essery
R.
,
Li
L.
&
Pomeroy
J.
1999
A distributed model of blowing snow complex terrain
.
Hydrol. Process.
13
(
1415
),
2423
2438
.
Fang
X.
&
Pomeroy
J. W.
2009
Modelling blowing snow redistribution to prairie wetlands
.
Hydrol. Process.
23
(
18
),
2557
2569
.
Gao
J.
,
Williams
M. W.
,
Fu
X. D.
,
Wang
G. Q.
&
Gong
T. L.
2012
Spatiotemporal distribution of snow in eastern Tibet and the response to climate change
.
Remote Sens. Environ.
121
,
1
9
.
Han
C.-T.
,
Chen
R.-S.
,
Liu
J.-F.
,
Yang
Y.
&
Qing
W.-W.
2010
A discussion of the separating solid and liquid precipitations
.
J. Glaciol. Geocryol.
32
(
2
),
249
256
.
Harder
P.
&
Pomeroy
J.
2013
Estimating precipitation phase using a psychrometric energy balance method
.
Hydrol. Process.
27
(
13
),
1901
1914
.
Harder
P.
&
Pomeroy
J. W.
2014
Hydrological model uncertainty due to precipitation-phase partitioning methods
.
Hydrol. Process.
28
(
14
),
4311
4327
.
Hedstorm
N. R.
&
Pomeroy
J. W.
1998
Measurements and modelling of snow interception in the boreal forest
.
Hydrol. Process.
12
(
10–11
),
1611
1625
.
Koivusalo
H.
&
Heikinheimo
M.
1999
Surface energy exchange over a boreal snowpack: comparison of two snow energy balance models
.
Hydrol. Process.
13
(
14–15
),
2395
2408
.
Kulie
M. S.
,
Milani
L.
,
Wood
N. B.
,
Tushaus
S. A.
,
Bennartz
R.
&
L'ecuyer
T. S.
2016
A shallow cumuliform snowfall census using spaceborne radar
.
J. Hydrometeorol.
17
(
4
),
1261
1279
.
Leavesley
G. H.
,
Lichty
R. W.
,
Troutman
B. M.
&
Saindon
L. G.
1983
Precipitation-runoff Modeling System (PRMS) – User's Manual
.
USGS Water-Resources Investigations Report No.83-4238
.
Li
H.
,
Wang
J.
&
Hao
X. H.
2012
Influence of blowing snow on snow mass and energy exchanges in the Qilian Mountainous
.
J. Glaciol. Geocryol.
34
(
5
),
1084
1090
.
Macdonald
M. K.
,
Pomeroy
J. W.
&
Pietroniro
A.
2010
On the importance of sublimation to an alpine snow mass balance in the Canadian Rocky Mountains
.
Hydrol. Earth Syst. Sci.
14
(
7
),
1401
1415
.
Martinec
J.
,
Rango
A.
,
Roberts
R.
,
Baumgartner
M. F.
&
Apfl
G. M.
1998
Snowmelt Runoff Model (SRM) User's Manual
.
University of Berne, Department of Geography
,
Berne
.
Masahiko
H.
&
Takanori
K.
1995
Estimation of snowmelt volume using air temperature and wind speed
.
Environ. Int.
21
(
5
),
497
500
.
Musselman
K. N.
,
Molotch
N. P.
&
Brooks
P. D.
2008
Effects of vegetation on snow accumulation and ablation in a mid-latitude sub-alpine forest
.
Hydrol. Process.
22
(
15
),
2767
2776
.
Musselman
K. N.
,
Pomeroy
J. W.
,
Essery
R. L. H.
&
Leroux
N.
2015
Impact of windflow calculations on simulations of alpine snow accumulation, redistribution and ablation
.
Hydrol. Process.
29
(
18
),
3983
3999
.
Niu
G.-Y.
&
Yang
Z.-L.
2004
Effects of vegetation canopy processes on snow surface energy and mass balances
.
J. Geophys. Res. Atmos.
109
(
D23111
),
1
15
.
Parviainen
J.
&
Pomeroy
J.
2000
Multiple-scale modelling of forest SS: initial findings
.
Hydrol. Process.
14
(
15
),
2669
2681
.
Pomeroy
J.
,
Gray
D.
,
Brown
T.
,
Hedstrom
N.
,
Quinton
W.
,
Granger
R.
&
Carey
S.
2007
The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence
.
Hydrol. Process.
21
(
19
),
2650
2667
.
Revuelto
J.
,
López-Moreno
J. I.
,
Azorin-Molina
C.
&
Vicente-Serrano
S. M.
2015
Canopy influence on SD distribution in a pine stand determined from terrestrial laser data
.
Water Resour. Res.
51
(
5
),
3476
3489
.
Riggs
G. A.
,
Hall
D. K.
&
Salomonson
V. V.
2006
MODIS snow products user guide to collection 5
.
Digital Media
80
(
6
),
1
80
.
Satterlund
D. R.
&
Haupt
H. F.
1970
The disposition of snow caught by conifer crowns
.
Water Resour. Res.
6
(
2
),
649
652
.
Schmidt
R. A.
&
Gluns
D. R.
1991
Snowfall interception on branches of three conifer species
.
Can. J. For. Res.
21
(
8
),
1262
1269
.
Strack
J. E.
,
Liston
G. E.
&
Pielke
R. A.
Sr.
2004
Modeling snow depth for improved simulation of snow–vegetation–atmosphere interactions
.
J. Hydrometeorol.
5
(
5
),
723
734
.
Strasser
U.
,
Warscher
M.
&
Liston
G. E.
2011
Modeling snow–canopy processes on an idealized mountain
.
J. Hydrometeorol.
12
(
4
),
663
677
.
Thériault
J. M.
,
Stewart
R. E.
&
Henson
W.
2010
On the dependence of winter precipitation types on temperature, precipitation rate, and associated features
.
J. Appl. Meteorol. Climatol.
49
(
7
),
1429
1442
.
Thotpe
A. D.
&
Mason
B. A.
1966
The evaporation of ice spheres and ice crustals
.
Br. J. Appl. Phys.
17
,
541
548
.
Walter
T. M.
,
Brooks
E. S.
,
McCool
D. K.
,
King
L. G.
,
Molnau
M.
&
Boll
J.
2005
Process-based snowmelt modeling: does it require more input data than temperature-index modeling?
J. Hydrol.
300
(
1–4
),
65
75
.
Wang
L.
,
Sun
L. T.
,
Shrestha
M.
,
Li
X. P.
,
Liu
W. B.
,
Zhou
J.
,
Yang
K.
,
Lu
H.
&
Chen
D. L.
2016
Improving snow process modeling with satellite-based estimation of near-surface-air-temperature lapse rate
.
J. Geophys. Res. Atmos.
121
,
12005
12030
.
Weber
M.
,
Bernhardt
M.
,
Pomeroy
J. W.
,
Fang
X.
,
Härer
S.
&
Schulz
K.
2016
Description of current and future snow processes in a small basin in the Bavarian Alps
.
Environ. Earth Sci.
75
(
17
),
117
127
.
Wigmosta
M. S.
,
Vail
L. W.
&
Lettenmaier
D. P.
1994
A distributed hydrology–vegetation model for complex terrain
.
Water Resour. Res.
30
(
6
),
1665
1679
.
Xia
K.
,
Li
Y.
&
Liu
W. P.
2009
Snow depths simulated by BATS_SAST model and its improvement
.
J. Glaciol. Geocryol.
31
(
5
),
871
879
.
Ye
H.
,
Yang
D.
&
Robinson
D.
2008
Winter rain on snow and its association with air temperature in northern Eurasia
.
Hydrol. Process.
22
(
15
),
2728
2736
.
Zhang
G.
,
Xie
H.
,
Yao
T.
,
Li
H.
&
Duan
S.
2014
Quantitative water resources assessment of Qinghai Lake basin using Snowmelt Runoff Model (SRM)
.
J. Hydrol.
519
(
Part A
),
976
987
.
Zhao
Q.
,
Liu
Z.
,
Fang
S.
&
Lu
Z.
2007
Improved snowmelt model based on EOS/MODIS remote sensing data
.
Arid Land Geogr.
30
(
6
),
915
920
.