Improving the snow physics of WEB-DHM and its point evaluation at the SnowMIP sites

In this study, the snow physics of a distributed biosphere hydrological model, referred to as the Water and Energy Budget based Distributed Hydrological Model (WEB-DHM) is significantly improved by incorporating the three-layer physically based energy balance snowmelt model of Simplified Simple Biosphere 3 (SSiB3) and the Biosphere-Atmosphere Transfer Scheme (BATS) albedo scheme. WEB-DHM with improved snow physics is hereafter termed WEB-DHM-S. Since the in-situ observations of spatially-distributed snow variables with high resolution are currently not available over large regions, the new distributed system (WEB-DHM-S) is at first rigorously tested with comprehensive point measurements. The stations used for evaluation comprise the four open sites of the Snow Model Intercomparison Project (SnowMIP) phase 1 with different climate characteristics (Col de Porte in France, Weissfluhjoch in Switzerland, Goose Bay in Canada and Sleepers River in USA) and one open/forest site of the SnowMIP phase 2 (Hitsujigaoka in Japan). The comparisons of the snow depth, snow water equivalent, surface temperature, snow albedo and snowmelt runoff at the SnowMIP1 sites reveal that WEB-DHM-S, in general, is capable of simulating the internal snow process better than the original WEB-DHM. Sensitivity tests (through incremental addition of model processes) are performed to illustrate the necessity of improvements over WEB-DHM and indicate that both the 3-layer snow module and the new albedo scheme are essential. The canopy effects on snow processes are studied at the Hitsujigaoka site of the SnowMIP2 showing that the snow holding capacity of the canopy plays a vital role in simulating the snow depth on ground. Through these point evaluaCorrespondence to: M. Shrestha (maheswor@hydra.t.u-tokyo.ac.jp) tions and sensitivity studies, WEB-DHM-S has demonstrated the potential to address basin-scale snow processes (e.g., the snowmelt runoff), since it inherits the distributed hydrological framework from the WEB-DHM (e.g., the slope-driven runoff generation with a grid-hillslope scheme, and the flow routing in the river network).


Introduction
Seasonal snow cover is an important component of land surface hydrology and is critical for simulation of water and energy budgets in cold climate regions.Snow with its high albedo, low roughness, relatively low thermal conductivity and considerable spatial and temporal variability, can greatly alter energy and water interactions among the atmosphere, vegetation and land.Snow has the ability to store and release water within the hydrological cycle.The appearance of snow cover may lead to a temporal shift in the runoff during the spring snowmelt period and is a significant parameter from the view of hydrological simulation.
To understand and represent the snow processes in land surface modeling, a large number of approaches have been used in many land surface schemes (LSSs) in diversified numerical expressions, ranging from simple degree -day models to physically based sophisticated multi-layer energy balance models (Brun et al., 2008).Many numerical studies have been carried out to develop and validate snow submodels of different complexity in LSSs of many climate and hydrological models (e.g., Verseghy, 1991;Blöschl et al., 1991;Douville et al., 1995;Tarboton and Luce, 1996;Yang et al., 1997;Loth and Graf, 1998a, b;Marks et al., 1999;Jin et al., 1999a, b;Sun et al., 1999;Sud and Mocko, 1999;Essery et al., 1999;Smirnova et al., 2000;Mocko and Sud, 2001;Sun and Xue, 2001;Xue et al., 2003;Yang and Niu, 2003;Dai et al., 2003;Zanotti et al., 2004;Sun and Chern, 2005;Liston and Elder, 2006;Hirai et al., 2007, Ellis et al., 2010;Dutra et al., 2010).Several snow-scheme intercomparison studies have been undertaken to gain an improved understanding of snow cover simulation in LSSs and to address issues related to the current state of snow modeling used by the atmospheric and hydrologic research community (e.g., the Project for the Intercomparison of Land-Surface Parameterization Schemes (PILPS) -Phase 2d (Slater et al., 2001) and Phase 2e (Bowling et al., 2003), the Snow Model Intercomparision Project Phase 1 (SnowMIP1; Etchevers et al., 2004) and Phase 2 (SnowMIP2; Rutter et al., 2009;Essery et al., 2009) and the Rhône-Aggregation LSS Intercomparison Project (Boone et al., 2004)).Many studies showed that snow accumulation processes were well represented by single-layer snow models but diurnal freeze and thaw cycles were not well captured by these models, resulting in errors in the simulation of snow surface temperature and snow melting in terms of timing and the total amount (Lynch-Stieglitz, 1994;Sun et al., 1999, Jin et al., 1999b;Slater et al., 2001;Luo et al., 2003;Xue et al., 2003), raising the importance of the development and application of multilayer energybalance-based snow models.On the other hand, uncertainties in the forcing data and initial conditions would have great impact in the snow process simulations while comparing the performances among different complexity of snow models (e.g., Slater et al., 2001;Feng et al., 2008).
At present, multilayer energy-balance snow parameterizations are mainly employed by the 1-D LSSs that are generally used in climate models.To our knowledge, only very few distributed hydrological models (DHMs) have included such sophisticated energy-balance snow schemes for studying the cold region processes (e.g., Cherkauer and Lettenmaier, 1999;Zanotti et al., 2004).In fact, a DHM with a multilayer energy-balance snow module can physically describe the snow accumulation, ablation and snowmelt runoff which are critically important for both the current water resources management practices and the climate change adaptation studies in cold and high mountain river basins.This study discusses the improvement of snow physics in a distributed biosphere hydrological model named as Water and Energy Budget based Distributed Hydrological Model (WEB-DHM; Wang et al., 2009a, b, c).WEB-DHM is developed by fully coupling Simple Biosphere 2 (SiB2; Sellers et al., 1996) with a hillslope hydrological model (Yang et al., 2002(Yang et al., , 2004)).It can realistically simulate the land surface and hydrological processes and provide consistent descriptions of water, energy and CO 2 fluxes at a basin scale.The snow physics of WEB-DHM is improved by adopting the ideas derived from the studies of different snow models and by incorporating the three-layer snow physics of Simplified Simple Biosphere 3 (SSiB3; Sun and Xue, 2001;Xue et al., 2003).SSiB3 is developed by coupling SSiB (Xue et al., 1991) with a three-layer version of the Simple Atmosphere-Snow Trans-fer (SAST; Sun et al., 1999) and has been successfully applied to simulate snow processes in cold regions (Xue et al., 2003;Durand and Margulis, 2006;Walisher et al., 2009).WEB-DHM with improved snow physics is hereafter termed WEB-DHM-S.Since the in-situ observations of spatiallydistributed snow variables with high resolution are currently not available over large regions, the new distributed system (WEB-DHM-S) is at first rigorously tested with comprehensive point measurements.This evaluation data comprise the observational datasets from four open sites of the SnowMIP1 (Col de Porte in the French Alps, Weissfluhjoch in the Swiss Alps, Goose Bay in Canada, and Sleepers River in USA) and one open/forest site of the SnowMIP2 (Hitsujigaoka in Japan).

Model description
A short review of the snow processes in WEB-DHM is given in Sect.2.1, while the snow processes in WEB-DHM-S are discussed in detail in Sect.2.2.Details of the hydrological and land surface submodels of WEB-DHM can be found in Wang et al. (2009a) and Sellers et al. (1996).

Snow processes in WEB-DHM
In WEB-DHM, the parameterization of the snow submodel is the same as that for SiB2 (Sellers et al., 1996).A single-layer bulk snow mass balance is considered with constant density (200 kgm −3 ), and the thermal regime of snow is not distinguished from that of soil.Attenuation of downward shortwave radiation through the canopy is considered with multiple scattering between the canopy and snow/ground but attenuation of radiation within the snow layer is ignored.Only the top 5 cm of the snow water equivalent is considered for variation of the heat capacity of the surface skin, which affects the surface energy balance in the case of a large snow mass.The snow surface temperature is represented by the average snowpack temperature, which tends to result in incorrect simulation of the surface energy budget, which in turn affects the overall accumulation and melting processes.Moreover, it does not consider the prognostic snow albedo.The dry snow albedo is given as a constant value of 0.8 for visible (VIS) shortwave radiation and 0.4 for near infrared (NIR) shortwave radiation.For melting snow, the snow albedo is simply set to 60% of the dry snow albedo.

Snow processes in WEB-DHM-S
In this section, the energy and mass budget equations along with snow parameterization are presented in detail.In WEB-DHM-S, the snow parameterizations for the canopy are kept the same as in WEB-DHM, but the single-layer snow scheme on the ground is replaced by the SSiB3 snow scheme when the snow depth is greater than 5 cm.Initially, the snowpack is divided into three layers that start with the same initial snow Hydrol.Earth Syst.Sci., 14, 2577Sci., 14, -2594Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/2577/2010/ temperatures.The top layer thickness is kept at a fixed depth of 2 cm regardless of the total snow depth to provide reasonable simulation of the diurnal changes in the snow surface temperature.The maximum thickness of the middle layer is kept at 20 cm, and the bottom layer represents the remaining body of the snowpack.A surface energy balance equation is formulated only for the top layer, which is influenced by the surface radiation budget and sensible and latent heat fluxes.The heat budget of the second and third layers is controlled by the heat conduction and the penetrating shortwave radiation.Over time, these three layers evolve differently through their own energy budgets and the heat exchanges between them.
Meanwhile, the mass budget for each layer is calculated accordingly by taking account of the precipitation, evaporation/condensation, compaction, liquid water retention, snowmelt runoff and infiltration into the underlying layers.When snow melts, meltwater in a layer increases, thereby increasing the layer-average density and mass.Any meltwater in a layer exceeding the liquid water holding capacity is delivered to the underlying layer.Water leaving the bottom snow layer is available for partitioning into soil water infiltration and/or surface runoff by the soil-vegetation-atmosphere transfer (SVAT) system.This snow scheme can produce a variable density profile.
The snow-covered surface albedo scheme is parameterized using a physically based prognostic snow albedo scheme of the Biosphere-Atmosphere Transfer Scheme (BATS) model (Dickinson et al., 1993;Yang et al., 1997), and the snow cover fraction is calculated using the formulations of Mocko and Sud (2001).Major differences between the snow processes in WEB-DHM and WEB-DHM-S are presented in Table 1.The soil model coupled with a three layer snow model in WEB-DHM-S is shown in Fig. 1.

Energy balance equations
The energy content of the snowpack is affected by the shortwave radiation penetration, heat conduction between sublayers, ground heat fluxes, the flux of advection due to precipitation, energy due to phase change and net radiation at the surface accompanied by sensible and latent heat fluxes.Specific enthalpy is used as the prognostic variable instead of snow temperature in the energy balance equation, which includes the internal energy of liquid water or ice as well as the energy of the phase change.It is assumed that liquid water at its melting point has zero enthalpy so that the phase change processes can be tackled easily.The same approach was also employed by Lynch-Stieglitz (1994), Tarboton and Luce (1996), Jin et al. (1999a), Sun et al. (1999)   Xue (2001).The energy budget equation for the canopy is the same as that in WEB-DHM.However, the canopy temperature T c is influenced by the snow surface enthalpy.The energy budget equation for the canopy is where C c (Jm −2 K −1 ) is the effective heat capacity for the canopy, R nc ,H c and λE c (Wm −2 ) are net radiation, sensible heat flux and latent heat flux for the canopy respectively, and ξ c (Wm −2 ) is the energy transfer due to phase changes in the canopy.The equation for enthalpy of each snow layer is where H (Jm −3 ) is the volumetric enthalpy of water, Z j is the snow depth of layer j and G sn (Wm −2 ) is the heat flux through the snow layer.H and G sn are defined as where R nsn (Wm −2 ), H sn (Wm −2 ), λE sn (Wm −2 ), G pr (Wm −2 ), K (Wm −1 K −1 ), T sn (K) and SW sn (Wm −2 ) are net radiation, sensible heat, latent heat flux, thermal energy from rain at the snow surface, thermal conductivity of snow, snow temperature and shortwave radiation flux absorbed by the snow layer respectively.Turbulent and radiative fluxes are calculated using the formulations of SiB2 (Sellers et al., 1996) except that the snow surface temperature is used instead of the average bulk snow temperature for the surface energy balance.f ice is the dry-snow mass fraction of the total mass in the snow layer, and h v (Jkg −1 ) is the latent heat of fusion for ice.C v (Jm −3 K −1 ) is the mean snow volumetric specific heat capacity, parameterized as a function of the bulk density of snow (ρ s ; kgm −3 ) and intrinsic density of ice (ρ i ; kgm −3 ) following Verseghy (1991): The thermal conductivity of snow K (Wm −1 K −1 ) is adopted from Jordan (1991).
where K i (2.29 Wm −1 K −1 ) and K a (0.023 Wm −1 K −1 ) are the thermal conductivities of ice and air respectively.The penetration of shortwave radiation flux into the snow layers is accounted for in this model.Hence, the shortwave energy available for the surface energy budget is completely different from that in WEB-DHM.The shortwave radiation SW sn at the snow layer is defined following Jordan (1991): where SW nsn = SW sntop (1-α s ).SW sntop (Wm −2 ) is the radiation incident on the snow surface, α s is snow albedo and β vis and β nir are extinction coefficients; β vis = 0.003795d −1/2 ρ s (Z j ) and β nir = 400.The grain size diameter d (m) is specified as a function of density following Anderson (1976).Thermal energy from rain (G pr ) can be calculated as where IF 0 (ms −1 ) is the infiltrated flux rate of rain at the snow surface, T rain (K) is the temperature of rainfall, ρ w (kgm −3 ) and C w (Jkg −1 K −1 ) are the density and specific heat capacity of water.For simplicity, T rain is considered as air temperature.Ground surface temperature (T g ) and deep soil temperature (T d ) are obtained by considering conductive heat flux at the snow/soil interface and the forcerestore model (Deardorff, 1978) of the heat balance in the soil surface.
where C g and C d are the effective heat capacity (Jm −2 K −1 ) for the soil surface and deep soil, τ d is the day length (s) and K(Z 1 ) is the effective thermal conductivity at the snow/soil interface.The prognostic equations of snow surface enthalpy and canopy temperature are solved simultaneously by calculating the temperature increments for the physics time step using an implicit backward numerical scheme.T sn (Z 2 ), T sn (Z 1 ), T g and T d are the variables with slow change which are solved explicitly using a forward numerical scheme.The final equations for solving T c and T sn (Z 3 ) are represented as where K eff (Wm −2 K −1 ) is the effective thermal conductivity of snow between the top and the middle snow layer and M snow (m) is the snow water equivalent (SWE).K eff is defined as Eqs. ( 11) and ( 12) are solved using "one step test method".Initially, it is tested by assuming that current state of the snow surface layer is in frozen state completely (f ice = 1) and then T c and T sn (Z 3 ) will be solved.The test is true if T sn (Z 3 ) is less than the freezing temperature (T f ).Otherwise, the state will be in either partially melted or completely melted state.In this case, T sn (Z 3 ) is assumed to T f and f ice is solved.The surface layer is in partial melting state if 0 < f ice < 1.The snow is melted completely if f ice < 0. In this case, the solutions for f ice and T sn (Z 3 ) should be f ice = 0 and T sn (Z 3 ) = T f .

Mass balance equations
The mass balance equation for the canopy is the same as in WEB-DHM.The mass balance for snow is represented by the change in liquid water and ice content in the snowpack.The relative change in snow mass is controlled by snowfall/rainfall, compaction, snow melting, runoff, infiltration into the underlying snow layer/soil and evaporation/sublimation at the snow surface.Neglecting the effect of water vapor diffusion and its phase change to mass distribution, the mass balance equations for the snow layer are other layers (j = 2, 1), (14) where M snow,j (m) corresponds to the SWE at snow layer j , P s (ms −1 ) is the rate of snowfall, IF j (ms −1 ) = min (O j ,P avs ), is the actual liquid water infiltration flux at the interfaces, R j (ms −1 ) is runoff from the lower interface and E sn (ms −1 ) is the combined evaporation and sublimation rate.O j is the liquid water outflow rate that will be drained to the underlying layer if the total liquid water in layer exceeds its liquid water holding capacity (C r ).Liquid snow mass fraction, f liq = (1-f ice ) is used to calculate the total amount of liquid water.P avs is the pores available in the layer.R j is calculated as the difference between IF j and O j .The liquid water holding capacity (C r ) is taken as a function of the snow layer density following Anderson (1976): where C rmin = 0.03, C rmax = 0.1, γ e = 200 kgm −3 and γ i (kgm −3 ) is bulk density of ice.The bulk density of ice for new snowfall is calculated following the formulation used in the CROCUS snow model (Brun et al., 1989;Brun et al., 1992): 50 , (16)   where T air is the air temperature (K) and u m is the wind speed (ms −1 ).

Snow compaction
Three snow compaction processes, namely destructive metamorphism, densification due to snow overburden and compaction due to snow melting, are included.The compaction process is critically important for the evolution of density and snow depth.The snow depth is decreased by the compaction and is increased by snowfall.These three components of snow compaction are parameterized following Anderson (1976).The empirical equation for destructive metamorphism is where γ i (kgm −3 ) and γ l (kgm −3 ) are bulk densities of ice and liquid water and C 3 and C 4 are empirical constants.After snow has undergone its initial settling stage, densification due to overburden proceeds at a slower rate.This compaction rate is a function of snow overburden pressure W s (Nsm −2 ), such that where η o (3.6 × 10 6 Nsm −2 ) is the viscosity coefficient, C 5 = 0.08 K −1 and C 6 = 0.023 m 3 kg −1 .The decrease in thickness of the snow sublayer due to melting is estimated as where h i is the dry-snow mass in a unit depth and dh i is the dry-snow mass that melts in the unit depth.Hence, total compaction over one time step is given by The rate of change in snow density caused by snow compaction is given by

Snow albedo
The snow albedo is parameterized using a physically based prognostic snow albedo scheme of the BATS model (Dickinson et al., 1993;Yang et al., 1997).The albedo is computed for VIS and NIR spectral bands with adjustments for illumination angle and snow age.The total snow albedo (α s ) is the weighted average of VIS and NIR albedos, which depends on the spectral ratio of the incident shortwave radiation.VIS and NIR albedos (α vis , α nir ) are defined as where α vd and α nird are the albedos of the diffused shortwave radiation in the VIS and NIR bands respectively, α vis0 (0.95) and α nir0 (0.65) represent fresh-snow albedos for the VIS and NIR bands, f zen is the correction term for a solar zenith angle larger than 60 • and f age is the snow aging factor accounting for the effect of grain growth due to vapor diffusion and the effect of dirt and soot.The snow albedo parameterization is very sensitive to α vis0 and α nir0 .These fresh-snow albedos can be parameterized depending upon the snow type and characteristics of the site.Details of f zen and f age can be found in Dickinson et al. (1993), Yang et al. (1997).

Dataset
Dataset for evaluation of models include four open sites of the SnowMIP1: Col de Porte (CDP), Weissfluhjoch (WFJ), Goose Bay (GSB) and Sleepers River (SLR).In addition, one open/forest site of the SnowMIP2: Hitsujigaoka (HSG) is selected for forest snow processes evaluation.Meteorological forcing data includes hourly air temperature, relative humidity, wind speed, precipitation amount, the snow/liquid fraction, downward shortwave radiation (DSR) and downward longwave radiation (DLR).The vegetation coverage parameter is set to zero for simulation at all the SnowMIP1 sites.Details about data and site characteristics are discussed here and a summary is given in Table 2.

Col de Porte (1996-1998)
CDP is a mid-range elevation site at 1340 m above mean sea level (a.m.s.l.), located in the northern French Alps (45.3  Goodison et al. (1998) and its phase was determined based on an air temperature relationship derived from comparisons of the Geonor gauge with snowfall observations.Evaluation data comprise hourly observations of snow surface temperature from a downwardlooking radiometer, hourly observations of snow depth from an ultrasonic sensor supported by weekly snow course observations of the SWE and snow depth, and the daily total of bottom runoff from a 5 m 2 lysimeter protected from lateral flow.Data from this site have been used to evaluate many SVAT snow schemes (e.g., Brun et al., 1992;Douville et al., 1995;Loth and Graf 1998a;Sun et al., 1999;Essery et al., 1999;Sun and Xue, 2001;Boone and Etchevers, 2001;Strasser et al., 2002;Xue et al., 2003;Essery and Etchevers, 2004;Brown et al., 2006;Li et al., 2009).

Weissfluhjoch (1992-1993)
WFJ is a high-elevation site at 2540 m a.m.s.l. with flat topography, located in the eastern Swiss Alps (46.83 • N, 9.81 • E) and managed by the Swiss Federal Institute for Snow and Avalanche Research.The average air temperature during the period of continuous snow cover is −2.9 • C. Rainfall does not occur from mid-October to mid-May.Snow continuously accumulates from mid-October until mid-April and then melts through May and June owing to temperatures above the melting temperature.Although this site is windier than CDP, drifting and blowing effects are weaker (Essery and Etchevers, 2004;Brown et al., 2006).Evaluation data comprise hourly observations of snow surface temperature from an infrared thermometer, hourly observations of snow depth from an ultrasonic sensor supported by weekly and sometimes biweekly snow pit observations of the SWE and snow depth, daily snow albedo and daily snowmelt runoff.Data from this site have been used in the assessment of many snow models (e.g., Fierz and Lehning, 2001;Lehning et al., 2002;Fierz et al., 2003;Essery and Etchevers, 2004;Etchevers et al., 2004;Brown et al., 2006).
The 15 years forcing and validation data do not correspond to the same site.Hourly air temperature, humidity, wind speed and precipitation were measured at the GSB airport site.Radiation measurements were made at the GSB Upper Air station located at the east end of the airport (53.30• N 60.37 • W).Incoming longwave radiation was estimated using observations of hourly air temperature, relative humid-ity, cloud type and opacity following Idso (1981) and Sellers (1965).Hourly precipitation rate data was derived from 6-hourly precipitation totals observed with a Niphershielded gauge corrected for wind-induced undercatch, wetting loss and trace precipitation amounts following Metcalfe et al. (1997) and Goodison et al. (1998).Mean daily temperature ranges from −16.4 • C in January to 15.8 • C in July, with a mean annual total snowfall of 434 mm.Daily snow depth observations were made manually using a ruler at the GSB airport site.The site is humid and is relatively windy compared to the other SnowMIP sites.Potential blowing snow conditions were encountered approximately 10% of the time during the December to April period (Brown et al., 2006).Data from this site have been used in the assessment of many snow models (e.g., Bélair et al., 2003;Essery and Etchevers, 2004;Brown et al., 2006;Gordon et al., 2006).

Sleepers River (1996-1997)
SLR is a low-elevation site at 552 m a.m.s.l., located in the northeastern Vermont (44.50 • N, 72.17 • W).The site is characterized by almost flat topography surrounded by northern hardwood forest.The average air temperature during the period of continuous snow cover is −4.5 • C. The snow cover lasts from the beginning of November to the end of April.
Precipitation was separated into snow and rain as a linear function of air temperature, with precipitation assumed to be all rain at temperatures above 2 • C and all snow below 0 • C. Evaluation data comprise hourly observations of snow depth from an ultrasonic sensor supported by weekly and sometimes biweekly snow pit observations of the SWE and snow depth, and daily snowmelt runoff.Snowmelt runoff data is not used for evaluation due to some uncertainties associated with the lysimeter data.Blowing effect is not seen this year as the wind speed is too low.Data from this site have been used in the assessment of many snow models (e.g., Anderson, 1976;Lynch-Stieglitz, 1994;Albert and Krajeski, 1998).
Center for the Hokkaido Region (Open site), about 2500 m away from the forest site.Precipitation rate is corrected for wind-induced undercatch (Yokoyama et al., 2003) and is partitioned between snow and rain following the approach of Yamazaki (2001) using the wet bulb temperature.Only the snow depth measurements are made available for the evaluation.Data from this site have been used in the study of canopy snow influence on water and energy balance above coniferous forest (e.g., Nakai et al., 1999a, b;Suzuki and Nakai, 2008, Rutter et al., 2009, Dutra et al., 2010).

Simulation results
The performance of the model is evaluated by comparing the simulated and observed SWE, snow depth, snow surface temperature, snow density, snow albedo and snowmelt runoff.
The bias error (BIAS) and root mean square error (RMSE) are used as evaluation criterion for the simulated results and are defined as where Xsim i and Xobs i are simulated and observed values at a given time step for n paired simulation and observation values.

Snow depth, SWE and snow density
Simulation results and observations of the snow depth, SWE and snow density for the two seasons at CDP, one season at WFJ and SLR are shown in Fig. 2. For WEB-DHM-S, the snow depth along with the realistic simulation of snow density is well reproduced at all sites whereas WEB-DHM is unable to capture the variability of snow depth because it assumes a constant snow density.At the CDP site in 1996-1997 (see Fig. 2a), the SWE is underestimated by WEB-DHM in the beginning of the accumulation season which undertakes its impact throughout the snow season.WEB-DHM-S capture the accumulation season well but the mid winter ablation in late January is found not enough to meet the observations.The possible reason may be due to the uncertainty in the precipitation phase.At the end of the melting season, WEB-DHM underestimates the SWE because of the early melting and the low snow albedo.In 1997-1998 (see Fig. 2b), the ablation prevailed at mid-March causing the continuous decrease in the SWE and the SWE is increased to about 0.25 m with significant snowfall in mid-April.Although both models are able to simulate the snow accumulation process well, the results show that the SWE is overestimated by both models in the mid-season from late January to mid-February.This overestimation is due to the failure in capturing the rapid decrease in the SWE during 21 January to 27 January.The uncertainty in the precipitation phase is not a major in this case.The reason may be the rapid increase in the albedo due to dust and fallen leaves.After the mid season ablation, the SWE is well simulated by WEB-DHM-S.In the mean time, the SWE simulated by WEB-DHM is very low as compared to the observations due to the low albedo and the increased rate of snowmelt also.
The results of statistical analysis of the simulation results are presented in Table 3.At the WFJ site, snow coverage lasts from mid-October to late June (see Fig. 2c).The results show that the SWE is underestimated by WEB-DHM in the accumulation season owing to the strong melt simulation in early November, and all the snow has melted by mid-May, whereas the SWE simulated by WEB-DHM-S during accumulation seasons is found to be in good agreement with the observed SWE.The snow depth simulated by WEB-DHM-S is found to be remarkably underestimated from early April to early June.Statistical analysis shows that WEB-DHM has less BIAS than WEB-DHM-S (see Table 3) but it does not mean that the WEB-DHM results are good.Indeed, there is a large overestimation of snow depth by WEB-DHM from January to mid-April and a large underestimation from mid-April to late June.At the SLR site, WEB-DHM-S overestimated the SWE, snow depth and snow density throughout the snow season as shown in Fig. 2d.This bias may be due to the misrepresentation of the precipitation phase as the total precipitation is divided into the rain and snow as a linear function of air temperature as discussed in the Sect.3. The results for the SWE simulation during the melting season for WEB-DHM are found better than that for WEB-DHM-S.
The time-slice evaluation of the model in simulating the first, maximum, minimum during the mid season, one prior to last and last SWE observations at the CDP, WFJ and SLR sites are presented in Table 4.For the maximum SWE observation in 1996-1997 at the CDP site, the results show that the SWE is slightly overpredicted by WEB-DHM-S (BIAS = 0.023) whereas largely underpredicted by WEB-DHM (BIAS = −0.069).Likewise, the minimum SWE observation during the mid season is simulated in the similar trend of the maximum SWE simulation.Both the models underpredicted these variables at CDP in 1997-1998.In this year, the maximum SWE is well simulated by both the models but WEB-DHM is found to have the large BIAS (−0.057) in simulating the minimum SWE at the mid season as compared to the BIAS (−0.003) for WEB-DHM-S.At the WFJ site, WEB-DHM has very large BIAS (−0.325) in simulating the maximum SWE.Overall, WEB-DHM-S shows better performance than WEB-DHM in simulating the maximum SWE, minimum during the mid season and one prior to the last SWE observations at the CDP and WFJ sites.However, the performance of WEB-DHM is better than WEB-DHM-S in simulating all these parameters at SLR.
Figure 3 shows the result for 15 year simulation at the GSB site in simulating the interannual variability of the snow depth.Both the models are found to be capable in multiyear simulations.The correlation coefficient for WEB-DHM and WEB-DHM-S is found to be 0.68 and 0.78, respectively.The snow depth is largely overestimated in most of the years.The reasons for this bias are not well understood.However this site is affected by blowing snow condition but Gordon et al. (2006) showed that the results are not much improved by incorporating the blowing snow physics in CLASS model.The comparison of the SWE and snow density are excluded in this study as the snow course measurements are made in a sparsely wooded area 4 km away from the snow depth measurement site.
The results for the snow density as shown in Fig. 2 reveal that WEB-DHM-S is able to capture the trend of the seasonal variation in the snow density.At the CDP site in 1996-1997, the snow density is well simulated throughout the snow season whereas in 1997-1998, it is overestimated in the mid-season (during mid February) owing to the overestimation of snowmelt.At the end of the melting season, the observed snow density has increased to 450 kgm −3 but the model fails to simulate this event owing to underestimation of the SWE during this period.The model output shows similar characteristics at the WFJ site.At SLR, the snow density is overestimated throughout the snow season due to  GSB (1969GSB ( -1984)). the overestimation of the snow depth and the SWE.In general, WEB-DHM-S is found to simulate the variability in the snow depth, SWE and snow density more accurately than WEB-DHM.

Snow surface temperature
Snow surface temperature is an important parameter of the land surface energy balance as it plays a vital role in the estimation of exchanges of moisture and heat fluxes between the snow surface and atmosphere.The results for the snow surface temperature simulations at CDP (1996( -1997( ), CDP (1997( -1998) ) and WFJ (1992WFJ ( -1993) ) are shown in Fig. 4a, b and c, respectively.The results indicate that the simulation performance of WEB-DHM-S is significantly improved as compared to that of WEB-DHM.WEB-DHM has large RMSE and BIAS (see Table 3) because the snow surface temperature is calculated as the averaged temperature for a single bulk layer of snow mass, and thus, the nighttime surface temperature is overestimated.It is found that the RMSE at CDP (1996)(1997) considerably reduce from 3.68 in WEB-DHM to 2.72 in WEB-DHM-S.In 1997-1998, the RMSE for WEB-DHM and WEB-DHM-S are found 3.22 and 2.07, respectively.At WFJ, there is a reduction of RMSE from 5.70 to 3.10 while employing WEB-DHM-S.The observed snow surface temperature at WFJ is available up to 3 May 1993 only whereas continuous snow cover exists till 30 June 1993.The statistical values of BIAS and RMSE for WEB-DHM at this site will increase if we analyze the results for the whole snowy period because snow melts out too early in the simulation of WEB-DHM.The results show that WEB-DHM-S still has some cold bias during the night at CDP while the model has warm bias during the day and night at WFJ (see Fig. 4a, b, c).The warm bias is due to the underestimation of snow albedo whereas the cold bias is associated with the deficiency in Monin-Obukhov similarity theory to calculate the turbulent fluxes in a highly stable condition and the uncertainty in the roughness length of the snow surface.However, the simulations can be improved by the inclusion of windless coefficient as discussed in Brown et al. (2006).

Snow albedo
The snow albedo observed at the WFJ site is used in the model evaluation.There are also snow albedo observations for the CDP site but they are not used in this study as the CDP albedo is underestimated owing to partial obstruction of the sensor's field of view (Etchevers et al., 2004).Fresh snow albedo in the VIS band is calibrated with a factor of 0.95 for the WFJ site and 0.87 for the CDP site.Figure 5 compares the observed daily mean albedo and the simulation results of WEB-DHM and WEB-DHM-S.The simulation results show that WEB-DHM-S is able to capture the seasonal evolution of snow albedo; however, there is a strong bias of 0.1 to 0.15 during the accumulation period, and thus, the results obtained are identical to those obtained using the CLASS model and those available through the SnowMIP (Essery and  (1996)(1997) from 11 November 1996 to 3 April 1997.Etchevers, 2004;Etchevers et al., 2004;Brown et al., 2006).The main reason behind this bias is that the observed albedo for new snow is around 0.95 whereas the simulated maximum albedo is 0.84.which contributes to the surface runoff and infiltration to the soil surface.The timing and total amount of snowmelt runoff is better simulated by WEB-DHM-S than by WEB-DHM.

Snowmelt runoff
In two seasons at CDP, WEB-DHM-S is found to capture the snowmelt runoff during the accumulation season, midablation season and final melting season.Although the results for WEB-DHM also show similar runoff behavior, they include large biases (underestimation) during the accumulation and final melting season in both years.Due to early melt for WEB-DHM, the runoff is overestimated from the beginning to the end of March in 1996-1997 and from the beginning to the middle of April in 1997-1998 (see Fig. 6a, b) (see Fig. 6a, b).At the WFJ site, the observations of snowmelt runoff are available only for a short period (27 April to 7 July 1993) and the simulation results of WEB-DHM-S have far better agreement with the observed runoff pattern than the simulation results of WEB-DHM.A large amount of snowmelt runoff is simulated by WEB-DHM during early April to early May owing to the early melting in the case of WEB-DHM (see Fig. 6c).A substantial improvement in snowmelt runoff simulation is achieved at both sites by WEB-DHM-S with less RMSE and BIAS (see Table 3).Comparison of simulated daily totals of snowmelt runoff with the available observed values at (a) CDP (1996CDP ( -1997)), (b) CDP (1997CDP ( -1998) ) and (c) WFJ (1992WFJ ( -1993)).

Sensitivity to incremental process representation
Different sets of simulations (see Table 5) are carried out at CDP and WFJ to illustrate the sensitivity to incremental process representation.WEB-DHM is taken as the control run (CTRL) simulation.At CDP, WEB-DHM with realistic albedo (CTRL A) improves the SWE simulation in the accumulation season but it fails to simulate the SWE in the melting season due to the overestimation of albedo (see Fig. 7a, b).The snow season is overpredicted by 35 days in 1996-1997and 17 days in 1997-1998. At WFJ, it (1996)(1997), (b) CDP (1997CDP ( -1998) ) and (c) WFJ (1992WFJ ( -1993) ) for different sets of simulations as shown in Table 5. CTRL A WEB-DHM + BATS albedo scheme CTRL B WEB-DHM + 3 layer snow scheme CTRL C WEB-DHM + 3 layer snow scheme + realistic albedo CTRL D WEB-DHM + 3 layer snow scheme + BATS albedo scheme = WEB-DHM-S NEW simulate the SWE in the accumulation and initial melting period whereas the SWE in the late melting season is well simulated (see Fig. 7c).It shows that the constant realistic albedo parameterization to WEB-DHM without its decay function is not able to improve the simulation capability at all.The inclusion of BATS albedo scheme into WEB-DHM (CTRL B) is able to improve the performance of WEB-DHM in simulating the SWE and albedo at both the CDP and WFJ sites but the snow depth is still poorly simulated due to the lack of prognostic simulation of snow density.Furthermore, the total amount and timing of the snowmelt runoff also has large biases as compared to the observations in CTRL B simulation (see Fig. 8).This indicates that the single layer snow model alone is not enough to simulate the overall process well.
WEB-DHM combined with 3 layer snow scheme (CTRL C) mainly improves the snow depth simulation.This test shows that the overall simulation results at the CDP site are improved in 1996-1997.But in the melting season of 1997-1998, the SWE results are worse than that of CTRL B due to low albedo in CTRL C. At WFJ, the SWE in the accumulation season is slightly improved but the snow is melted too early as occurred in CTRL.CTRL C with increased albedo value (CTRL D) lengthens the snow cover days due to the high albedo without decay function in melting season.This implies that both, the 3 layer snow physics and the new albedo scheme are critically important.At the end, WEB-DHM-S (NEW in Table 5 and Fig. 7) incorporates both 3 layer snow scheme and BATS prognostic albedo scheme for accurate simulation of overall snow processes.

Effect of canopy on snow processes
The effect of canopy on snow processes is evaluated at the Hitsujigaoka (HSG) forest site of the SnowMIP2.Only WEB-DHM-S is used for simulation as snow parameterization for the canopy is kept the same as in WEB-DHM.The model is run blindly using its default parameters for the needleleaf-evergreen trees following SiB2 (Sellers et al., 1996).The zero plane displacement height (0.63 h) and roughness length for the canopy (0.13 h) is taken from Suzuki et al. (2008) where h is the vegetation height  CDP (1997CDP ( -1998) ) and (c) WFJ (1992WFJ ( -1993) ) for different sets of simulations as shown in Table 5.
Only the snow depth measurements are available for model evaluation and the observed snow depth is given as the averaged value of 28 stake measurements in this area.Firstly, the simulation is conducted using the default Sat-cap1 with the value of 0.3 mm.The result as shown in Fig. 9 demonstrates that the snow depth is overestimated from the beginning to the end of January.This bias is mainly due to low value of Satcap1 as Suzuki and Nakai (2008) reported the maximum daily Satcap1 with the value of 6.9 mm.Thus, further simulations are carried out for Satcap1 with the values of 3 mm and 6 mm to test the impact of canopy interception over the snow processes beneath the canopy.With increase in Satcap1 (3 mm and 6 mm), the snow depth under the canopy is reduced throughout the snow season (see Fig. 9).For Satcap1 with the value of 6 mm, the snow depth in accumulation season (mainly in January) is improved but is underestimated from the mid of February to the end of snow season.The evaporation from the canopy snow is found to be 8.3%, 19.2% and 25.2% of the total precipitation for Satcap1 with the values of 0.3 mm, 3 mm and 6 mm, respectively.The RMSE value for the snow depth simulation is increased from 0.093 to 0.108 and 0.124 while increasing Satcap1 from 0.3 mm to 3 mm and 6 mm.The reason for poor performance of the model for realistic Satcap1 is unclear, may be due to an error in the precipitation amount and its phase.We set Satcap1 to the value of 6 mm to simulate the response of the snow depth to the vegetation coverage (Vcover).It is found that the snow depth increases for the decreased Vcover (see Fig. 10).It is obvious that decrease in Vcover causes less interception by the canopy and more snow falls to the ground surface increasing the snow albedo beneath the canopy.But this may not follow at every site and more sites should be validated beforehand to draw solid conclusions.Currently, the model does not include mass releases from the canopy due to melt drip and drop of the snow due to the strong winds and bending of branches which may enhance the poor performance of the model.
In the mean time, we simulated the snow depth at the open site with meteorological data obtained at NAHRC.Since radiation data are not available at NAHRC, they are obtained from the forest site (see Table 2).Figure 11 shows the observed snow depth at the open and forest sites compared with the snow depth simulated by WEB-DHM-S at the open site.From the beginning to the end of the mid ablation season (late February), the snow depth at the open site was observed with high peaks compared to that at the forest site.Afterward, the variability of the snow depth at these two sites is found quite different.Snow is melted out too early at the open site as compared to that at the forest site.The quite contrast may be due to the variability in the precipitation amount and its phase.The simulated result shows that the snow depth at the open site is highly overestimated after early March and thus the snow cover days are overpredicted by 15 days.The model also fails to capture the maximum snow depth.The reasons for these biases are unclear, may be due to the high sublimation or the problem with the forcing data (especially the radiation measurements) and the reasons will be carefully investigated in future.

Conclusions
This study has presented the improvements in the snow physics of WEB-DHM by incorporating a three-layer physically based energy balance snowmelt model of SSiB3 and the BATS albedo scheme.WEB-DHM with improved snow physics is termed WEB-DHM-S.The three-layer snow model in WEB-DHM-S adds more features to the original WEB-DHM to simulate the snow processes more accurately.The snow processes include the variability of snow density, snow depth and SWE, liquid water and ice content in each layer, prognostic snow albedo, diurnal variation in the snow surface temperature, thermal heat due to conduction and liquid water retention.
Datasets from four open sites (CDP, WFJ, SLR and GSB) of the SnowMIP1 and one open/forest site (HSG) of the SnowMIP2 were used for model evaluation.The simulation results of snow depth, SWE, surface temperature and snowmelt runoff revealed that WEB-DHM-S is capable of simulating the internal snow process more accurately than the original WEB-DHM.Snow albedo is better parameterized in WEB-DHM-S than in WEB-DHM.Although WEB-DHM-S is capable of capturing an albedo trend similar to that observed, it still has a strong bias of 0.1 to 0.15 in the albedo value during the accumulation period and hence needs the improvements of the albedo scheme to account for the effect of snow type and dynamic evolution of grain size.Different sensitivity tests are performed to understand the effect of incremental process representations in the model.It is found that both the schemes (the 3-layer snow scheme and the BATS albedo scheme) are critically important for improving the WEB-DHM.The canopy effect on snow processes is studied at Hitsujigaoka site of the SnowMIP2 showing that the snow holding capacity of the canopy plays a vital role in simulating the snow depth on ground.More forest sites will be evaluated in future studies for detailed understanding of the forest snow processes.Through these point evaluations and sensitivity studies, the WEB-DHM-S has demonstrated the potential to address basin-scale snow processes (e.g., the snowmelt runoff), since it inherits the distributed hydrological framework from the WEB-DHM (e.g., the slope-driven runoff generation with a grid-hillslope scheme, and the flow routing in the river network).In next studies, the WEB-DHM-S can be further coupled with a frozen soil scheme (e.g., Wang et al., 2010) and a glacier model to improve the integrated water resources management in cold and high elevated river basins.

Fig. 1 .
Fig. 1.The soil model coupled with a three-layer snow model as described in WEB-DHM-S.

Figure 6 Fig. 4b .
Figure6compares the observed snowmelt runoff and simulation results of WEB-DHM and WEB-DHM-S at the CDP and WFJ sites.Although the snowmelt runoff measurements for the CDP site are available for the whole simulation period, the runoff comparison is made for the snow season only.The total snowmelt is computed as the sum of melt in each layer

Fig. 9 .
Fig. 9. Comparison of observed and WEB-DHM-S simulated snow depth for different maximum canopy snow storage (Satcap1) at the HSG forest site.

Fig. 10 .
Fig. 10.The response of WEB-DHM-S simulated snow depth to different vegetation cover (vcover) at the HSG forest site.

Fig. 11 .
Fig. 11.Comparison of WEB-DHM-S simulated snow depth at HSG open site with observed snow depth at HSG open and forest sites.

Table 1 .
Major differences of snow processes in WEB-DHM and WEB-DHM-S.

Table 2 .
Meteorological characteristics of study sites.
-1997 whereas the snow cover lasts to the beginning ofMay in 1997May in  -1998..Winter air temperatures are not particularly low and rainfall can occur at anytime during the snow season.The site is not windy and is relatively humid.Precipitation was measured with a Geonor gauge with correction for undercatch following • N, 5.77 • E) and managed by Météo-France.The site is characterized by flat topography with loamy soil covered with short grass.The soil generally does not freeze.Continuous snow cover is recorded from the end of November to the beginning of April in 1996 • N, 142.38 • W).The site is mostly flat with sandy soil.It has a cool temperate climate and the snowpack is maritime type.The average air temperature during the period of continuous snow cover is −0.96 • C. Vegetation includes approximately 7m high todo fir.Vegetation coverage is set to 100% for simulation and effective leaf area index (LAI) is set to 3. Canopy snow was present most of the time from the middle of December 1997 to the middle of March 1998 and snow beneath the canopy was present from the middle of Decem- ber 1997 to the middle of March 1998(Suzuki and Nakai,  2008).Radiation measurements were taken at roof of the research center, about 500 m away from the forest site.Precipitation was measured at the National Agricultural Research www.hydrol-earth-syst-sci.net/14/2577/2010/ Hydrol.Earth Syst.Sci., 14, 2577-2594, 2010

Table 3 .
BIAS and RMSE for the SnowMIP1 sites.

Table 4 .
BIAS in simulating the first, maximum, minimum during the mid season, one prior to last and last SWE observations at the CDP, WFJ and SLR sites.

Table 5 .
Different sets of simulations.