Interactive comment on “ Frozen soil parameterization in a distributed biosphere hydrological model ” by L .

Parameterization of frozen soil processes in distributed hydrological models is very important for cold region hydrology but researches in this aspect are seldom reported. This paper introduced a new frozen soil parameterization scheme for distributed hydrological modeling and tested it in a cold region alpine watershed with successful results, and thus is an important contribution to the hydrological modeling in cold region. The paper is also clearly presented and well organized


Introduction
Frozen soil (comprising permafrost and seasonally-frozen soil) process is critically important in the land surface hydrology of cold regions, since the freeze-thaw cycle significantly modulates the soil hydraulic and thermal characteristics that directly affect the water and energy cycles in the soil-vegetation-atmosphere transfer (SVAT) system.
At present, the improved modeling of the frozen soil process in a land surface scheme has been recognized an indispensable task for more reliable estimates of soil moisture and temperature profiles, particularly in the winter of Northern Hemisphere.The Project for Intercomparison of Land surface Parameterization Schemes phase 2(d) (PILPS 2(d)) has shown that the models with a frozen soil parameterization generally simulated realistic soil temperature during winter than those without a frozen scheme (Luo et al., 2003).Up to now, many studies have made efforts on improving the frozen soil parameterization in land surface modeling (e.g., Bonan et al., 1996;Slater et al., 1998;Koren et al., 1999;Smirnova et al., 2000;Li and Koike, 2003;Poutou et al., 2004;Woo et al., 2004;Niu and Yang, 2006;Zhang et al., 2007;Nicolsky et al., 2007;Luo et al., 2009).
By contrast, the frozen soil process is often inadequately represented or even neglected in most distributed hydrological models for basin-scale simulations, with only very few exceptions (e.g., Cherkauer and Lettenmaier;1999;Zhang et al., 2000;Stocker-Mittaz et al., 2002;Tian et al., 2006;Mou et al., 2008;Ye et al., 2009).In fact, the representations of the frozen soil process can be indispensable in distributed hydrological modeling for the understanding of the water and energy cycles in some Northern-Hemisphere L. Wang et al.: Frozen soil (c) discretization from a model grid to a number of geometrically symmetrical hillslopes, and (d) description of the water moisture transfer from atmosphere to river.Where R sw and R lw are downward solar radiation and longwave radiation, respectively, H is the sensible heat flux, and λ is the latent heat of vaporization.Here, the land surface submodel is used to describe the transfer of the turbulent fluxes (energy, water, and CO 2 ) between atmosphere and ground surface for each model grid; while the hydrological submodel simulates both surface and subsurface runoff using grid-hillslope discretization, and then simulates flow routing in the river network.river basins, because many large rivers in Northern Hemisphere (e.g., the Yellow River (see Tang et al., 2008) and the Heihe River (see Gao et al., 2008)) originate from high and cold mountain regions.As a result, the hydrological simulations in the upper cold subbasins of these river basins with a spatially-distributed manner become crucial for improving integrated water resources management.
This study aims at improving the performance of a distributed biosphere hydrological model (WEB-DHM; Wang et al., 2009a, b) in the simulations of the frozen soil process, which has been formulated in a previous study (Li and Koike, 2003).This frozen soil scheme will be modified and incorporated into the WEB-DHM for better descriptions of the frozen soil process.By using yearlong continuous observations (soil moisture, soil temperature, and discharge) from the cold region hydrology experiment of WA-TER (Watershed Allied Telemetry Experimental Research; Li et al., 2008Li et al., , 2009)), the newly-developed WEB-DHM with Hydrol.Earth Syst.Sci., 14, 557-571, 2010 www.hydrol-earth-syst-sci.net/14/557/2010/ L. Wang et al.: Frozen soil parameterization in a distributed biosphere hydrological model 559 the revised frozen scheme has been rigorously evaluated in a small cold river basin (Binggou) at the upper area of the Heihe River.

Model description
The WEB-DHM (Water and Energy Budget-based Distributed Hydrological Model) was a distributed biosphere hydrological model, which can give consistent descriptions of water, energy and CO 2 fluxes at a basin scale (see Fig. 1).It can efficiently simulate hydrological processes of largescale river basins while incorporating subgrid topography (see Wang et al., 2009a).
In the study, a new version of the WEB-DHM has been developed by modifying and incorporating a frozen soil parameterization (Li and Koike, 2003).First, the original WEB-DHM is briefly reviewed in Sect.2.1.Details about the hydrological submodel were discussed in Wang et al. (2009a); while the formulations of the land surface submodel can be found in Sellers et al. (1996a).Second, the frozen soil parameterization in WEB-DHM is presented in detail.

General model structure
As illustrated in Fig. 1, the general model structure can be described as follows: 1.A digital elevation model (DEM) is used to define the target basin, which is then divided into sub-basins (see Fig. 1a).Within a given sub-basin, a number of flow intervals are specified to represent time lags and the accumulating processes in the river network.Each flow interval includes several model grids (see Fig. 1b).
2. For each model grid with one combination of land use type and soil type, the land surface submodel independently calculates turbulent fluxes between the atmosphere and land surface (see Fig. 1b and d).The vertical distributions of water for all the model grids, such as ground interception storage and soil moisture profile, can be obtained through this biophysical process.
3. Each model grid is subdivided into a number of geometrically symmetrical hillslopes (see Fig. 1c).A hillslope with unit length is called a basic hydrological unit (BHU) of the WEB-DHM.For each BHU, the hydrological submodel is used to simulate lateral water redistributions and calculate runoff comprised of overland, lateral subsurface and groundwater flows (see Fig. 1c  and d).Overland flow is described by Manning's equation, and lateral subsurface flow and groundwater discharge are simulated using Darcy's law (Wang et al., 2009a).The runoff for a model grid is the total response of all BHUs within it.
4. For simplicity, the streams located in one flow interval are combined into a single virtual channel.All the flow intervals are linked by the river network generated from the DEM.All runoff from the model grids in the given flow interval is accumulated into the virtual channel and led to the outlet of the river basin.The flow routing of the entire river network in the basin is modeled using the kinematic wave approach.

Soil model
Two different soil subdivision schemes are used in describing land surface processes and hydrological processes.
In the calculation of land surface processes, the three-layer soil structure for the unsaturated zone is the same as that in SiB2.The depth of the first layer (D 1 ) is defined as 5 cm, while the root depth (D 1 + D 2 ) could be defined according to vegetation type by SiB2 default.The thickness of the deep soil zone (D 3 ) changes with fluctuation of the water table and is equal to the depth of the groundwater level minus the thickness of the upper two layers.
In the simulation of soil water flow, a multiple-sublayer soil structure is employed to describe the unsaturated zone.In the model, the non-uniform vertical distribution is represented using an assumption of exponentially decreasing hydraulic conductivity with increasing soil depth given by (Beven, 1982;Cabral et al., 1992;Robinson and Sivapalan, 1996), where k surface and k z are hydraulic conductivities at the soil surface and depth z, and f is a decay factor.The surface layer with a depth D 1 is kept as the first layer.The root zone and deep soil zone are uniformly subdivided into several sublayers.As shown in Fig. 2, the multiple-sublayer structure is employed to calculate vertical interlayer flows and lateral runoff.The vertical interlayer flows in the unsaturated zone are described using a one-dimensional Richards equation (see Wang et al., 2009a).
In this soil model, the van Genuchten equation (van Genuchten, 1980) is used as the soil hydraulic function.

Soil hydraulic properties
Total volumetric water content of the j -th soil layer (θ j ; m 3 m −3 ) is defined as Where, θ liq,j and θ ice,j are the liquid water content and the ice content (m 3 m −3 ) of the j -th soil layer; ρ ice and ρ liq are the density (kg m −3 ) of ice and liquid water, respectively.For the j -th soil layer, the unfrozen water content (θ liq,j ) is assumed as a simple power function of soil temperature (T soil,j ) (see Nakano et al., 1982;Romanovsky and Osterkamp, 2000;Li and Koike, 2003)   where, a and b are two empirical coefficients associated with soil type.Then the changing rate of liquid soil moisture (or ice) to soil temperature can be derived as The soil hydraulic conductivity K j and soil matric potential ψ j at the j -th soil layer in the unsaturated zone are described using a modified version of the van Genuchten's equation (van Genuchten, 1980): Where, θ s is porosity and θ r is residual water content; α, n are empirical parameters in van Genuchten's equation with m = 1 − 1 n; the hydraulic conductivity reduction factor for the j -th soil layer (f ice,j ) is defined as a function of soil temperature at that layer f ice,j = exp −10 • T f − T soil,j and 0.05 ≤ f ice,j ≤ 1; (6) and K sat,j is the saturated hydraulic conductivity at depth z (defined as the location of the center of the j -th soil layer), which is measured downward in a direction normal to the soil surface (m).K sat,j is represented using the assumption of an exponential increase in hydraulic conductivity with increasing soil depth (that is the decay factor f < 0).Similar to K j , the groundwater hydraulic conductivity K G is also formulated considering the frozen soil effect Where, the reduction factor f ice,G (= exp −10 • T f − T B and 0.05 ≤ f ice,G ≤ 1) is expressed as a function of the temperature of the bottom soil layer (T B ); K G0 is the groundwater hydraulic conductivity without considering frozen soil (m s −1 ).

Soil thermal properties
For the soil temperatures at the ground surface and the deep soil, the force-restore model (Deardorff, 1977) of the heat balance in the soil surface is kept, but the effective heat capacities of the soil surface and the snow-free soil are modified to represent the latent heat of fusion or the change of soil thermal conductivity where, C g and C d are the effective heat capacity (J m −2 K −1 ) for the soil surface and the snow-free soil; T g and T d are the temperatures for the soil surface and the deep soil, respectively (K); Rn g is the absorbed net radiation by soil surface (W m −2 ); H g is the sensible heat flux from soil surface (W m −2 ); E g is the bare soil evaporation rate (kg m −2 s −1 ); λ is the latent heat of vaporization (J kg −1 ); τ d is daylength (s); ξ gs is the energy transfer due to phase changes in snow on ground (W m −2 ).
The soil thermal conductivity (H s,new ; W m −1 K −1 ) is calculated following Li and Cheng (1995) (11) Hydrol.Earth Syst.Sci., 14, 557-571, 2010 www.hydrol-earth-syst-sci.net/14/557/2010/ L. Wang et al.: Frozen soil parameterization in a distributed biosphere hydrological model 561 Where, H s is soil thermal conductivity without considering frozen soil; H s,max is the maximum heat conductivity after freezing; θ liq,min is the minimum liquid water content after freezing.In the study, θ liq,min = 0.05; while H s,max = 2.5 considering the existence of gravels in the soil.
The new equations of C d and C g are described as follows where d s is the effective depth (m) that feels the diurnal change of temperature (Stull, 1988); C soil is the volumetric heat capacity of soil (J m −3 K −1 ) and ∂T g represents the apparent heat capacity of soil freezing in the surface soil layer, and ∂θ liq,1 water store (m); C w is the volumetric heat capacity of water (J m −3 K −1 ).It should be mentioned that when the soil surface temperature is just below the freezing point, the phase changing rate ∂θ liq,1 ∂T g reaches its maximum value, making the effective heat capacity a very large value and therefore hampering the heat transfer.In order to solve this problem, C g = d s C soil +min 0.05, M gw + M gs C w is set for the transition zone from T f − 0.01 to T f in this study.
The depth of seasonal frost penetration is determined by the soil temperature profile, which is solved with Stefan solution (Yershov, 1990).The Stefan solution assumes a linear soil temperature profile in the frozen soil column, and for simplicity the soil temperature at 5 cm (T soil,D 1 ) is assumed as Where, η is an empirical factor and η = 0.5 is used in this study.
The frost and thaw depths (ζ f and ζ t ) is simulated following Li and Cheng (1995), and can be expressed in equations of the approximation Stefan solution (The Institute of Geocryology, 1974;Li and Koike, 2003) as follows.
where, κ f , and κ t are the thermal conductivity of freezing and thawing soils (W m −1 K −1 ), respectively; τ h is the time length (s) and τ h = 3600 s in the study.
After determining the position of freezing front, the sublayer soil temperature in the root zone and deep soil are estimated by a simple function of frost depth.
If z ≤ ζ f , the soil temperature at a given depth z is If z > ζ f , the soil temperature at a given depth z is where D s is the top soil depth (see Fig. 2).

Datasets for the study area
The cold region hydrology experiment is one of key experiments within the Watershed Allied Telemetry Experimental Research (WATER; Li et al., 2008Li et al., , 2009)).The Binggou watershed (Fig. 3), located at 100 • 12 -100 • 18 E and 38 • 1 -38 • 4 N, is one of the three foci experimental areas where cold region hydrology experiments were carried out within the framework of the WATER.The Binggou watershed, located in the upper reaches of the Heihe River Basin (Fig. 3), is a high mountain drainage system with an area of 30.48 km 2 .The elevation is from about 3440 to 4400 m (Fig. 3).The watershed has the obvious vertical-zonality natural landscape (Yang et al., 1992).In the altitude from 3440 to 4000 m, there is mainly alpine meadow; while in the altitude from 3440 to 3700 m, the shrubs and the fish-scale shape sod coexist.In the region above 4000 m, there is mainly the non-vegetation's alpine desert, and the exposed decency rock debris quite grows having high water-permeability.The longterm mean annual temperature is about −5.8 • C in the mean altitude (3900 m) of the watershed (Yang et al., 1992).In the watershed, the permafrost distributes at the region higher than 4000 m, with the air temperature lower than 0 • C in 9 months (September to next May); while the discontinuous permafrost dominates the region lower than 4000 m, with the air temperature lower than 0 • C in 7 months (October to next April) (Yang et al., 1993).The longterm mean annual precipitation is about 686 mm, with about 74% rainfall and 26% snowfall (Zhang and Yang, 1991).The snowfall prevails from October to April; the rainfall concentrates on July and August; while sleet occurs in May, June and September (Yang et al., 1992).The mean depth of the seasonable snowpack is about 0.5 m, with a maximum of 0.8-1.0m (Yang et al., 1993).
The datasets of the Binggou Watershed, as used in WEB-DHM, are as described below.
DEM and land use were provided by the WATER project (Li et al., 2008(Li et al., , 2009)).The model simulation adopted a grid L. Wang et al.: Frozen soil parameterization in a distributed biosphere hydrological model The precipitation, relative humidity, air temperature, and wind speed, as well as air pressure, downward longwave and shortwave radiation, obtained from the station were summated into the hourly time series and taken as the forcing data for the whole watershed because the basin is small and distributed data is not available.The surface air temperature inputs were modified with a lapse rate of 6.5 K km −1 , considering the elevation differences between the model grids and the DY station.However, the altitudinal effect on relative humidity was assumed negligible.At the basin outlet, the Binggou stream gauge was newly built in 2008 and discharge data has been obtained from 17 January to 20 November 2008 for model evaluation, with the frequency of several times a day.The method of discharge measurements is briefly described as follows.First, a rectangle cross section was built using concrete at the selected point.Second, the water depth (which was used to derive the area of the cross section) and the flow velocity at the cross section were measured by the local staff several times a day (not hourly).Third, the flow velocity and the area of the cross-section were used to calculate the discharge at the cross-section.
The dynamic vegetation parameters are Leaf Area Index (LAI) and the Fraction of Photosynthetically Active Radiation (FPAR) absorbed by the green vegetation canopy, and can be obtained from satellite data.Global LAI and FPAR MOD15 BU 1 km data sets (Myneni et al., 1997) were used in this study.These are 8-daily composites of MOD15A2 products and were obtained using the Warehouse Inventory Search Tool (WIST) of NASA.
All simulations were carried out with 250 m spatial resolution and hourly time step.

Model evaluations at the Binggou watershed
For the Binggou watershed, the hydrological processes are not only controlled by the hydro-meteorological conditions of the land surface, but also by the underlying frozen soil.The spring snowmelt occupies about 30% of the annual runoff; while the residual comes from rainfall and groundwater (Zhang and Yang, 1991).In the lower area of the watershed, soil starts thawing around April and stops by August; while in the upper regions, soil starts thawing around late May.In the watershed, the thickness of the active frozen soil layer is about 1.0-1.5 m in the lower regions, and is greater than 3.0 m in the upper mountain regions (Yang et al., 1993).
Hydrol.Earth Syst.Sci., 14, 557-571, 2010 www.hydrol-earth-syst-sci.net/14/557/2010/In the early spring (April to May), with the increase of air temperature, snowmelt occurs from the lower regions to the mountain areas.However, during this period (April to May), the air temperature exceeds 0 • C only at noon, but drops to below 0 • C at night.Consequently, much of the snowmelt water freezes again at night before its departure from snowpack.Therefore, the snowmelt runoff in the early spring is small (around 15% of annual runoff; Zhang and Yang, 1991).From May to June (late spring), the air temperature increases to above 0 • C stably, and the snowmelt runoff becomes very large (greater than 25% of annual runoff; Zhang and Yang, 1991).This is also attributed to the little permeability of the underlying seasonal frozen soil layers which have thawed only in upper soil layers (Kane and Stein, 1983).In summer, snow and seasonal frozen soil layers disappear, and thus rainfall becomes the major source for river discharges.But the permanent frozen soil layers still exist, which prohibit the water infiltration to deeper layers.Heavy rainfall events in summer will usually result in severe flash floods in this watershed, along with landslides and debris flows (Yang et al., 1993).

Model calibration
The vegetation static parameters including morphological, optical and physiological properties were initially defined following Sellers et al. (1996b).In the summer 2008, by using the WEB-DHM without the frozen scheme, the land surface parameters were optimized using the observed radiation fluxes at the DY station in July; and then the two van Genuchten parameters (α and n) were optimally obtained by the calibration of the July soil moistures at upper layers (5, 10 and 20 cm depths) at the DY station.In the cold season from 21 November 2007 to 20 April 2008, by using the WEB-DHM with the frozen scheme, the parameters (a and b) used in the frozen soil scheme were optimized through matching the simulated and observed soil temperatures at 5 cm depth (T soil,D 1 ).Third, by using the WEB-DHM with the frozen scheme, the other soil hydraulic parameters were optimally obtained by the calibration of the discharges at the basin outlet in July and August that covers the annual largest flood peak in 2008.

Parameters optimized through the WEB-DHM without the frozen scheme
First, the land surface parameters were optimized using the observed surface radiation fluxes at the DY station in July.
For the DY station, land is covered by the SiB2 biome 9 ("Agriculture/C3 grassland").The canopy cover fraction, the height of canopy top, and the height of canopy bottom, as well as the root depth (Dr), the top soil depth (Ds), and the ground roughness length (z s ) have been designed as 0.2, 0.05 m, 0.005 m, 0.25 m, 1.25 m, and 0.001 m according to the field obervations in Li et al. (2009).The soil re-flectance to visible radiation for the SiB2 biome 9 was optimized as 0.15 using observed upward shortwave radiation in July 2008; while the surface emissivity has been kept as 1.0 following Seller et al. (1996b).Other time-invariant vegetation parameters were set following Sellers et al. (1996b).Soil hydraulic properties have been kept equal to values derived from FAO (2003) during the calibration of land surface parameters.
Second, the two van Genuchten parameters (α and n) were optimized by using the measured soil moistures at upper layers at the DY station.The porosity and the residual soil moisture were set as 0.585 and 0.017 according to the yearlong soil moisture measurements from 21 November 2007 to 20 November 2008 in the DY station.The van Genuchten parameters α and n regulate the soil hydraulic function which controls the soil water transport.They were optimized as 0.1 and 2.1, respectively, by comparing the simulated and observed soil moistures at the upper soil layers (5, 10, and 20 cm) in July 2008 at the DY station; while keeping the original values of K surface obtained from FAO (2003).

Parameters optimized through the WEB-DHM with the frozen scheme
First, at the DY station, the frozen soil parameters (a and b) are optimized through matching the simulated and observed T soil,D 1 in the cold season (from 21 November 2007 to 20 April 2008); while d s is set as 0.6 m, according to the measured winter soil temperature profiles at the DY station (Li et al., 2009).Second, at the basin-scale, the other soil hydraulic parameters were further optimized to obtain good reproduction of the flood event that occurred at the Binggou stream gauge during the summer (July and August) of 2008.These parameters include the saturated hydraulic conductivity for soil surface K surface , the soil anisotropy ratio anik, the groundwater hydraulic conductivity (without considering frozen soil) K G0 , and the hydraulic conductivity decay factor f .The optimization was done using a trial and error method by matching the simulated and observed flood peaks and tails.It should be mentioned that the assumption of an exponential increase in hydraulic conductivity with increasing soil depth is used in the WEB-DHM with the frozen scheme (f < 0).
The basin-averaged values of the land surface and soil hydraulic parameters, as well as the frozen soil parameters used in the Binggou watershed are listed in Table 1.

Calibration results
The bias error (BIAS) and root mean squared error (RMSE) are used as evaluation criterion for the simulated results, where BIAS and RMSE are defined as www.hydrol-earth-syst-sci.net/14/557/2010/ Hydrol.Earth Syst.Sci., 14, 557-571, 2010 where x oi is the observation, x si is the simulation, and N is the total number of time series for comparison.
Figure 4 shows the simulated and observed hourly upward solar radiation from 1 to 31 July 2008 at DY station with the BIAS and the RMSE, by using the WEB-DHM without the frozen scheme.The observed soil moisture and temperature profiles were used to initialize the model at the first hour of 1 July 2008; while the initial water table depth was assumed as same as the initial depth of the unsaturated zone (Ds=1.25 m).After the calibration of the soil reflectance to visible radiation, the diurnal cycles of the upward solar radiation are well represented by the calibrated WEB-DHM without the frozen scheme.The BIAS and RMSE for the simulated upward shortwave radiation at the DY station are −3.8W m −2 and 32.6 W m −2 , respectively.It should be mentioned that the measurements of upward longwave radiation in the station were found erroneous for all periods, and was not used for model evaluation in the study.
Figure 5 illustrates the hourly evolutions of precipitation and the simulated and observed hourly volumetric liquid soil moisture at 5, 10, 20, 40, 80, and 120 cm in July 2008 at the DY station, by using the WEB-DHM without the frozen scheme.Reasonable responses of soil moisture at the upper layers (5, 10, and 20 cm) to the rainfall events are reproduced with high accuracies (Fig. 5b-d).The BIAS for the simulated soil moistures at 5, 10, and 20 cm are −0.003,0.011, and −0.012; while their RMSE values are 0.031, 0.029 and 0.030, respectively.The soil moisture at 120 cm is also accurately simulated by the WEB-DHM without the frozen scheme (see Fig. 5g), since the soil at this depth was still in frozen in July at the station.The RMSE values for the  simulated soil moisture at 40 and 80 cm (particularly 80 cm) are relatively large, which can be attributed to the thawing activity.It reveals that even in the summer simulations at the DY station, a frozen soil scheme is essential.Figure 6a gives the hourly simulated and observed soil temperature at 5 cm in the DY station from 21 December 2007 to 20 April 2008, by using the WEB-DHM with the frozen scheme.After the calibration of frozen soil parameters (a and b), the soil temperature at 5 cm in the cold seasons is well reproduced, with the BIAS and RMSE values equal to 0.84 K and 1.68 K, respectively.
Figure 6b plots the calibrated hourly hydrograph from July to August 2008 at the Binggou gauge, by using the WEB-DHM with the frozen scheme.Here, the WEB-DHM with the frozen scheme was initialized from a previous simulation starting from 21 November 2007.It is shown that after calibration, in general the model reproduces both the peak and base flows well.The difference of timing and magnitudes between the observed and simulated peak flows is possibly attributed to the sparse density of the meteorological sites used in the study (only the DY station).It should be mentioned that the discharges at the Binggou gauge were measured discontinuously and irregularly, and thereby the evaluation criterions (e.g., the Nash-Sutcliffe model efficiency coefficient (Nash and Sutcliffe, 1970) and the bias error are not estimated.

Model validation
The calibrated WEB-DHM with the frozen scheme was then used for a yearlong simulation from 21 November 2007 to 20 November 2008, to check its performance.It is obvious that the WEB-DHM with the frozen scheme predicts the snow depth much better than the WEB-DHM without the frozen scheme.7 2007.11.21 2008.2.19 2008.5.19 2008.8.17 2008.1117 2008.3.17 2008.5.16 2008.7.15 2008.9.13 Discharge (m  17 2008.3.17 2008.5.16 2008.7.15 2008.9.13 Discharge (m November 2007 to 20 November 2008, simulated by the WEB-DHM without and with the frozen scheme.Results showed that the WEB-DHM with the frozen scheme generally gives more realistic yearlong (in particular the thawing periods) soil moisture profile than those simulated by the WEB-DHM without the frozen scheme.For the soil moisture at surface layer (0-5 cm), root zone (5-25 cm), and deep soil layer (25-125 cm), the results by the WEB-DHM with the frozen scheme obtains the RMSE of 0.068, 0.043, and 0.088, respectively; while the ones by the WEB-DHM without the frozen scheme gets the RMSE of 0.109, 0.068, and 0.086, respectively.The overestimation of soil moisture at deep soil layer (25-125 cm), from May to July 2008 by the WEB-DHM with the frozen scheme is possibly attributed to the large gravels distributed in deep soil layers.Without considering the large gravels, the excessive recharges to the deep soil layer from the unconfined aquifer may have been simulated with the WEB-DHM with the frozen scheme.The noises in the simulated liquid soil moistures are caused by the changes of the soil temperatures from below freezing to above freezing.seasons) hourly discharges at Binggou gauge are further evaluated.

Discharges at the
Figure 11 displays the hourly hydrographs simulated by the WEB-DHM without and with the frozen scheme at the Binggou gauge from 17 January to 20 November 2008.Results show that the streamflows simulated by the WEB-DHM with the frozen scheme agree fairly well with the observations; while those calculated by the WEB-DHM without the frozen scheme have large difference from the observations, with the overestimation of baseflows (January-April), the underestimation of snowmelt flows (April-May), and the overestimation of peak flows (July-August).The improved simulation results by the WEB-DHM with the frozen scheme can be attributed to the following reasons: 1.The consideration of reduction factor of the groundwater hydraulic conductivity (f ice,G ) largely improves the model's performance from January to April (also see Fig. 12). Figure 12b demonstrates the logarithmic hourly hydrographs simulated by using the WEB-DHM with the frozen scheme at the Binggou gauge from 17 January to 20 November 2008, which further confirms the good performance of the new system in simulating base flows.Without the treatments of the frozen soil effect on the groundwater hydraulic conductivity, excessive baseflows have been obtained from January to April 2008 (see Fig. 12a).
2. In the spring season (April-May), more snowmelt runoff calculated by the WEB-DHM with the frozen scheme is due to the treatments of the hydraulic conductivity reduction for each soil layer (f ice,j ) as a function of soil temperature, which has resulted in less infiltration during the snow melting.
3. The poorer flood peaks obtained by the WEB-DHM without the frozen scheme, are possibly caused by the overestimation of the soil hydraulic conductivities at the deep soil layer and the unconfined aquifer during the thawing process.

Concluding remarks
In this study, the distributed biosphere hydrological model WEB-DHM was improved by incorporating a frozen soil parameterization.The WEB-DHM with the frozen scheme was then applied to the Binggou watershed for evaluation using the in-situ observations from WATER.
After calibrating land surface parameters, soil hydraulic parameters, and frozen soil parameters, the WEB-DHM with the frozen scheme was then used for a yearlong validation from 21 November 2007 to 20 November 2008, to check the model's applicability in the continuous simulation.Results show that the WEB-DHM with the frozen scheme has given much better performance than the WEB-DHM without the frozen scheme, in the simulations of soil moisture profile at the DY station and the discharges (base and peak flows as well as snowmelt runoff) at the basin outlet in the yearlong simulation.
In the literatures regarding to the frozen soil parameterization in land surface modeling, there are mainly three ways to simulate the ice content or the unfrozen water content within the soil.These methods estimate the ice content depending on the heat available (e.g., Slater et al., 1998;Takata and Kimoto, 2000;Dai et al., 2003), or using the freezing point depression equation (e.g., Koren et al., 1999;Smirnova et al., 2000;Niu and Yang, 2006;Zhang et al., 2007), or relying on empirical equations (Pauwels and Wood, 1999;Li and Koike, 2003).The empirical equations-based frozen soil parameterization by Li and Koike (2003) is used in the study, to improve the original WEB-DHM for simulating frozen soil dynamics, largely due to its simple structure and low computation costs.Furthermore, the study by Li and Koike (2003) has shown the good compatibility between the empirical frozen soil parameterization and SiB2; while the original WEB-DHM uses a hydrologically-improved SiB2 (Wang et al., 2009c) as its land surface submodel.Although the method in the study estimates the ice content based on empirical equations, the calibrated WEB-DHM with the frozen scheme has demonstrated acceptable accuracies in simulating point-scale frozen soil dynamics and basin-scale integrated streamflows.Meanwhile, the frozen scheme using empirical equations also results in difficulty in achieving an Hydrol.Earth Syst.Sci., 14, 557-571, 2010 www.hydrol-earth-syst-sci.net/14/557/2010/ optimal set of the frozen soil parameters a and b, as the new model (the WEB-DHM with the frozen scheme) is very sensitive to the two parameters (particularly b) in simulating the soil temperature and moisture profiles.Different from Li and Koike (2003) that formulated frozen soil process in a 1-D land surface model (SiB2), this study has modified and incorporated the frozen soil scheme into a distributed biosphere hydrological model (WEB-DHM).The newly-developed WEB-DHM with the frozen scheme has made it possible to simulate the basin-scale cold-region land surface hydrological processes in a spatially-distributed manner while considering the topographically-driven lateral flows.It can be used as a model operator in the catchmentscale land surface hydrological data assimilation system in cold regions, to improve modelling of soil moisture and the surface energy budget as well as streamflows.

Symbol Definition
C d effective heat capacity for the snow-free soil (J m −2 K −1 ) C g effective heat capacity for the soil surface (J m −2 K −1 ) C soil volumetric heat capacity of soil (J m −3 K −1 ) C w volumetric heat capacity of water (J m −3 K −1 ) d s effective depth that feels the diurnal change of temperature (m) E g bare soil evaporation rate (kg m −2 s −1 ) f hydraulic conductivity decay factor f ice,G reduction factor for the groundwater hydraulic conductivity f ice,j hydraulic conductivity reduction factor for the j -th soil layer H g sensible heat flux from soil surface (W m −2 ) H s soil thermal conductivity without considering frozen soil H s,new soil thermal conductivity considering frozen soil K G groundwater hydraulic conductivity (m s −1 ) K G0 constant groundwater hydraulic conductivity without considering frozen soil effect (m s −1 ) K j soil hydraulic conductivity at the j -th soil layer (m s −1 ) K sat,j saturated hydraulic conductivity for the j -th soil layer (m s −1 ) K surface saturated conductivity at the soil surface (m s −1 ) L f latent heat of fusion (J kg ice content of the j -th soil layer (m 3 m −3 ) θ j total volumetric water content of the j -th soil layer (m 3 m −3 ) θ liq,j liquid water content of the j -th soil layer (m 3 m −3 ) θ r residual soil water content (m 3 m −3 ) θ s porosity (m 3 m −3 ) λ latent heat of vaporization (J kg −1 ) ξ gs energy transfer due to phase changes in snow on ground (W m −2 ) ρ ice density of ice (kg m −3 ) ρ liq density of liquid water (kg m −3 ) ψ j soil matric potential at the j -th soil layer Fig. 1.Overall structure of WEB-DHM: (a) division from basin to sub-basins, (b) subdivision from sub-basin to flow intervals comprising several model grids, (c) discretization from a model grid to a number of geometrically symmetrical hillslopes, and (d) description of the water moisture transfer from atmosphere to river.Where R sw and R lw are downward solar radiation and longwave radiation, respectively, H is the sensible heat flux, and λ is the latent heat of vaporization.Here, the land surface submodel is used to describe the transfer of the turbulent fluxes (energy, water, and CO 2 ) between atmosphere and ground surface for each model grid; while the hydrological submodel simulates both surface and subsurface runoff using grid-hillslope discretization, and then simulates flow routing in the river network.

Fig. 2 .
Fig. 2.Soil model of the WEB-DHM.Two different soil subdivision schemes are used for describing land surface and hydrological processes, respectively.The three-layer soil structure used in SiB2 is retained to represent the unsaturated zone in the calculations of land surface processes.The unsaturated zone is divided into multiple sublayers when simulating water flows within it and water exchanges with groundwater aquifer.

Fig. 4 .
Fig. 4. Simulated and observed hourly upward shortwave radiation (unit: W m −2 ) at the DY station in July 2008, by using the WEB-DHM without the frozen scheme.

Fig. 5 .
Fig. 5. Hourly precipitation (a), and the simulated and observed hourly volumetric liquid soil moisture at 5, 10, 20, 40, 80, and 120 cm (b-g) at the DY station in July 2008, by using the WEB-DHM without the frozen scheme.

Fig. 6 .
Fig. 6.The hourly simulated and observed soil temperature at 5 cm in the DY station from 21 December 2007 to 20 April 2008 (a), and the calibrated hydrograph at the Binggou gauge in July and August 2008 (b), by using the WEB-DHM with the frozen scheme.

Fig. 7 .Fig. 8 .
Fig. 7. Hourly observed (interpolated from the observations of soil temperature profile) and simulated thaw depth at the DY station, from 1 May to 31 August 2008, by using the WEB-DHM with the frozen scheme.

Fig. 9 .FigureFigure 9
Fig. 9. Hourly snow depth at the DY station from 21 November 2007 to 20 November 2008, simulated by the WEB-DHM without and with the frozen scheme.Here, the snow depth is assumed as five times of the snow-water equivalent, and the large fluctuations of the observed snow depth at the station were caused by strong wind blowing.

Fig. 11 .
Fig. 11.Linear hourly hydrographs simulated by the WEB-DHM without and with the frozen scheme at the Binggou gauge from 17 January to 20 November 2008.

Fig. 12 .
Fig. 12. Logarithmic hourly hydrographs simulated by using the WEB-DHM with the frozen scheme, neglecting (a) and considering (b) the frozen soil effect on the groundwater hydraulic conductivity at the Binggou gauge from 17 January to 20 November 2008.

Table 1 .
Basin-averaged values of the parameters used in the Binggou watershed.

station from 1 May to 31 August 2008
Figure 7 displays the hourly observed (interpolated from the observations of soil temperature profile) and simulated thaw depth at the DY station, from 1 May to 31 August 2008, by using the WEB-DHM with the frozen scheme.It is shown that in general the model estimates the thaw depth with an acceptable accuracy, with the BIAS and RMSE values equal to −0.15 m and 0.20 m, respectively.