Extension of the Representative Elementary Watershed approach for cold regions: constitutive relationships and an application

The Representative Elementary Watershed (REW) approach proposed by Reggiani et al. (1998, 1999) represents an attempt to develop a scale adaptable modeling framework for the hydrological research community. Tian et al. (2006) extended the original REW theory for cold regions through explicit treatment of energy balance equations to incorporate associated cold regions processes, such as snow and glacier melting/accumulation, and soil freezing/thawing. However, constitutive relationships for the cold regions processes needed to complete these new balance equations have been left unspecified in this derivation. In this paper we propose a set of closure schemes for cold regions processes within the extended framework. An energy balance method is proposed to close the balance equations of melting/accumulation processes as well as the widelyused and conceptual degree-day method, whereas the closure schemes for soil freezing and thawing are based on the maximum unfrozen-water content model. The proposed closure schemes are coupled to the previously derived balance equations and implemented within the Thermodynamic Watershed Hydrological Model (THModel, Tian, 2006) and then applied to the headwaters of the Urumqi River in Western China. The results of the 5-year calibration and 3-year validation analyses show that THModel can indeed simulate runoff processes in this glacier and snow-dominated catchment reasonably well, which shows the prospects of the REW approach and the developed closure schemes for cold regions processes. Correspondence to: F. Tian (tianfq@tsinghua.edu.cn)


Abstract. The Representative Elementary Watershed (REW) approach proposed by
represents an attempt to develop a scale adaptable modeling framework for the hydrological research community. Tian et al. (2006) extended the original REW theory for cold regions through explicit treatment of energy balance equations to incorporate associated cold regions processes, such as snow and glacier melting/accumulation, and soil freezing/thawing. However, constitutive relationships for the cold regions processes needed to complete these new balance equations have been left unspecified in this derivation. In this paper we propose a set of closure schemes for cold regions processes within the extended framework. An energy balance method is proposed to close the balance equations of melting/accumulation processes as well as the widelyused and conceptual degree-day method, whereas the closure schemes for soil freezing and thawing are based on the maximum unfrozen-water content model. The proposed closure schemes are coupled to the previously derived balance equations and implemented within the Thermodynamic Watershed Hydrological Model (THModel, Tian, 2006) and then applied to the headwaters of the Urumqi River in Western China. The results of the 5-year calibration and 3-year validation analyses show that THModel can indeed simulate runoff processes in this glacier and snow-dominated catchment reasonably well, which shows the prospects of the REW approach and the developed closure schemes for cold regions processes.

Introduction
Prediction of hydrological responses at the catchment scale in a changing environment, such as due to natural and human induced climate and landscape changes, is a grand challenge in the hydrological sciences, and is of crucial importance for sustainable water management and hazard mitigation. The current generation of physically based, distributed hydrological models for the upscaling of point-scale balance equations, as set out by Freeze and Harlan (1969), have the drawback that there is no a priori perception about how the micro-scale processes interact with each other and therefore cannot account for the self-organized features that emerge at the macro-scale, i.e. catchment scale, as a result of these interactions (McDonnell et al., 2007). One of the most obvious emergent features frequently recognized in the literature is preferential flow which results from the interactions of subsurface flow, soil and biota, and cannot be easily handled using traditional micro-scale process descriptions. Models such as SHE (Abbott et al., 1986), THIHMS-SW (Ni et al., 2008), among many others, are often criticized for their excessive complexity, parameter identifiability and thus uncertainty in comparison with limited data availability and excessive computational demands.
to these balance equations to make the system of governing equations determinate (Reggiani et al., 1999), the REW approach has the potential to represent the net effects of the spatial and temporal heterogeneity of micro-scale processes (Lee et al., 2007), including the effects of self-organizing features that emerge at the macro-scale without any conceptual jump or inconsistency. For these reasons we argue, as do Reggiani and Schellekens (2003) and Reggiani and Rientjes (2005), that the REW approach can serve as an alternative blueprint compared to the classical Freeze and Harlan (1969) blueprint, providing a key building block to a general modeling framework when combined with the alternative downward approach to modeling (Sivapalan et al., 2003).
However, the set of balance equations arising from the REW approach, on their own, are indeterminate, i.e. they have more unknowns than equations. In other words, the states of the system and the flux exchanges cannot be determined just through the use of the conservation principles which are quite common in continuum mechanics. They have to be supplemented by additional equations, including geometric and constitutive relationships, to make the final set of balance equations determinate. This issue is usually called the closure problem (Zhang et al., 2005) which lies at the heart of the REW approach, as pointed out by Beven (2002Beven ( , 2006. In the original work of Reggiani et al. (1998), each REW is divided into five sub-regions (zones), i.e. an unsaturated zone (u-zone), a saturated zone (s-zone), a concentrated overland flow zone (c-zone), a saturated overland flow zone (o-zone), and the main channel reach (r-zone). Remaining within this system definition, several researchers have proposed different closure schemes and have applied them to hypothetical and real watersheds. Reggiani et al. (1999) developed a set of constitutive relationships through exploitation of the second law of thermodynamics, and subsequently demonstrated the applicability of the REW approach and the corresponding closure schemes in hypothetical (Reggiani et al., 2000(Reggiani et al., , 2001 and real watersheds (Reggiani and Rientjes, 2005). Based on Reggiani et al.'s (1999) closure scheme, Zhang et al. (2005) introduced an interception component and developed a multi-layer soil column definition for the REW formulation and applied the resulting REWASH model to Geer River basin, and improved the representation for subsurface storm flow by incorporating the macropore domain and applied to Alzette River basin (Zhang et al., 2006). Lee et al. (2007) and Tian (2006) surveyed the variety of methods in use for developing the constitutive relationships, and classified them into two primary types, namely, the experimental method and the statistical method. Of these, the former is based on direct regression analysis of laboratory and/or field experimental results, and the latter is based on the upscaling of micro-scale process understanding to the macro-scale, similar to statistical mechanics. With the help of the fully physically based, distributed hydrological model, CATFLOW (Zehe et al., 2001), Lee et al. (2007) proposed a closure scheme that was different from Reggiani et al. (1999Reggiani et al. ( , 2005 by explicitly considering spatial heterogeneity of soil properties. Tian (2006) obtained similar results by coupling the Monte Carlo simulation approach and the 1-D Richards' equation. On the basis of the closure relations developed through up-scaling of the CATFLOW model, Lee et al. (2005) developed the CREW model and applied it to the Weiherbach watershed in Germany and Susannah Brook watershed in Australia.
Despite all this progress, the original REW approach cannot fully account for energy related processes, especially in cold regions which exhibit a complex hydrological regime. In fact, cold regions cover nearly half of the global land area (Yang et al., 2000), at least during the cold season, and influence the mass and energy exchanges between landscape and atmosphere to a large extent, and the water that melts from the glaciers and snowpacks accounts for a significant component of daily water supply in cold regions (McManamon et al., 1993;Williams and Tarboton, 1999). Hydrological processes in cold regions are strongly influenced by energy storage and transfer processes, which, so far, cannot be represented within the Reggiani et al. (1998Reggiani et al. ( , 1999 original formulation of the REW approach. By re-defining the structure and composition of the REW, Tian et al. (2006) have extended the REW approach for cold regions. In their revised formulation, each REW is partitioned into six surface sub-regions and two subsurface sub-regions. Vegetation, snow, soil ice, and glacier ice are added to the existing system that included water, gas, and soil matrix. As a result, energy related processes, i.e. evaporation/transpiration, snowpack and glacier accumulation/depletion, and soil freezing/thawing, can be modeled in a physically reasonable way. The Thermodynamic Hydrological Model (THModel) has been developed by adopting Lee et al.'s (2007) closure scheme within this new extended REW formulation and applied to the semi-arid experimental watershed (Chabagou River basin) in China with good results (Tian, 2006;Tian et al., 2008). However, up till now the closure schemes incorporated in the model completely exclude the energy balance equations and constitutive relationships appropriate for cold regions processes have not been developed within the REW framework. This paper will serve as a companion paper to Tian et al. (2006) and aims to develop parsimonious parameterizations for the snow/glacier melting and accumulation, soil freezing and thawing by taking full advantage of energy balance equations, and to demonstrate the capability of a generalized REW-based cold regions model for long-term streamflow modeling in high mountainous cold regions.
The remainder of this paper is organized as follows. In Sect. 2 we give a brief review of special REW balance equations for cold regions as well as the numerical aspect of the THModel. Following this, in Sect. 3, we will propose appropriate constitutive relationships for the glacier zone, snow covered zone, and soil ice within the unsaturated zone, which will make the coupled balance equations fully determinate. This will be followed up, in Sect. 4, by a case study involving the application of the new THModel to a headwater catchment of the Urumqi River Basin located in northwest of China, and a discussion of the results of model application. This will be followed, in Sect. 5, by a summary of the main results and conclusions.

Brief review of THModel
In the THModel (Tian et al., , 2008, the REW is divided into a surface layer and a subsurface layer (see Fig. 1). Six sub-regions (or zones) are defined in the surface layer, i.e. a bare soil zone (b-zone), a vegetated zone (v-zone), a snow covered zone (n-zone), a glacier covered zone (g-zone), a sub-stream-network (t-zone), and the main channel reach (rzone), while two sub-regions are defined in the sub-surface layer, i.e. an unsaturated zone (u-zone) and a saturated zone (s-zone). The ice phase is introduced into both the u-zone and s-zone to allow soil freezing and thawing to be modeled. The mass, momentum, energy balance equations for each zone are derived in a systematic and extensible way by Tian et al. (2006). All the balance equations are coupled together and solved simultaneously for either temperate or cold regions. The constitutive relationships for temperate regions processes within THModel framework have been worked out and tested at Chabagou watershed (Tian et al., 2008) by following Lee et al. (2007). In this paper we will develop the constitutive relationships for cold regions processes and couple the existing balance equations and constitutive relationships to get a solvable mathematical system. Although the balance equations for all the zones are already there, in the interest of consistency we will here list the special balance equations for cold regions including those for g-zone, n-zone, and u-zone, and make some minor revisions in them for later use (see Table 1). In frozen areas the saturated zone is considered to be beneath the frozen soil layer and mass exchange rate between the s-zone and other sub-regions is rather slow (Lin, 1980). As a first attempt, the balance equations for s-zone are, therefore, ignored in this paper. For u-zone, as shown in Table 1 (Eq. 7), the exchange terms between the u-zone and its environment and among different phases are highly complex. To confine the problem to a manageable level, we assume that (a) the depth and area ratio of u-zone are constant in the frozen areas, (b) water exchange with neighboring REWs or the external world is negligible, (c) water exchange with the g-zone is negligible, and (d) water exchange with s-zone is negligible. Also, we provide here an alternative simpler way to account for the saturation excess runoff by introducing the mass exchange term from the u-zone to the t-zone e ut l instead of changing the area of the t-zone as done in Tian et al. (2008). Given the above assumptions, we obtain the simplified mass balance equation for u-zone after adding Eqs. (7) and (8) in Table 1, which is shown as Eq. (9) in Table 1. For details about the balance equations and the established closure relationships please refer to Tian et al. (2006Tian et al. ( , 2008, and about the symbols in the following equations please refer to the Nomenclature. The numerical part of THModel is implemented with the help of a widely-used ODE solver, CVODE (Cohen and Hindmarsh, 1996), within which the iterative scheme for solving a set of ODEs is based on the Backward Differentiation Formulas (BDE) coupled with Newtonian Iteration for 568 L. Mou et al.: Application of the REW approach for cold regions   nonlinear equations and the preconditioned GMRES algorithm for linear equations.

Constitutive relationships for the energy related processes in cold regions
Special processes in cold regions can be classified into two types. One is the glacier/snow melting and accumulation, which may contribute to annual or seasonal runoff and hence greatly influence the overall water balance, and the other is soil freezing and thawing which may alter soil properties and hence can dramatically influence water and energy transfers in the unsaturated zone. Both of them are driven by relevant energy processes, which are primarily thermal processes. These contribute to the intimate and complex coupling between water and heat transfer and storage processes. Generally speaking, both snow and glacier melting processes are subject to similar driving forces, i.e. the solar radiation and sensible heat flux. We can close the balance equations of the g-zone and n-zone by explicitly considering the energy transfer processes in a rigorous manner, which is called energy balance method in this paper (see Sect. 3.1 below). It is, however, always difficult in practice to apply the energy balance method to model snow melt for long-term, large-scale applications owing to the large heterogeneity of snow distribution in space and time, the complex structure of snowpack and limited data in cold, remote regions. A simpler Table 1. Balance equations for glacier, snow and unsaturated zone in THModel .

No. Zone
Note Balance equation Liquid phase, mass balance d dt (ρ n l ε n l y n ω n )=e nT l +e nu l +e nt l +e n lg +e n ln 5 n-zone Solid phase, mass balance d dt (ρ n n ε n n y n ω n )=e nT n +e n ng +e n nl 6 Heat balance y n ω n c n d dt T n =l ′ lg (e n lg +e n ng )+l il e n nl +R n n ω n +Q nT +Q nu 7 Liquid phase, mass balance Simplified mass balance ρ l y u ω u d dt ε u l =e ub l +e uv l +e un l +e ut l −ρ i y u ω u d dt ε u i 10 Heat balance ω u y u c u d dt T u =Q ub +Q uv +Q un +Q ut +l il ρ i y u ω u d dt ε u i Note: (1) All heat balance equations are for the whole sub-regions; (2) In this paper we use T to denote the temperature other than θ in Tian et al. (2006) (not to be confused with the vector symbol T to denote the momentum term); (3) Heat exchange term between g-zone and u-zone is omitted for the trivial water exchange term and usually minor area fraction of g-zone, and the accompanying heat fluxes with runoff from g-zone and n-zone are small and thus omitted. method using the widely-used degree-day formula is, therefore, then introduced to provide an alternative approach to close the n-zone balance equations (see Sect. 3.2 below). Mass and energy exchange terms among different subregions are illustrated generally in Fig. 2 and the detailed one for the n-zone is shown in Fig. 3.
3.1 Energy balance method for closing glacier and snow melting/accumulation processes Mass and energy terms are intimately coupled during glacier and snow melting/accumulation processes, which can be clearly seen from the mass balance equations for g-zone and n-zone and the corresponding heat balance equations in Table 1. In order to close the balance equations, we must specify all the redundant unknowns in a physically reasonable way. In the interest of simplicity, we will take the g-zone as an example in the following section while the same method can be applied to the n-zone without too many changes. For the balance equations in the g-zone, i.e. Eqs. (1-3) in Table 1, we take the variables y g , T g , e gt l as independent unknowns, and assume that the variation of g-zone area (ω g ) is rather small and thus ω g itself can be considered as a constant over hydrological time scales, and the water exchange term between the g-zone and the u-zone (e  Precipitation may fall as either rainfall (e gT l ) or snowfall (e gT i ) which depends on the vertical distribution of temperature in the atmosphere. We adopt here, however, a simple yet popular temperature threshold method, i.e. if air temperature is above the threshold then precipitation falls as rain, otherwise precipitation falls as snow.
(2) Areal average net radiation intensity (R g n ) The net radiation varies with extraterrestrial radiation, sunshine duration and terrain features in the basin. The areal average net radiation can be estimated by GIS software (Kumar et al., 1997). For simplicity, we calculate using Eq. (1) as follows: where R g n is the local net radiation intensity absorbed at a representative position with mean elevation in the g-zone (MJm −2 day −1 ), c g rn is a coefficient which incorporates the influence of terrain factors including slope, aspect, and hillshade on solar radiation as well as the spatial heterogeneity within the REW. R g n is affected by both surface solar radiation and long-wave radiation released from the ground, which can be calculated as follows.
where R g s is the local surface solar radiation (MJm −2 day −1 ), and R g nl is local net long-wave radiation (MJm −2 day −1 ), α g is the albedo of the glacier.
The surface solar radiation is determined by the extraterrestrial radiation, sunshine duration, and also the geographic location. The formula recommended by Shuttleworth (1993) and Allen et al. (1998) is adopted in this paper as follows: where R a is the daily extraterrestrial radiation (MJm −2 day −1 ), n sun is the actual duration of sunshine (hour), N sun is the maximum possible duration of sunshine or daylight hours (hour), a s is a regression constant indicating the fraction of extraterrestrial radiation reaching the earth on overcast days, i.e. n sun =0, a s +b s is the fraction of extraterrestrial radiation reaching the earth on a clear day, i.e. n sun =N sun , d r is the inverse relative distance Earth-Sun, ω s is the sunset hour angle, ϕ is the latitude, δ is the solar declination.
The albedo α g depends on the characteristics of glacier surface including wetness, grain size, snow age, snow depth, and snow density (Anderson, 1976), although a good relationship between albedo and daily air temperature could be found. For example, Kang et al. (1992) proposed the following empirical formulae to estimate daily averaged albedo on the basis of observed air temperature in the headwaters of the Urumqi River, which will be adopted in our case study (Sect. 4).
where α snow is snow surface albedo, α ice is glacier surface albedo, T a is the daily mean air temperature ( • C). When it snows the glacier is covered with snow and Eq. (5) can be used to calculate α g , or Eq. (6) can be used. Also from Shuttleworth (1993) and Allen et al. (1998), the long-wave radiation R g nl is estimated as follows: where σ is Stefan-Boltzmann constant, i.e. 4.903×10 −9 (MJK −4 m −2 day −1 ), T a,max is the daily maximum absolute temperature (K), T a,min is the daily minimum absolute temperature (K), e a is the actual vapor pressure (kPa), and R g so is clear-sky radiation (MJm −2 day −1 ).
(3) Heat exchange rate with atmosphere (Q gT ) Q gT can be divided into two components. One is the energy flux associated with rainfall or snowfall which is denoted by Q gT 1 , the other is the sensible heat flux caused by turbulent air motions above the glacier surface, which is denoted by Q gT 2 . This is slightly different from Tian et al. (2006), which includes only the latter term.
where c pl , c pi , c pa are, respectively, the specific heats of water, ice and air at constant pressure (kJ/(m 3 ·K)), ρ l and ρ i are, respectively, the density of water and ice (kg/m 3 ), T gs is the glacier surface temperature, which can be replaced by T g , D g h is the turbulent transfer coefficient estimated using Eqs. (10-13) below (Price and Dunne, 1976): where D g hn ,D g hs , D g hu are the turbulent transfer coefficient under neutral, stable, and unstable conditions, respectively, κ is von Karman's constant (≈0.40), u is wind speed (m/s) measured at height z ′ (m), Z is depth of glacier and its covers above ground surface (m), d is the zero plane displacement height (m), z 0 is the aerodynamic roughness height over the glacier surface (m), and R i is called the Richardson number. Glacier surface evaporation is the combined phase transition term either by evaporation from water or sublimation from ice, which will contribute to the latent heat flux across the glacier surface. Usually it is hard to distinguish between evaporation and sublimation and therefore they can be calculated together.
where q s is the saturation value of specific humidity of the glacier surface at T gs , q a is the specific humidity of atmosphere (kg kg −1 ), ρ a is the density of air (kg/m3), l ′ lg is the combined latent heat of vaporization and sublimation (kJ/kg), l lg and l ig are the latent heats of vaporization and sublimation, respectively (kJ/kg), the meanings of other symbols are similar to those of Eqs. (8) and (9).
It should be noted that the above equations for sensible and latent heat flux are point-scale equations and are subject to the effects of spatial heterogeneity due to the spatial variability of temperature, wind speed, and glacier depth. In the interest of simplicity, areal averaged sensible and latent heat fluxes are replaced by the ones expressed by Eqs. (9) and (15), estimated at mean elevation of watershed and multiplied by coefficients c gT 2 and c g lg , respectively, which is similar to Eq. (1). More sophisticated schemes explicitly accounting for the effects of the actual heterogeneity present could be adopted or developed, this is left for future research.
When the total energy supply exceeds the negative heat storage of the g-zone, i.e. l ′ lg (e g lg +e g ig )+R g n ω g +Q gT >y g ω g c pi 0−T g , glacier melting will occur and the temperature of g-zone will be kept at zero, i.e. T g ≡0 • C. We thus obtain the following equation: −l il e g il = max l ′ lg (e g lg +e g ig )+R g n ω g +Q gT +y g ω g c pi T g , 0 where l il is the latent heat of melting.
3.2 Degree-day method for closing snow melting and accumulation equations The energy balance methods discussed in Sect. 3.1 for the g-zone can also be applied to the n-zone without great changes. The two exceptions are (a) the net solar radiation is different for diverse albedo of the snow surface, and (b) the infiltration of snowmelt water into the u-zone cannot be neglected. However, the data required for this method is not easy to collect, especially in remote regions, due to the broad extent and spatial heterogeneity of snowpack. We therefore propose a more parsimonious approach to close the snow balance equations based on the widely-used but conceptual degree-day method (Finsterwalder and Schunk, 1887). Using the degree-day method we can omit the heat balance equation, and consequently the only independent unknown is (y n ) and dependent unknowns are ω n , e nT l , e nT n , e n lg , e n ng , e nu l , e n nl , e nt l , of which e nT l and e nT n are, respectively, rainfall and snowfall separated from precipitation by the threshold temperature method.
(1) Area fraction of snow zone (ω n ) ω n is the area ratio of snow covered zone which can be related to snow water equivalent with the help of the snow cover depletion curve (Luce et al., 1999), which can be expressed as follows: where W a is the snow water equivalent of the n-zone, and W a,max is the maximum snow water equivalent of the n-zone. The snow depletion curve embodies the snow cover change resulting from snowmelt, and further influences the area of the remaining zones. Here we assume the following principles: when the snow-pack expands, the areas of the b-zone, v-zone, and t-zone will decrease in that order; conversely, when snowpack contracts, the areas of the three zones will increase in the reverse order. Moreover, the area of the r-zone will keep constant whenever the n-zone expands and contracts. The specific shape of the curve can be determined by snow course observations or by empirical formulae (Luce et al., 1999).
(2) Snow surface evaporation (e n lg +e n ng ) Snow surface evaporation is the total evaporation rate from snow surface including evaporation and sublimation, which is calculated by the following formula: e n lg + e n ng = c n lg × E 0 × ρ l × ω n where E 0 is the evaporation rate obtained from the evaporation pan, c n lg is the coefficient of snow surface evaporation, which is assumed to be equal to 0.5, as suggested by Lai and Ye (1991).
(3) Infiltration (e nu l ) e nu l is the mass exchange rate between the n-zone and uzone, i.e. infiltration. The infiltration process is greatly influenced by soil properties such as soil moisture and soil temperature. The experimental formula developed by  is adopted in this paper to estimate cumulative infiltration, which can then be transformed into infiltration rate (e nu l ): where INF is the cumulative infiltration volume (mm), c INF is the coefficient accounting for spatial heterogeneity which is set to 1.0 in our case study (Sect. 4), T u K is the initial soil temperature (K), t 0 is the infiltration time (hour), S u 0 is soil surface saturation, S u I is the initial moisture content including water and ice for the frozen soil, which can be expressed as: where ε u is the soil porosity of u-zone. Equation (20) is also applied to estimate the infiltration from the b-zone and v-zone. The difference among different sub-regions lies in the soil surface saturation S u 0 which is assumed here the area ratio of the surface sub-region to u-zone. For example, for infiltration from snow-pack to the u-zone, S u 0 is defined as S nu 0 =ω n ω u , while for the v-zone and b-zone, the soil surface saturation should be S vu 0 =ω v ω u and S bu 0 =ω b ω u , respectively.
(4) Snowmelt (e n nl ) The degree-day equation first proposed by Finsterwalder and Schunk (1887) and widely used is adopted after a minor modification for snowmelt modeling.
where ρ n is the mass density of snowpack, ω n is the area ratio of n-zone, a nl and b nl are the coefficients of snow melting, and T 0 is the critical temperature ( • C).
(5) Runoff generation (e nt l ) The residual term from snowmelt minus infiltration and evaporation will stay in the snowpack until it exceeds the liquid water-holding capacity, estimated from snow density. The density of new snow is proposed to decrease with temperature, while in the case of old snow it is determined by overburden pressure (Anderson, 1976;Jordan, 1991). Besides for short-term calculations, we should also estimate the infiltration rate within the snowpack and the time taken by the melting water to travel through the snowpack to reach the soil surface (Sun et al., 1999).

Soil freezing and thawing
In the balance equations for the u-zone, i.e. Eqs. (9) and (10) in Table 1, we take the variables ε u l , ε u i , T u as independent unknowns, and assume the depth (y u ) and area ratio (ω u ) of the u-zone to be constant. The dependent unknowns to be specified by appropriate closure relationships are e ub l , e uv l , e un l , e ut l , Q ub , Q uv , Q un , Q ut , of which the infiltration terms, e ub l , e uv l , and e un l , have been discussed in Sect. 3.2. However, the total number of unknowns still exceeds the number of available equations so we should bring in one more new equation.
Owing to the matric and osmotic potentials, unfrozen soil water is maintained in a balanced state with ice. Assuming the existence of balance between water potential and vapor pressure over pure ice surface, the maximum unfrozen-water content model is applied to couple the mass and energy balance equations in the u-zone (Hu et al., 1992. where ε u is the soil porosity of the u-zone, T u and T u k are the averaged temperatures of the u-zone in • C and in K, respectively, c is the density of solute in soil, (mol/kg), R is the universal gas constant (8.3143J/mol/K), g is the gravitational acceleration (m/s 2 ), ψ b is the air entry value of the soil matrix potential (m), and B is the coefficient related to soil texture.
According to Eq. (22), the change of unfrozen water content with soil temperature can be expressed as a function of soil temperature, i.e. dε u l dT u =f (T u ). As a result, the mass and heat balance equations of u-zone in Table 1 can be deduced as follows: dt T u =Q ub +Q uv +Q un +Q ut +l il e ub l +e uv l +e un l +e ut l −ρ l y u ω u f (T u ) dT u dt −→ ω u y u c u d dt T u =ω u y u c u × Q ub +Q uv +Q un +Q ut +l il (e ub l +e uv l +e un l +e ut l ) ω u y u [c u +l il ρ l ×f (T u )] −→ ×ρ l y u ω u ρ l y u ω u dε u l dt =ρ l y u ω u f (T u ) Q ub +Q uv +Q un +Q ut +l il (e ub l +e uv l +e un l +e ut l ) ω u y u [c u +l il ρ l ×f (T u )] (24) ρ l y u ω u d dt ε u l =e ub l +e uv l +e un l +e ut l −ρ i y u ω u d dt ε u i Eq. (24) −−−−→ ρ i y u ω u d dt ε u i = e ub l +e uv l +e un l +e ut l −ρ l y u ω u f (T u ) Q ub +Q uv +Q un +Q ut +l il (e ub l +e uv l +e un l +e ut l ) ω u y u [c u +l il ρ l ×f (T u )] (1) Heat exchange rate with above sub-regions (Q uv , Q ub , and Q un ) Q uv , Q ub , and Q un are the heat exchanges between the uzone and the v-zone, b-zone, and n-zone, respectively, which Headwaters of the Urumqi River China The Urumqi River include the heat fluxes associated with infiltration and the heat conduction terms between them. Take the vegetation zone for example, the heat exchange rate is calculated by the following equation: where the first term on r.h.s. represents heat flux associated with mass exchange, the second one accounts for heat conduction, T v is the averaged temperature of v-zone, and D u is the averaged coefficient of thermal conductivity (kW/(m·K)).
(2) Runoff generation (e ut l ) e ut l is the saturation excess runoff which can be calculated as follows.
where ε u c is the soil moisture at field capacity.

Study area: headwaters of the Urumqi River
The development and testing of the constitutive relationships are performed at the headwaters of the Urumqi River (see Fig. 4). The Urumqi River is located within the Tianshan Mountain, in the north-west of China, and originates from Peak Tianger II, located at an elevation of 4486 m, and flows northward to Urumqi City of China. The headwaters of the river are situated in a permafrost region, seasonally covered by snow, and drain an area of 28.9 km 2 with an average elevation of 3860 m. It contains 7 glaciers with a total area of 5.6 km 2 , of which the biggest one is No. 1 Glacier with an area of 1.84 km 2 . The headwaters of the Urumqi River are one of the most heavily instrumented and well studied experimental watersheds in the cold regions of China. There are 3 streamflow gauging stations, i.e. Total Control station (the outlet station), No. 1 Glacier station, and Kongbingdou station, and 1 meteorological station named Daxigou located at the elevation of 3539 m, and many other routine observations of glacier mass and energy balance, glacier movement, and frozen hydrologic issues conducted over many years by different researchers (Yang et al., 2000;Shi et al., 2000).
Runoff in the Urumqi River is contributed from rainfall as well as glacier and snow melting. According to Daxigou meteorological station and the Total Control hydrological station, the mean annual precipitation is about 400 mm, of which over 75% happens during the summer season from May to August. Mean air temperature during the summer season is about 3 • C which leads to melting of glacier and snowpack (see Fig. 5 for the mean monthly precipitation and positive accumulated air temperature). It is reported that runoff volume during the summer season comprised 90% of yearly runoff volume (Yang and Han, 1994). Our case study will focus on the runoff modeling during summer season. Ground surface temperature (K) Difference between ground surface and air temperature (K) Fig. 6. Relationship between ground surface temperature and difference between ground surface and air temperature.

Data pre-processing
(1) Air temperature Air temperature decreases with elevation at rate of −0.68 • /100 m in the headwaters of the Urumqi River according to Ding et al. (1998). We transform the air temperature at Daxigou meteorological station to average air temperature at mean elevation of the study area.
(2) Ground surface temperature The temperature of the t-zone, b-zone, and v-zone are assumed to be equal, and denoted as ground surface temperature. The MODIS images of the watershed in 2002 are used to derive the relationship between daily mean air temperature and temperature difference between air and ground surface, as shown in Fig. 6. The temperature measured by remote sensing image is always at some specific time, e.g. 10:30 a.m. (local time) for the image shown in Fig. 6, which should be transformed to the daily mean temperature with the help of a correspondence curve relating the point value and the daily mean value. (3) Precipitation The precipitation data cannot be used until error corrections are made for use in mountainous regions (Yang et al., 1988). The type of precipitation depends on the air temperature, which can be expressed as follows: According to Kang and Ohmura (1994), the critical temperatures T 1 , T 2 for Daxigou meteorological station are 5.5 • C and 2.8 • C, respectively. The temperature data is, however, the average one at mean elevation of the watershed as mentioned above, and therefore both critical temperatures should be subjected to some calibration.
(4) Snow depletion curve We adopt the snow cover depletion curve obtained by Luce et al. (1999) The initial ice content varies among years and is estimated by averaged daily maximum air temperature from January to April, illustrated in Fig. 7. to some extent (Aziguli et al., 1999). The calibration period is 1990-1994, and the validation period is 1995-1997. The main climatic characteristics during the period of model application are listed in Table 2.
In Sect. 3 we proposed two different closure schemes for g-zone and n-zone. In our case study, the energy balance method (Sect. 3.1) is adopted for the g-zone and the degreeday method (Sect. 3.2) is adopted for the n-zone. We use a single REW for the whole study area (28.9 km 2 ) while allowing for the spatial heterogeneity inherent in the various hydrological processes to be represented by the corresponding parameters in the closure relationships.
The geographic data such as drainage area, slope, channel length etc. are extracted from a DEM with 25 m×25 m resolution (Table 3). The main types of landscape are glacier, snow, swamp, gravel desert, and high mountainous meadow. The soil types underlying the study area include cold desert soils, alpine meadow soil and peat bog soil (Kang et al., 1997).
The parameters subject to calibration are listed in Table 5 which include threshold temperature and coefficients relating to partitioning of precipitation and snow melting degree- day model, and the coefficients accounting for spatial heterogeneity of energy input. We represent the seasonal frozen layer with initial soil ice contents instead of layer depths. The initial soil ice content of each year is determined by the mean daily maximum air temperature from January to April (see Table 4 and Fig. 7). Three standard indices, i.e. Nash-Sutcliffe efficiency coefficient (R 2 ), index of volumetric fit (IVF), and relative error (R E ), are selected to guide mannual calibration, see Eqs. (30)-(32).
where Q oi , Q si are observed and simulated daily discharges, respectively, Q o is the mean observed discharge, n is the total number of time steps. During the 5-year calibration period, we obtained an R 2 value falling within 0.64∼0.80 except for 1992 (0.24), IVF within −0.12∼0.073, and R E within 0.24∼0.41 (see Table 6). The highest R 2 and the lowest I V F and R E are obtained simultaneously in the dry and warm year (i.e. 1990, see Table 2). The abnormal R 2 value for 1992 may be due to the abnormal data in that year, e.g. the discharge exhibits an opposite trend in comparison to precipitation for the two days identified by two dotted lines in Fig. 8, even while having similar air temperatures; the gauged streamflow increased when air temperature decreased in the left circle, while the opposite trend happened in the right circle.
The model predictions during the 5-year calibration period are presented in various panels of Figs. 9 and 10a. Figure 9a Table 5. Calibrated parameters.

Symbol
Value In Eq. Symbol Value In Eq. shows the measured precipitation rates over the simulation period, i.e. May-August in each of the 5 years. Figure 9b presents the corresponding observed air temperature, which is the key to determine the type of precipitation, either snowfall or rainfall, as discussed before. In fact, during the simulation period about 25% of the precipitation is in the form of snow, which is a significant component of the precipitation input. Figure 9c presents the modeled snow water equivalent (SWE). Note that the model simulations are carried out only in a piece-wise fashion, during the May-August period, which means that the initial conditions of the various storage terms needed to be specified. In the case of snow storage, we set the initial storage at a notional amount of 120 mm. In reality, of course, snow accumulation mostly happens during the winter period, although the accumulation might continue into spring also in some years and the maximum snowpack may indeed occur in spring or early summer, as shown in Fig. 9c.
One can also see that as the temperature goes up and remains above the freezing temperature in the May-August period, the snow that is accumulated over the previous months and even after May begins to melt. The difference between the initial snow storage at the beginning of May and the final  storage at the end of August represents (approximately) the net volume of snowmelt in each year. The total volume of snowmelt water in each year will be proportional to the net energy available for the melting (we call it current energy content): to first order, we assume that the energy available can be approximately measured by the average temperature over the May-August time period. This is supported by analysis of total solar radiation measured over the simulation period -it showed a very good reasonable relationship to average air temperature over the same period.
On the other hand, not all of this snowmelt contributes to runoff in the river, since some of the snowmelt can infiltrate into the soil and runoff generation rate from snowmelt is only what is in excess of the local infiltration capacity. In a cold region such as this, infiltration capacity of the soil is dependent not only on soil moisture content but also soil ice content; the latter can tend to make the soil impermeable and reduce the infiltration capacity substantially. Because much of the soil freezing occurs during the cold season, the magnitudes of the soil water and ice contents at the beginning of May depends on the energy that had been available over that period (we call this antecedent energy content): as before, we could assume to first order the relative soil ice and water contents are governed approximately by the average temperature over the January-April period. For example, from Fig. 7 we know that a linear relationship exists between the initial ice content and the average daily maximum air temperature from January to April. The soil ice content can thus be seen as some kind of memory of the antecedent energy content. In other words, the infiltration capacity of the soil is dependent on, and has a memory of the antecedent energy content, as measured by the average temperature over the previous January-April period. In summary, then, the actual amount of snowmelt runoff in each year will be dependent, to first order, on the average air temperature over the May-August period, and the average temperature over the antecedent January-April period. Figure 9d presents the predicted soil ice and water contents over the May-August period. These are predicted by means of an energy balance model of unsaturated zone. The general trends are that as the atmosphere warms the soil ice content decreases while the soil water content increases at the same rate, confirming that the changes represent essentially a simple phase change. However, the initial conditions for both ice and water contents are very different between years, which reflects the extent of freezing that happened in the an-tecedent period. For example, Fig. 9d shows that the first two years (1990)(1991) have relatively higher initial soil ice contents and thus the corresponding soil ice curves occupy a higher position during the May-August period compared to the subsequent three years (1992)(1993)(1994). As indicated by Eq. (19), more ice in the soil will decrease the infiltration capacity of water into the soil column. Therefore, the higher soil ice content will mean lower infiltration capacity and thus higher surface runoff generation ratio, which especially impacts the runoff generation from snow melting greatly due to its moderate melting rate. Figure 9e presents rates of runoff generation predicted by the model, broken up into glacier melting, direct rainfall, and snow melting. Our detailed analysis of these rates shows that, averagely, glacier melting contributes roughly 325 mm to streamflow per year, and is relatively invariant between years. On the other hand, runoff due to rainfall (30 mm) and the melting of snow in the n-zone (30 mm) contribute much smaller volumes per year, although they contribute much more to both the intra-annual and inter-annual variabilities of streamflow, as shown in Fig. 9e. The contribution of snow melting from the glacier zone is less than 10 mm, and relatively constant between years. All of these are assumed to contribute to surface runoff (overland flow). Due to the existence of permafrost during much of the year, lateral subsurface flow is generally weak and is assumed to be negligible. Figure 9e shows, as mentioned above, that glacier melting is the most dominant component except during several intense rainfall events. This is attributable to the lower (almost zero) infiltration capacity of the glacier zone compared to the snow covered zone, the large area ratio associated with glacier in the study catchment. The volume of glacier melt is thus dependent only on the current energy content (i.e. average temperature over the May-August period). Snow melt accounts for a much smaller portion of the total runoff volume as well as the flood peaks, which is attributable to the higher infiltration capacity as well as the lower area ratio associated with the snow pack, when compared to the glacier zone. The daily contribution of snowmelt to the total runoff is seen to be small during the last three years (1992)(1993)(1994) compared to the first two years (1990)(1991), even while exhibiting similiar average air temperatures (see also Table 2). This can be explained by their lower soil ice contents in the last three years. Therefore, one can see that the influence of energy input on the runoff generation is twofold: the current energy input provides the available energy for ice/snow melting which is the source of runoff, while the antecedent energy input determines the infiltration capacity which is a critical threshold to balance the storage and release functions of the watershed. Finally, Fig. 9e also indicates that the rainfall-runoff component is a significant contributor to flood peaks, even if it may not be significant in volumetric terms. This feature implies a dual process control on the volume of runoff and the associated flood peaks in the Urumqi catchment and other similar catchments in cold regions, the understanding of which is crucial to make predictions of the impacts of climate change (i.e. changes in rainfall and/or temperature), with implications for water resources utilization and flood hazard management. Figure 10a shows the correlation between observed discharge and the simulated one via a 1:1 straight line, and we can see that there is no significant systematic deviation for different flow regimes, e.g. high flow regime and low flow regime.

Model validation and sensitivity analysis
Using the parameters obtained by calibration over the period 1990-1994, the model is then validated for the year 1995-1997. The results are shown in Table 6 and Figs. 10b and 11 from which similar conclusions can be drawn from the validation period as from the calibration period (Sect. 4.3). For example, in Fig. 11d, e the years 1996 and 1997 exhibit the opposite extremes for soil ice contents and snow melting induced runoff which more clearly expose the influence of the antecedent energy condition on runoff generation. Note that, as mentioned before, the results from the year 1996 were not used for model validation and also interpretation of the results. This year experienced extreme precipitation in the form of rainfall, and may have changed the landscape, especially the river hydraulic properties, rendering the observations in this and subsequent years subject to considerable uncertainty, and the measures of good fit for the years 1996 and 1997 would not be highly reliable. Note that the R 2 value (0.49) obtained for 1997 is not as good as that for 1995 (0.73).
Earlier, we discussed the first order controls on the annual runoff volume, and argued that the latter may be dependent, broadly, on current energy content (as measured by the average air-temperature over the May-August period), and also on the antecedent energy content (as measured by the average air temperature over the January-April period). Figure 12 articulates this functional relationship in a simple and compact manner. It relates the observed daily mean discharge from 1990 to 1997 (averaged over the May-August period, but excluding the year 1996 for data absence) and the observed daily mean temperatures (representing the energy input) during the antecedent periods (January-April, Fig. 12a)  as well as the simulation periods (May-August, Fig. 12b). The results presented in Fig. 12 confirm (i) the negative correlation between runoff and antecedent energy content (due to its indirect impact on the infiltration capacity of frozen soil and thus losses), in contrast to (ii) the positive correlation between runoff and current energy content (due to the direct impact on melt rates). The outlying point for 1997 in Fig. 12 may be an indication of the possible change of landscape properties that may have been caused by the 1996 flood.
Model sensitivity analysis is simply carried out by one-ata-time perturbation approach on all 9 calibrated parameters in Table 5. The effects of parameter changes on runoff are defined by the relative change in runoff volume (δ V ) and the relative change in Nash-Sutcliffe efficiency (δ R 2 ) shown in Eqs. (33) and (34) (Zhang and Savenije, 2005) with the calibration dataset.
where Q i+ , Q i− , and Q im refer to the discharge at the time step i with the parameter value varied by +50%, −50%, and the manually calibrated parameter value, respectively. The larger δ V , the more sensitive the model is to the parameter under study.
where R 2 + , R 2 − , and R 2 m are the Nash-Sutcliffe efficiency values with the parameter value varied by +50%, −50%, and the manually calibrated parameter value, respectively.
Similiar to δ V , the larger δ R 2 , the more sensitive the model is to the parameter under study.
The sensitivity analysis results (see Table 7) show that runoff simulations are highly sensitive to precipitation partitioning and spatial heterogeneity parameters. The former one, i.e. T 1 , determine how much water in streams comes from melting and how much from direct runoff, and the latter one, i.e. c g rn , influence glacier melting available energy greatly. Our sensitivity analysis results can be explained if we recall the fact that the case study area has a short residence time and that glacier melting is dominant.

Summary and conclusion
As stated by many researchers (Beven, 2002(Beven, , 2006Zehe and Sivapalan, 2007;Reggiani and Schellekens, 2003;Lee et al., 2007), the closure problem is one of the most crucial issues standing in the way of the REW approach becoming an alternative blueprint for distributed hydrological modeling. Although great progress has been made in theoretical and applied aspects of the REW approach in recent years (Reggiani et al., 1998(Reggiani et al., , 1999Reggiani and Rientjes, 2005;Tian, 2006;Tian et al., 2006Tian et al., , 2008Lee et al., 2007;Zhang et al., 2006), energy balance equations have been excluded in these applications. The development and testing of appropriate constitutive relationships for processes occurring in cold regions continue to hamper its application to cold and even temperate regions.
In this paper, within the extended framework of the REW approach provided by Tian et al. (2006), we proposed a set of closure schemes for cold region processes. These were classified into two different types, i.e. glacier and snow melting/accumulation, and soil freezing/thawing. A rigorous energy balance method has been proposed to close the balance equations of melting/accumulation processes, along with the conceptual degree-day method. The closure schemes for soil freezing and thawing are based on the maximum unfrozenwater content model. We applied the proposed closure schemes to the headwaters of the Urumqi River and obtained the promising results. Although in this paper we focus on runoff simulation, the model also produces reasonable results on the variations of internal variables, e.g. snow water equivalent, soil water and ice content, etc. These could be further diagnosed and/or validated by using more observed datasets.
Our results show that the volume of streamflow in the headwaters of Urumqi River basin is totally dominated by the glacier melting due to the large area ratio covered by the glaciers which points to considerable uncertainty about water resource availability under possible climate changes, and the peak value is controlled by both glacier melting and significant rainfall events. However, this conclusion may be restricted within the headwater areas. When the model application is expanded to the mid-and downstream areas, the area ratio covered by the glaciers decreases and the importance of snow melting could increase accordingly. At the same time, the average elevation will be lower and the air temperature should increase and hence the runoff contributed from direct rainfall should also become more important. Also, the model results show the complex influences of energy input on runoff generation, i.e. the antecedent energy input determines the infiltration capacity and the current energy input decides the snow/ice melting amount which points to the importance of seasonal variability of energy input on runoff generation.
It should be noted that the constitutive relationships proposed in this paper for cold regions processes may not occur throughout the year in all cold regions. For example, in seasonally frozen areas soil may behave as frozen soil in the cold season and as normal soil in the warm season, and for seasonally snow-covered areas, snow cover accumulates during the cold season and melts and disappears during the warm season. The simulation of the switching pattern can be accomplished by switching the closure schemes in different seasons, in which case the transition between accumulation and melting should be modeled very carefully. Also, sensitivity analysis suggested that the spatial heterogeneity of energy input, which influences the runoff simulation greatly, should be considered more carefully in the mountainous areas.
In our case study, the runoff contributed by glacier melting turned out to be the most important component, and for this reason in this paper we focused on the melt processes of glacier and deal with the other processes in a relatively simpler manner. Although the closure schemes developed in this paper will not prevent us from focusing on snow melt processes, more detailed analyses should be carried out to investigate snow accumulation and redistribution, and the depletion of snow cover will become dominant in most cold catchments. These difficult issues and the generalization of the REW framework to deal with these are left for future research.

Latin symbols a nl
Coefficient of daily snow melting P Daily precipitation ML −2 T −1 q a , q s Actual humidity and saturated humidity Q gT , Q nT Heat exchange rate with atmosphere of g-zone and n-zone, respectively  Relative change in Nash-Sutcliffe efficiecy Horizontal projected area of watershed L 2 γ r , γ t Slope of r-zone slope of t-zone ω s Sunset hour angle ω g , ω n Time-averaged horizontal projected area ratio of g-zone,n-zone, respectively ω u , ω b Time-averaged horizontal projected area ratio of u-zone,b-zone, respectively ω v , ω t Time-averaged horizontal projected area ratio of v-zone,t-zone, respectively Note: M is the dimension of mass, L is the dimension of length, T is the dimension of time, and is the dimension of temperature. mol is the dimension of mol.