Distributed modeling of landsurface water and energy budgets in the inland Heihe river basin of China

A distributed model for simulating the land surface hydrological processes in the Heihe river basin was developed and validated on the basis of considering the physical mechanism of hydrological cycle and the artificial system of water utilization in the basin. Modeling approach of every component process was introduced from 2 aspects, i.e., water cycle and energy cycle. The hydrological processes include evapotranspiration, infiltration, runoff, groundwater flow, interaction between groundwater and river water, overland flow, river flow and artificial cycle processes of water utilization. A simulation of 21 years from 1982 to 2002 was carried out after obtaining various input data and model parameters. The model was validated for both the simulation of monthly discharge process and that of daily discharge process. Water budgets and spatial and temporal variations of hydrological cycle components as well as energy cycle components in the upper and middle reach Heihe basin (36 728 km2) were studied by using the distributed hydrological model. In addition, the model was further used to predict the water budgets under the future land surface change scenarios in the basin. The modeling results show: (1) in the upper reach watershed, the annual average evapotranspiration and runoff account for 63% and 37% of the annual precipitation, respectively, the snow melting runoff accounts for 19% of the total runoff and 41% of the direct runoff, and the groundwater storage has no obvious change; (2) in the middle reach basin, the annual average evapotranspiration is 52 mm more than the local annual precipitation, and the groundwater storage is of an obvious declining trend because of irrigation water consumption; (3) for the scenario of conservation forest construction in the upper reach basin, although the evapotranspiration from interception may inCorrespondence to: Y. Jia (jiayw@iwhr.com) crease, the soil evaporation may reduce at the same time, therefore the total evapotranspiration may not increase obviously; the measure of changing the farmland to pasture land in the middle reach basin has obvious effects on decreasing evapotranspiration, increasing the discharge at Zhengyixia, and decreasing the storage deficit; reducing the irrigation surface water use in the middle reach basin has obvious functions on increasing the discharge to downstream but the groundwater exploitation increasing should be restricted to prevent the groundwater table decline.


Introduction
The study on hydrological cycle at the river basin scale is of great significance to solve resources and environmental issues, which provides an important scientific basis for the research of water resources, flood forecasting, water environment, as well as water ecology.With the changes of global environment and land surface, the hydrological cycle has been greatly changed.The traditional hydrological statistical methods and lumped hydrological models cannot meet the practical demand (Wang et al., 2003).Physicallybased distributed hydrological simulation could profoundly analyze every component of water cycle process and their inter-relationship in a holistic manner, which can help understanding watershed runoff variation and water resources evolution trend in the changing environment, impacts on hydrological cycle from human activities, and hydrological prediction in ungauged basins (PUBs).Therefore, it is of great significance to research and develop distributed hydrological models (Abbott et al., 1986).
Most land surface models consider the water and energy processes only in the vertical direction, whereas many researchers (Koster et al., 2000;Walko et al., 2000) have noticed that it is very important to reflect the effect of horizontal Published by Copernicus Publications on behalf of the European Geosciences Union.
redistribution of moisture and water on the partition of land surface water and energy fluxes.Therefore, the combination of land surface models and physically-based distributed hydrological models is a research front in the field, just as the VIC (Variable Infiltration Capacity) model does.The WEP (Water and Energy transfer Processes) model (Jia et al., 2001(Jia et al., , 2003) ) applied in the paper is just one of these kinds of models which combine the land surface modeling of vertical water and energy fluxes, and that of the horizontal hydrological processes.Moreover, the modeling of hydrological processes in WEP is physically-based, and it includes not only the overland flow and river flow routings but also the multi-layered groundwater simulation and water utilization to fully reflect the effect of horizontal redistribution of moisture and water on the partition of land surface water and energy fluxes, which makes it different from the VIC model.
Confronting the complexity of land surface hydrology and the existing key issues like heterogeneity/scale, our knowledge of land surface processes is still very poor, and the case study of a specific basin is still meaningful and helpful at present.This paper showing the efforts of a detailed case study in the inland, arid and hilly Heihe river basin is believed to be of referential values.
Since 1960s, especially 1980s, high-intensive human activities have destroyed the natural system in the Heihe river basin, and caused serious ecological and environmental problems, including reduction of flow discharged to downstream, gradual disappearance of some seasonal rivers, degradation of forest and grassland, and frequently occurring of dust storms (Wang et al., 2002).Rational water resources allocation and effective regulation are urgently needed to realize the harmonization between human and natural system.For the practical demand of water resources assessment, allocation and regulation, it is highly desired to develop runoff forecasting models and study the hydrological cycle mechanism and water resource evolution rules in the Heihe river basin.Besides inland mountainous hydro-meteorology, snow melting, soil vegetation and frozen soil, frequent conversion between river water and groundwater, impacts from high-intensive human activities in the middle reach basin and ecological water demand in lower reach regions are also involved in the hydrological system in the Heihe river basin.This is a rather complicated "human-nature" system.
Meaningful explorations on hydrological simulation in the Heihe river basin have been conducted by many scholars and institutions (Cao et al., 2000;CAS-CAREERI, 2000;Kang et al., 2002;Xia et al., 2003;Huang et al., 2004;Zhang et al., 2005;Wu et al., 2005).A conceptual model for monthly runoff, HBV(Hydrologiska Byråns Vattenbalansavdelning) model (Lindstrom et al., 1997) was used to simulate mountainous runoff process by dividing the upper mountainous reaches into two basic landscape zones, including mountainous snow and frozen soil zone, and mountainous vegetation zone (Kang et al., 2002).Distributed Time Variant Gain Model (DTVGM) was developed according to rainfall-runoff nonlinear theory and conceptual simulation method, which was applied in the upper reach basin (Xia et al., 2003).Distributed hydrological model SWAT was applied to simulate the runoff by dividing the upper reach main river into 157 hydrological response units (HRU) (Huang et al., 2004).Simulation of groundwater was conducted in the middle reach basin (Zhang et al., 2005;Wu et al., 2005).Despite these researches have already achieved good results, the existed models are commonly conceptual -based or the adopted spatial and temporal steps are too large.They are difficult to give a detailed description on the physical mechanism of hydrological cycle and accompanied processes, and to effectively analyze the role of every hydrological component on the whole water cycle process, e.g. the spatial and temporal distributions of rainfall, topography, land surface, soil, geology and water consumption, etc.In addition, they are also difficult to predict water resources evolution trend under the impacts of future land surface change and human activity.
In addition, in high-elevation basins like the Heihe river basin, the consideration of frozen soil on land surface hydrology is quite important, especially on river base flow.Experiencing the difficulty in reproducing the base flows in the hydrographs, we suggested an equation to reflect the temperature influence to the saturated hydraulic conductivity in the frozen soil.
The spatial calculation unit of WEP model is a square or rectangular grid.The vertical structure within a grid cell is shown in Fig. 1a, and the horizontal structure within a watershed is shown in Fig. 1b.Each grid cell in the vertical direction, from top to bottom, includes nine layers, namely an interception layer, a depression layer, three upper soil layers, a transition layer, an unconfined aquifer and two confined aquifers.State variables include depression storage on land surface and canopies, soil moisture content, land surface temperature, groundwater tables and water stages in rivers, etc.To consider the subgrid heterogeneity of land use, the mosaic method (Avissar and Pielke, 1989) is used which reflects composition of different land uses within a grid cell.The areal average of water and heat fluxes from all land uses in a grid cell produces the averaged fluxes in the grid cell.Land use is at first divided into three groups, namely a water body group, a soil-vegetation group and an impervious area

Evapotranspiration
Evapotranspiration in a grid cell consists of interception of vegetation canopies, evaporation from water body, soil, urban cover and urban canopy and transpiration from the dry fraction of leaves with the source from the three soil layers.The averaged evapotranspiration E is expressed as: Where Fw, Fsv and Fu are the area fractions of water body, soil-vegetation and impervious area respectively(%), Ew, Esv and Eu are the evaporation or evapotranspiration from them respectively(mm/day). group.The soil-vegetation group is further classified into bare soil, tall vegetation (forest or urban trees) and short vegetation (grass or crops).The impervious area group consists of impervious urban cover and urban canopy.In addition, for the convenience of describing soil evaporation, grass or crop root water uptake and tree root water uptake and reflecting surface soil moisture content change with changing root depth, the surface soil of soil-vegetation group is divided into 3 layers.Runoff routings on slopes and in rivers are carried out by applying one-dimensional kinematical wave approach from upstream to downstream.Numerical simulation of multilayered aquifers is performed for groundwater flows in mountainous and plain areas separately with the consideration of groundwater exchange with surface water, soil moisture and stream flow (Jia et al., 2003).

Evapotranspiration
Evapotranspiration in a grid cell consists of interception of vegetation canopies, evaporation from water body, soil, urban cover and urban canopy, and transpiration from the dry fraction of leaves with the source from the three soil layers.The averaged evapotranspiration E is expressed as: Where F W , F SV and F U are the area fractions of water body, soil-vegetation and impervious area, respectively (%), E W , E SV and E U are the evaporation or evapotranspiration from them, respectively (mm/day).The evaporation from the water body group is calculated with the Penman equation: where RN is the net radiation (MJ/m 2 /day), G is the soil heat flux (MJ/m 2 /day), is the gradient of saturated vapour pressure to temperature (Pa/ • C), δ e is the air vapour pressure deficit (Pa), r a is the aerodynamic resistance (sm −1 ), ρ a is the air density (k gm −3 ), C p is the air specific heat (J kg −1 K −1 ), λ is the latent heat of vaporization (MJ/kg), and γ is the psychometric constant (Pa/ • C).
The evaporation from the impervious area group is taken as the smaller one of current depression storage and the potential evaporation.The maximum depression storage of impervious area is assumed as 2 mm in this study.
The evapotranspiration from the soil-vegetation group is calculated as follows: Where E i is the interception of vegetation (mm/day), E tr is the transpiration from the dry part of vegetation leaves with numbers 1 and 2 representing tall vegetation and short vegetation, respectively (mm/day), and E s is the evaporation from soils (mm/day).
The computation of interception is referred to the Noilhan and Planton (1989) model that is an interception reservoir method, and the transpiration is expressed as: (4) Where Veg is the fraction of tall (or short) vegetation in the soil-vegetation group, δ is the fraction coefficient of the foliage covered by a water film (Noilhan and Planton, 1989), E P M is the Penman-Monteith transpiration (Monteith, 1973), r c is the canopy resistance (s/m), and the others as denoted above.
The aerodynamic resistance and canopy resistance are two important parameters for the evapotranspiration calculation.In the model, the aerodynamic resistance under neutral atmospheric conditions is calculated according to turbulent diffusion theory; the Monin-Obukhov similarity theory is used to modify the computation of aerodynamic resistance under unstable and stable atmospheric conditions; the canopy resistance is calculated by the Dickinson equation (Dickinson et al., 1984).
Evaporation from soils is assumed to come only from the top layer.It is usually estimated by multiplying the potential evaporation (based on the Penman equation) with an evaporation coefficient, which is called the potential method hereafter.However, the potential method may cause theoretical inconsistency of heat flux partition on soil surface because the net radiation and soil heat flux corresponding to the saturated vapour pressure of soil are used in the Penman equation while the actual soil may be unsaturated.Based on the energy balance on the soil surface, aerodynamic diffusion equations of latent and sensible heat fluxes and the wetness function concept (Lee and Pielke, 1992), we derived the following modified Penman equation to compute actual soil evaporation directly: where β is the wetness function; θ is the moisture volume content of surface soil ; θ f c is the moisture holding rate of surface soil; θ m is the soil moisture volume content corresponding to single molecule suction; and the other notations are the same as mentioned above.

Infiltration
Considering infiltration into a vertical uniform soil column when the surface is ponded, Green and Ampt proposed an infiltration model by assuming there is a wetting front which separates saturated soil above from soil below and by using Darcy's law.Compared with the other infiltration models, the Green-Ampt model has the advantages of simplicity, physically-based characteristics and measurable parameters.Mein and Larson (1973) extended it to model infiltration into uniform soil during a steady rain and Moore and Eigel (1981) extended it to model infiltration into two-layered soil profiles during steady rains.Jia and Tamai (1997) suggested a generalized Green-Ampt model for infiltration into multilayered soil profiles during unsteady rains on the basis of Moore and Eigel (1981).The generalized Green-Ampt model is adopted in this study.

Surface runoff
In the water body group, surface runoff is estimated as precipitation minus evaporation.In the impervious area group, surface runoff can be obtained by doing balance analysis of depression storage, precipitation and evaporation on land surfaces.
In the soil-vegetation group, surface runoff consists of two parts, namely the infiltration excess (Hortantype runoff) during heavy rainfall periods and the saturation excess (Dunnetype runoff) during the other periods.A heavy rainfall period Hydrol.Earth Syst.Sci., 13,[1849][1850][1851][1852][1853][1854][1855][1856][1857][1858][1859][1860][1861][1862][1863][1864][1865][1866]2009 www.hydrol-earth-syst-sci.net/13/1849/2009/ is defined as a period during which the rainfall intensity is larger than the saturated soil hydraulic conductivity.The infiltration excess R1 ie is solved by applying the generalized Green-Ampt model to infiltration in three soil layers during heavy rainfall periods.The equations to compute the infiltration excess are as follows: Where H s is the depression storage on soil surface (mm), H s max is the maximum depression storage (mm), P is the rainfall, E 0 is the evaporation, f is the infiltration rate calculated with the generalized Green-Ampt model.In Eq. ( 8), E 0 equals the potential evaporation and transpiration is neglected during heavy rainfall periods.More equations about surface runoff are referred to in Jia et al. (2001).

Subsurface runoff
The subsurface runoff is calculated according to the land slope and the soil hydraulic conductivity: Where R2 is the subsurface runoff from the unsaturated soil layers, L is the channel length in the grid, Slope is the land surface obliquity from a horizontal line and the others as denoted above.

Groundwater flow and groundwater outflow
Taking account of the recharge from unsaturated soil layers and lifted groundwater as source terms, a two-dimensional numerical simulation of multilayered aquifers is performed for groundwater flow to consider the interactions between surface water and groundwater by using the following Boussinesq equations (Zaradny, 1993).Its expression is referred to in Jia et al. (2001).

River flow and overland flow
The river flow and overland flow are routed by using the kinematic wave method: Where A is the area of lateral section, Q is the discharge, q L is the lateral inflow of unit channel length, n is the Manning roughness, R is the hydraulic radius, S 0 is the longitudinal slope of the river bed.

Snow melting process and frozen soil consideration
Although energy balance method provides a good physical basis for describing snow melting process, simple and practical Temperature-index approach is used to simulate daily or monthly snow melting process since too many parameters and data are needed to solve energy balance equation.In this study, the Temperature-index approach is adopted: Where SM is the melt snow and ice (mm/day), M f is the degree-day factor (mm/ • C/day), T a is the temperature index ( • C), T 0 is the threshold temperature for melt ( • C) S is the water equivalent of snowpack SW is the water equivalent of snowfall, E is the sublimation from the snowpack (mm).
The degree-day factor is generally considered as a model calibration parameter since it changes with elevation and seasons as well as land surface conditions changed.The value is usually 1∼7 mm/ • C/day.The value of bare soil is larger than that of grass land, while the latter is larger than that of forest.Daily average temperature is usually selected as the temperature index.The threshold temperature for melting is usually −3∼0 • C. In addition, threshold temperature for distinguishing rainfall and snowfall is needed (usually 0∼3 • C).
The temperature influence to the saturated hydraulic conductivity in the frozen soil is computed as: Where k f is the saturated hydraulic conductivity; k 0 is the saturated hydraulic conductivity when frozen soil melt; T a is the average daily temperature; T c is the threshold temperature for melting; a, b are constants.Through model calibration based on the observed river base flow, it is found that T c =−5 • C, a=0.05, and b=0.25 in the Heihe basin.

Anthropogenic components
The water use in every grid is deduced by using population and water use per capita.The water use per capita is decided according to statistics of water use in a watershed.In addition, water use leakage is deduced from water use and the leakage rate of water supply system.The sewerage is equal to water use subtracted by leakage, and it is set as one part of the lateral inflow to the channel.The groundwater lift is divided into the domestic water, industrial water and the irrigation water.The domestic water is calculated according to the annual regional domestic water lift and the population distribution, and the industrial water is deduced based on the GDP distribution.The irrigation water is calculated based on the annual lift, irrigation area and irrigation rule.
Y. Jia et al.: Distributed modeling of landsurface water and energy budgets

Energy processes
The energy balance equation on land surface is expressed as: Where RN is the net radiation (MJ/m 2 /day), Ae is the anthropogenic energy source (MJ/m 2 /day), lE is the latent heat flux (MJ/m 2 /day), H is the sensible heat flux (MJ/m 2 /day), G is the heat flux into soil (MJ/m 2 /day), RSN is the net shortwave radiation (MJ/m 2 /day), and RLN is the net long-wave radiation (MJ/m 2 /day).

Short-wave radiation
Because of unavailability of the direct observation data of the incoming short-wave radiation, it is estimated from the sunshine hours based on Shimazaki (1996).It at first deduces the direct short-wave radiation and the diffusion short-wave radiation based on the sunshine levels in current and prior hours, and then adds them up to obtain the total short-wave radiation after considering the solar zenith angle.The net short-wave radiation on each land use is obtained after considering the albedo and the sheltering factors.Its expression for each land use is referred to in Jia et al. (2001).

Long-wave radiation
The long-wave radiation is calculated by following Kondo (1994) who separately deduce the downward longwave radiation from atmosphere to landsurface and the upward long-wave radiation from landsurface to atmosphere based on the air temperature and the landsurface temperature.Net long-wave radiation is equal to the downward long-wave radiation subtracted by the upward one.Its expression for each land use is referred to in Jia et al. (2001).

Anthropogenic energy consumption
Statistic energy consumption indices on various types of urban land use are used to consider the impact of human activities on energy balance in an urban area (Kawamata, 1994).
Half of the energy consumption is assumed to emit to land surfaces and the other half to the air.

Latent heat flux
In the hydrological processes part, the computation of evapotranspiration has been described in detail.The latent heat flux λ E can be obtained by multiplying the evapotranspiration with the latent heat of the water λ.

Sensible heat flux
The sensible heat flux H can be expressed as: where ρ a is the density of the air (kg m −3 ), C P is the specific heat of the air (J kg −1 K −1 ), T is the air temperature ( • C), T s is the surface temperature ( • C) and r a is the aerodynamic resistance (sm −1 ).

Heat flux into soil and surface temperature
Based on the energy balance equation we get the heat flux into soil G: The force-restore method (FRM) is used to solve the surface temperature of different land covers.Hu and Islam (1995) suggested an optimal parameter α, which not only ensure minimum distortion of FRM to sinusoidal diurnal forcing but also makes distortion to higher harmonics negligible.They are followed in this research with the equations as follows: where G is the heat flux into soil (MJ/m 2 /day), T s is the surface temperature ( • C), T d is the deep soil temperature (approximated as the daily average of T s , unit: • C), δ is the considered soil depth (selected as the thickness of top soil layer, unit: m), d 0 is the damping depth of the diurnal temperature wave (m), k h is the soil heat conductivity (m/s), c h is the soil volumetric heat capacity (MJm −3 • C −1 ), ω=2π/τ , and τ =86 400.The soil thermal properties depend on the water content and the mineral composition of the soil.The soil heat capacity c h and the soil heat conductivity k h are referred to in Kondo (1994).
3 Model calibration and validation

Input data preparation
The distribution of the rivers and main hydro-meteorological stations in upper and middle reaches of the Heihe river basin is shown in Fig. 2  Based on the DEM of 1:25 000, basic geographic information such as elevation/topography (Fig. 3), slope and overland flow direction are obtained, and the basin is further divided into sub-basins on the basis of stream links deduced from the DEM and observed river information.On the basis of the aforementioned geographic information, flow accumulation is calculated for flow routing calculation.The river section is divided combining with water diversion channels and control stations, and basic information including river length, gradient, top width and bottom width can be obtained.The upper and middle reaches of the Heihe river basin are divided into 6 sub-basins (as shown in Fig. 2) and 70 river sections.Limited to the river channel data, in this study, the river flow routings are only performed for the Heihe mainstream, Babao river and Liyuan river while the river flow routings for other tributaries and mountainous rivers are approximated as overland flow routings in the upper reach watershed.It is noticed that in the flat middle reach basin down Yingluoxia, due to less precipitation, stronger evaporation capacity and soil infiltration capacity, there is little surface runoff generated.Moreover, most river channels disappear in lower reaches due to irrigation water withdraw in the middle reach.Therefore, to simplify the calculation, except keeping the flow routing calculation for the Heihe river mainstream, overland flow routing calculation for the flat middle reaches is omitted in the study, but the cumulative surface runoff is considered as a lateral inflow of the mainstream.
The meteorological data used in the model include daily rainfall (snowfall), temperature, wind speed, sun shine, humidity and maximum evaporation from 1981 to 2002 at 9 meteorological stations within or near the upper and middle reaches of the Heihe river basin: Gaotai, Zhangye, Ye-13/34 reaches due to irrigation water withdraw in the middle reach.Therefore, to simplify the calculation, except keeping the flow routing calculation for the Heihe river mainstream, overland flow routing calculation for the flat middle reaches is omitted in the study, but the cumulative surface runoff is considered as a lateral inflow of the mainstream.

Figure 3 Topography in the upper and middle reach basin
The meteorological data used in the model include daily rainfall (snowfall), temperature, wind speed, sun shine, humidity and maximum evaporation from 1981 to 2002 at 9 meteorological stations within or near the upper and middle reaches of the niugou, Qilian, Shandan, Tuole, Jinta, Dingxin and Jiuquan.The hydrological data used in the model include precipitation, water level and discharge at 11 hydrological stations, among which five stations, Zhamashike, Qilian, Yingluoxia, Gaoya and Zhengyixia, are located at the main river of the Heihe river, two stations, Liyuanbao and Gangoumen, are located at Liyuan river.Monthly discharge data from 1981 to 2002 as well as daily discharge data from 1990 to 2002 are collected.The Thiessen method is used to estimate the meteorological data for each grid, however, the elevation effects on rainfall and air temperature are considered.Based on the correlation analysis of the elevation and the annual average daily values  of meteorological factors, it is found that air temperature and rainfall have close correlations with the elevation (R 2 =0.983 and 0.850, respectively) while sunshine hours, wind speed and relative humidity have low correlation coefficients (R 2 =0.586, 0.002 and 0.381, respectively).Therefore, for every grid cell in the Thiessen polygon controlled by a meteorology gauge station, its air temperature is revised based on a lapse rate of −0.00546 • C/m and the difference of its elevation from the gauge station elevation.The rainfall is also revised in a similar way though it increases with the elevation increase (Eq.24), while the elevation effects on the remaining meteorological factors are neglected.Figure 4 shows the distribution of total precipitation in 2000 in the upper and middle reach basin deduced in the aforementioned approach.taking grid as basic units and the area percentage of each category in each grid is calculated.Figure 5  For 37 irrigation areas having detailed annual water use data , the spatial interpolation of water use (including woodland and grassland irrigation water use) is conducted according to irrigation area distribution and the water use data and is further distributed to the time scale of 10 days according to precipitation, cropping pattern and irrigation rules.For the other irrigation areas, the irrigation quota is set as 7500 m 3 /ha.The domestic water use in every grid is deduced from the distribution of population and daily water use capita.Urban domestic water use quota for current level is 140 L/D, while rural is 50 L/D. Te population in every grid is deduced from the county population in the GIS platform.

Parameter estimation
The main parameters in the model include soil parameters, groundwater aquifer hydraulic conductivity and specific yield, vegetation parameters, roughness of overland and river channel, and infiltration coefficient of river bed.Most parameters in the model are not needed to be calibrated, but 15/34 pattern and irrigation rules.For the other irrigation areas, the irrigation quota is set as 7500m 3 /ha.The domestic water use in every grid is deduced from the distribution of population and daily water use capita.Urban domestic water use quota for current level is 140L / D, while rural is 50L / D. The population in every grid is deduced from the county population in the GIS platform.for several important parameters including saturated soil hydraulic conductivity, groundwater aquifer hydraulic conductivity and specific yield, infiltration coefficient and thickness of river bed as well as the Manning roughness, some adjustments are conducted comparing simulated discharge with observed values during selected calibration period.

Soil
The relation of soil moisture content θ and suction S (negative pressure head) is expressed as (Haverkamp et al., 1977): where θ s=saturated moisture content; θ r=residual moisture content; α, β=constants.
The relation of unsaturated hydraulic conductivity k (θ) and soil moisture content θ is expressed as (Mualem, 1978): where k s =saturated hydraulic conductivity.
In this study, soil moisture characteristics parameters are calculated (Table 1) based on following soil moisture movement experimental research results in the Heihe river basin and adjacent regions: soil moisture movement in the irrigated areas has been studied in the middle reaches of the Heihe river basin (Cao et al., 2000); the soil moisture parameters have been studied in the diversifolia area in the lower reaches of the Heihe river basin (Zhu et al., 2002); the soil moisture movement parameters have been studied in the Qinwangchuan irrigated area (Wang et al., 2002).
The soil heat capacity c h (106 Jm-3K-1) (Kondo, 1994) and the soil heat conductivity k h (Wm −1 K −1 ) (Chung and Horton, 1987) can be expressed as: k h = 0.243 + 0.393 θ + 1.534 θ 0.5 (28) In the above two equations, c hm is the heat capacity of mineral composites, c hw the water heat capacity, θ the soil moisture content and θ s the saturated moisture content.

Groundwater aquifer
Groundwater aquifer hydraulic conductivity and specific yield in the Zhangye Basin are deduced from groundwater simulation and geological exploration data.The saturated soil layer parameters in mountainous areas are deduced from runoff simulation process, hydraulic conductivity is set as 65.5 m/month, and specific yield is set as 0.05.

Vegetation
Vegetation parameters change with seasons.Main parameters of four vegetation types including vegetation coverage rate (veg), leaf area index (LAI), vegetation height (h c ), root depth (lr) and minimum stomatal resistance (r smin ) are shown in Table 3. Aerodynamic parameters for vegetation and other surface covers are shown in Table 4 where z om , z ov , z oh are the roughness heights of momentum, vapor and heat, respectively, d is the display height, and h c is the canopy height.These vegetation and aerodynamic parameters are used to compute the canopy resistance and the aerodynamic resistance in Eq. (5).

Flow routing parameters
Manning roughness used in river flow routing calculation is deduced according to the National Hydrological Yearbook of China, while the Manning roughness used in overland flow

Snow melting parameters
The degree-day factor and the threshold temperature for melt are also model calibration parameters.After model adjustment, the degree-day factor is set as 1 mm/ • C/day for woodland, 2 mm/ • C/day for grassland, 2 mm/ • C/day for bare land, 5 mm/ • C/day for urban land use and 1 mm/ • C/day for glacier, respectively; the threshold temperature for melting is set as 0 • C; the threshold temperature for distinguishing rainfall and snowfall is set as 1 • C. Model calibration and validation For model calibration, the calibration period is 1996 to 2000, the upper and middle reaches are divided into 36 728 units by the square grid cell of 1 km×1 km, and a simulation from 1981 to 2002 is performed with a time step of 1 day.The model is calibrated using daily and monthly discharge processes observed at the Yingluoxia hydrological station, and the main calibration parameters include saturated soil moisture hydraulic conductivity, groundwater aquifer hydraulic conductivity and specific yield, infiltration coefficient of riverbed material as well as the Manning roughness.There are three calibration rules: (1) minimizing the average annual discharge error during the simulation period; (2) maximizing the Nash-Sutcliffe efficiency coefficient; (3) maximizing the correlation coefficient between simulated discharge and observed values.The parameters were adjusted using the "trial and error" method according to the above three rules: firstly, determine the initial value according to physical properties, experimental data and reference data; then adjust the parameters.Specifically, saturated soil moisture hydraulic conductivity, groundwater aquifer hydraulic conductivity and specific yield as well as riverbed infiltration coefficient were adjusted according to www.hydrol-earth-syst-sci.net/13/1849/2009/ Hydrol.Earth Syst.Sci., 13, 1849-1866, 2009    : 1982: to 1995: , and Period 2: 2001: to 2002.In other words, there are two validation periods due to the lack of observed monthly data between them.Similarly, due to the lack of observed daily data during 1982 to 1989, the two validation periods for daily discharges process were, Period 1: 1990 to 1995, which is different from that of monthly discharge process, and Period 2: 2001 to 2002.

Monthly discharge process
The monthly discharge process at the Yingluoxia hydrological station in the upper reaches of the Heihe river basin is well simulated with a high accuracy (Fig. 6).The simulation results are shown in Table 5.It can be found that the model is well validated in the two validation periods based on the relative error of annual average discharge, the Nash -Sutcliffe efficiency coefficient and the correlation coefficient.
Although the monthly discharge processes at other hydrological stations such as Zhamashike are not used to calibrate the model, the comparison of simulated results with observed values is conducted, which shows that the simulation results are close to the observed values (Fig. 7).

Daily average discharge process
The comparison of simulated daily discharge processes with observed values at the Yingluoxia hydrological station is shown in Fig. 8, the simulation results are shown in Table 6.From Fig. 8 and Table 6 we can find that the model is well validated in the two validation periods.Meanwhile, it is also found that the simulation accuracy of daily discharge process is lower than that of monthly discharge process.This is because: (1) the precipitation data at the daily scale are more spatial-temporally variable than those at the monthly scale, so the shorter time scale, the more obvious deviation between the spatial precipitation data obtained by spatial interpolation Hydrol.Earth Syst.Sci., 13,[1849][1850][1851][1852][1853][1854][1855][1856][1857][1858][1859][1860][1861][1862][1863][1864][1865][1866]2009 www.hydrol-earth-syst-sci.net/13/1849/2009/ 20/34 Although the monthly discharge processes at other hydrological stations such as Zhamashike are not used to calibrate the model, the comparison of simulated results with observed values is conducted, which shows that the simulation results are close to the observed values (Fig. 7).

Daily average discharge process
The comparison of simulated daily discharge processes with observed values at the method and the actual observed values.This deviation will be possibly passed to the simulation of discharge process, thus the simulation accuracy at daily scale will be lower than that at monthly scale; (2) the actual dynamic mechanism of basin hydrological cycle is very complex, the mathematics description of the mechanism and method dissolving the spatial variety of land surface and soil properties value in the model could be further improved.However, it is necessary to perform long time series of hydrological simulation at daily or even shorter time scale, the reasons are not only for the need of predicting daily discharge process, but also for the need of predicting future land surface and climate change impacts, which is determined by the dynamic mechanism of hydrological cycle and accompanied processes.(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) 14.61 14.10 −3.5% 0.78 0.89 4 Model application

Water and energy budgets and spatial and temporal variations
A continuous simulation from 1981 to 2002 is performed with a grid size of 1 km and a time step of 1 day using the distributed hydrological model WEP-HeiHe to study detailed processes of water budgets and hydrological cycle components.As an example, the water budget of 1999 in the upper reach watershed (with an area of 9999 km 2 ) is shown in Fig. 9. From the figure we can see that the average evaporation in 1999 in the region is 258 mm, of which 123 mm comes from land surface and vegetation interception, and 135 mm comes from soil evaporation and vegetation transpiration; among the total 156 mm runoff, surface runoff, subsurface runoff and groundwater flow accounts for 70 mm, 24 mm and 62 mm, respectively, and about 2 mm of the total runoff is taken by farmland irrigation.According to the observed data at the Yingluoxia hydrological station, the runoff amount in 1999 is 1610 million m 3 , which amounts to 161 mm runoff depth, while the simulated annual runoff depth is 154 mm, the relative error is 4.4%.The daily simulation results show that, in the upper reach watershed, the annual average evapotranspiration and runoff during the period from 1982 to 2002 accounts for 63% and 37% of the annual precipitation (428 mm), respectively.47% of mountainous runoff directly comes from surface runoff, and the other 53% comes from mountainous subsurface runoff and groundwater flow.The snow melting runoff accounts for 19% of the total runoff and 41% of the direct runoff.It should be noted that the snow melting runoff here include all the surface runoff when rain falls as snow, not only glacier melting runoff.
From the view of annual variation, the evaporation and runoff in upper reaches increases as precipitation increases, but the increase rate of evaporation is less than that of runoff; storage change reflects the objective laws of water storage in flood years and water consumption in drought years (Fig. 10).Table 7 gives the annual simulation results of water budgets in the lower plain area from Yingluoxia to Zhengyixia, i.e. the middle reach basin with an area of 19014 km 2 (not including the Liyuan river watershed) from 1982 to 2002.It can be seen that the annual average precipitation is 240mm, 56% of that in upper reach basin; the annual average runoff is 76 mm, 47% of that in upper reach basin; and the annual average evapotranspiration is 292mm, 108% of that in upper reach basin.For the annual average runoff, 30% is from surface runoff and 70% is from groundwater and Table 7 gives the annual simulation results of water budgets in the lower plain area from Yingluoxia to Zhengyixia, i.e. the middle reach basin with an area of 19 014 km 2 (not including the Liyuan river watershed) from 1982 to 2002.It can be seen that the annual average precipitation is 240 mm, 56% of that in upper reach basin; the annual average runoff is 76 mm, 47% of that in upper reach basin; and the annual average evapotranspiration is 292 mm, 108% of that in upper reach basin.For the annual average runoff, 30% is from surface runoff and 70% is from groundwater and spring.The annual average irrigation water is 145 mm (i.e.2.75 billion m 3 , and 10% is from groundwater), mainly utilized in the middle reach Zhangye Basin and Shandan-Minle Plain.Although the great amount of irrigation water use has a contribution in the local runoff generation, it leads to the annual average evapotranspiration 52 mm more than the local annual precipitation (240 mm).The deficit generally equals to the difference between the upper reach inflows and the outflow at Zhengyixia.
www.hydrol-earth-syst-sci.net/13/1849/2009/ Hydrol.Earth Syst.Sci., 13, 1849-1866, 2009  It can be seen that the annual average net radiation in the middle reach basin is slightly larger than that in the upper reach basin because of the topography influence, the sensible heat flux is similar, and the latent heat flux is also a little larger because of a large amount water use in the agricultural land.
The distribution of actual evapotranspiration, surface runoff, groundwater outflow and sensible heat flux in the upper and middle reach basin in 2000 is shown in Fig. 11.From the figure we can see that, surface runoff mainly generates in upper mountainous reaches, runoff in valleys and low-lying areas mostly come from the groundwater outflow from soil layers in mountainous areas.There is nearly no runoff generated in the lower reach basin except Linze low-lying areas and bare rock regions.Groundwater is recharged by river water seepage from Yingluoxia to Zhangye, while groundwater overflows to river from Zhangye to Zhengyixia.In addition, the distribution of evapotranspiration shows that the Zhangye Basin, the Minle-Shandan irrigation areas as well as artificial oases are the main evapotranspiration zones.Limited to data obtained, the distributions of soil moisture content, groundwater level and snow simulated by the model have not been validated yet.In the future, ground and remote sensing observed data will be used to validate the simulation results in order to support basin water resources management.

Water budgets prediction under the future land surface change
The distributed hydrological model developed in the study is based on the hydrological mechanism of the Northwest inland basin of China.The model parameters can be deduced from the physical properties of land surface and do not change with the changing precipitation.Therefore, the model may be used to predict the water budgets and water resources under the future land surface change, and to conduct scenario analysis under continuous deteriorated conditions and estimate the effects of adopting countermeasures.

Impacts of conservation forest construction in the upper reach basin
Assuming that the land covered with sparse forest in the upper reach basin with an area of 249.6 km 2 , which accounts for 2.5% of the entire area of the upper reach basin, is changed to forest, using the conditions of precipitation and Hydrol.Earth Syst.Sci., 13, 1849-1866, 2009 www.hydrol-earth-syst-sci.net/13/1849/2009/  Table 7 gives the annual simulation results of water budgets in the lower plain area from Yingluoxia to Zhengyixia, i.e. the middle reach basin with an area of 19014 km 2 (not including the Liyuan river watershed) from 1982 to 2002.It can be seen that the annual average precipitation is 240mm, 56% of that in upper reach basin; the annual average runoff is 76 mm, 47% of that in upper reach basin; and the annual average evapotranspiration is 292mm, 108% of that in upper reach basin.For the annual average runoff, 30% is from surface runoff and 70% is from groundwater and meteorology in 1999, the effects on water budgets in the upper reach basin simulated by the model are shown in Table 9.The simulation results show that, due to the construction of conservation forest in the upper reach basin, although the evapotranspiration from interception increases, the soil evaporation reduces at the same time, therefore the total evapotranspiration does not increase obviously; the runoff outflow slightly decreases, while the storage change of soil moisture and groundwater slightly increase mitigating the deficit of storage change in drought years such as 1999.It should be noted that the more important effects of conservation forest construction are to preserve the soil and water, decrease soil erosion in the upper reach and avoid sand sedimentation in the lower reach, thus improve the ecological environment in the basin.

Impacts of changing farmland to pasture land in the middle reach basin
Assuming that 20% of the farmland with an area of 0.37 million hectares in the middle reach in current level is changed to pasture land, using the conditions of precipitation and meteorology in 1999, the effects on water budgets in the middle reach basin simulated by the model are shown in Table 10.
It can be seen that, due to the measure of changing the farmland to pasture land in the middle reach basin, irrigation water consumption is greatly reduced, evapotranspiration and the storage change deficit decreases obviously, while the discharge at Zhengyixia increases obviously.Therefore, the measure is urgently needed in the inland basin with serious water shortage problems.

Impacts of reducing irrigation surface water consumption in the middle reach basin
To guarantee the water demand of downstream oases and lakes, the Heihe river administrative bureau began to control the water diversion from the Heihe river in the middle reach basin in recent years.Compared with that in 1990s, the groundwater table in the middle reach basin in 2004-2005 declined because of the reduced water diversion from the Heihe river and the increased groundwater exploitation, and the decreased amount of groundwater directly transformed from the diverted surface water was estimated to account for 45% of the total groundwater decrease, thus readjusting the industrial structure and saving the agricultural water consumption are significant measures to hold all environments stably (Wei et al., 2008).Limited to data obtained, the distributions of soil moisture content, groundwater level and snow simulated by the model have not been validated yet.In the future, ground and remote sensing observed data will be used to validate the simulation results in order to support basin water resources management.The distributed hydrological model developed in the study is based on the hydrological mechanism of the Northwest inland basin of China.The model parameters can be deduced from the physical properties of land surface and do not change with the changing precipitation.Therefore, the model may be used to predict the water budgets and water resources under the future land surface change, and to conduct scenario analysis under continuous deteriorated conditions and estimate the effects of adopting countermeasures.

Impacts of conservation forest construction in the upper reach basin
Assuming that the land covered with sparse forest in the upper reach basin with an area of 249.6 km 2 , which accounts for 2.5% of the entire area of the upper reach basin, is changed to forest, using the conditions of precipitation and meteorology in 1999, the effects on water budgets in the upper reach basin simulated by the model are shown in Table 9.The simulation results show that, due to the construction of conservation forest in the upper reach basin, although the evapotranspiration from interception increases, the soil evaporation reduces at the same time, therefore the In this study, we set a scenario of reducing the irrigation surface water use by 50% in 2002 while keeping other conditions including the groundwater exploitation unchanged to study the hydrological impacts.The average irrigation quota in the middle reach basin is as high as 8025 m 3 /ha and about 90% is from surface water diversion according to statistics, a reduction of 50% is assumed as an ideal scenario, which is very strict but still possible considering the case in some irrigation districts of high efficient water use in Zhangye city in the basin.The comparison of monthly discharge at Zhengyixia under history and scenario conditions in 2002 is shown in Table 11.The results show that reducing the irrigation surface water use in the middle reach basin has obvious functions on increasing the discharge at the Zhengyixia.However, it should be mentioned that strict water-saving measures should be pursued so as not to cause the increase of groundwater exploitation, otherwise the groundwater table decline issue will be sharpened as pointed out by Wei et al. (2008).

Conclusions
In this study, a physically and GIS based distributed hydrological model was developed in the Heihe river basin.A simulation of 22 years was carried out with a grid size of 1 km and a time step of 1 day in the upper and middle reach basin, and the model was validated using observed runoff data.The validation results show that the model has high simulation accuracy.The characteristics of the model are listed as following: (1) different from the existing methods of runoff simulating and flow routing, the model simulates all the processes of hydrological cycle, and couples the simulation of hydrological cycle process with energy cycle process; (2) to consider the subgrid heterogeneity of land use within a calculation unit, the mosaic method is used by dividing the land use into 8 categories, and water and heat fluxes for each land use is calculated separately, thus greatly improve the calculation accuracy using large grids in large basin, while only one dominant type of land use is considered in most of the existing models; (3) the model developed in the study has a quick calculation speed (the calculation time is about 15 min for 1-year simulation of 36 728 grid units in the upper and middle reach of the Heihe river basin using a PC computer with CPU2.2GHz).In addition, the parameters in the model are all physically-based and can be measured or deduced, and the model was well validated.Therefore, the model has the potential to be applied to PUBs.
The results of applying the distributed hydrological model to the upper and middle reach of the Heihe river basin show that, the model can be used to effectively reveal the temporal and spatial variations of hydrological cycle and water budgets components as well as energy cycle components.The results of applying the model to predict the water budgets under the future land surface change show that, due to the Hydrol.Earth Syst.Sci., 13,[1849][1850][1851][1852][1853][1854][1855][1856][1857][1858][1859][1860][1861][1862][1863][1864][1865][1866]2009 www.hydrol-earth-syst-sci.net/13/1849/2009/ construction of conservation forest in the upper reach basin, although the evapotranspiration from interception increases, the soil evaporation reduces at the same time, therefore the total evapotranspiration does not increase obviously; the runoff outflow slightly decreases, while the storage change of soil moisture and groundwater slightly increase at the same time.The measure of changing the farmland to pasture land in the middle reach basin has obvious effects on decreasing evapotranspiration, increasing the discharge at Zhengyixia as well as decreasing the storage change deficit.Reducing the irrigation surface water use in the middle reach basin has obvious functions on increasing the discharge to downstream but the groundwater exploitation increasing should be restricted to prevent the groundwater table decline.However, an important issue one must keep in mind is that our efforts of the hydrologic-response simulations would be impacted by the issue of equifinality, which refers to more than one parameter combination providing an equally good (or poor) representation of the integrated hydrologicresponse (Ebel, 2006).In this study, the WEP model is based on water balance, energy balance and hydro-physical process and parameters are physically-based.The domains of parameter values and their combinations would be constrained with the physics.Without the physics, there would no such constraints, and any combinations can lead to similar effects (Savenije, 2001;Loague and VanderKwaak, 2004).With the help of spatial technologies like GIS, the use of distributed land surface information for the estimation of some WEP parameters would partially keep the parameter set within reasonable bounds in our rainfall-runoff modeling efforts.However, it can't still be concluded that such a modeling effort has already eliminated the effects of equifinality.Therefore, as commented by many researchers, it is necessary to further discuss and research the equifinality issue.
In order to reduce the impact of equifinality, it would be better to further validate the distributed variables (e.g.ET, soil moisture, groundwater table, etc.).However, it is difficult for most models to validate in all their details.As stated by Beven (2005): Many hydrological variables are of obvious spatial heterogeneity.ET, soil moisture, groundwater table and energy fluxes, for example, might be predicted as an average over a model grid element and over a certain time step; the same variable might be measured at a point and a limit scale in space.Equifinality also relates to errors of input data and boundary conditions.Even so, new data sources for observation of hydrological processes like remote sensing data may alleviate some of the problems facing the validation and operational use of hydrological models, because they can provide systematic land surface observations (e.g., rainfall, snow, soil moisture, vegetation, surface temperature, energy fluxes, and land cover large areas, etc.) consistently over the large scale (Pan, 2007).Therefore, the equifinality of the results in this study will be expected to study with the aid of remote sensing data.

Figure 1
Figure 1 Structure of WEP model: (a) vertical structure within a grid cell and (b) horizontal structure

Figure 2 Fig. 2 .
Figure 2 Main rivers and hydro-meteorological stations in the upper and middle reach basin Fig. 2. Main rivers and hydro-meteorological stations in the upper and middle reach basin.

Fig. 3 .
Fig. 3. Topography in the upper and middle reach basin.

Figure 4
Figure 4 Distribution of total precipitation in 2000 in the upper and middle reach basin Fig. 4. Distribution of total precipitation in 2000 in the upper and middle reach basin.
shows the land use distribution in 2000 as an example.The area percentage of each category in each grid in other years from 1981 to 2002 is obtained by linear interpolation method based on that of 1985, 1995 and 2000.The soil type is generalized into 3 types, sand, loam and clay.The aquifer thickness of Zhangye Basin in the middle reach is deduced from the geological survey data.The thickness of soil layer in the upper mountainous reaches of Yingluoxia as well as Liyuanhe Baisn is calculated in accordance with the law of diminishing exponentially with the distance from the valley, and the thickness of soil layer in the other regions is set as 2 m.

Figure 5
Figure 5 Land use (2000) distribution in the upper and middle reach basin Fig. 5. Land use (2000) distribution in the upper and middle reach basin.

Figure 6 Fig. 6 .
Figure 6 Comparison of simulated results with observed values at Yingluoxia station Fig. 6.Comparison of simulated results with observed values at Yingluoxia station: (a) Calibration period, (b) Validation period 1, (c) Validation period 2.

Figure 7
Figure 7 Comparison of simulated results with observed values at Zhamashike station

Fig. 7 .
Fig. 7. Comparison of simulated results with observed values at Zhamashike station.

Figure 9 Figure 10
Figure 9 Water budgets in the upper reach basin in 1999

Fig. 9 .
Fig. 9. Water budgets in the upper reach basin in 1999.

Figure 9
Figure 9 Water budgets in the upper reach basin in 1999

Figure 11
Figure 11 Spatial variations of hydrological processes and sensible heat flux in the upper and middle reaches basin in 2000

Fig. 11 .
Fig. 11.Spatial variations of hydrological processes and sensible heat flux in the upper and middle reaches basin in 2000: (a) Evapotranspiration, (b) Surface runoff, (c) Groundwater outflow, (d) Sensible heat flux.
. The rivers in the basin include mainstream of the Heihe river, east branch Babao river in the upper reach, Liyuan river in the middle reach and several small rivers in mountainous areas.For model validation, geographic data including geological features, soil vegetation, hydrological and meteorological data and land use data as well as relevant socio-economic data (population and GDP for every Hydrol.Earth Syst.Sci., 13, 1849-1866, 2009 www.hydrol-earth-syst-sci.net/13/1849/2009/ 12/34 model validation, geographic data including geological features, soil vegetation, hydrological and meteorological data and land use data as well as relevant socio-economic data (population and GDP for every county in 1981-2002) and water use data (in 1981-2002 for every county and irrigation area) in the basin are collected and packed up in the GIS platform for use.

Table 2 .
Thermodynamic properties parameters of soil and other mediums.

Table 4 .
Aerodynamic parameters for surface covers.

Table 5 .
Simulation results of monthly discharge at Yingluoxia station.

Table 6 .
Simulation results of daily discharge process at Yingluoxia station.

Table 7 .
Simulation results of water budgets in the middle reach basin from 1982 to 2002 (Unit: mm).Qin=River inflows from Yingluoxia and Liyuanpu, E=Evapotranspiration, Qout=River outflow at Zhengyixia, DS=Storage change, I rr=Irrigation water,R=Total runoff, Rs=Surface runoff, Rg=Groundwater runoff and spring outflow, and DS=P +Qin-E-Qout)

Table 8
gives the annual simulation results of energy budgets in the upper and middle reach basins from 1982 to 2002.

Table 8 .
Simulation results of energy budgets in the upper and middle reach basins from 1982 to 2002 (Unit: MJ/m 2 /year).RN =net radiation, lE =latent heat flux, H =sensible heat flux, and G=soil heat flux.)

Table 9 .
Impacts of conservation forest construction on water budgets in the upper reach basin (Unit: mm/year).

Table 10 .
Impacts of changing farmland to pasture land on water budgets in the middle reach basin (Unit: mm/year).