Correcting basin-scale snowfall in a mountainous basin using a distributed snowmelt model and remote-sensing data

. Adequate estimation of the spatial distribution of snowfall is critical in hydrologic modelling. How-ever, this is a well-known problem in estimating basin-scale snowfall, especially in mountainous basins with data scarcity. This study focuses on correction and estimation of this spatial distribution, which considers topographic effects within the basin. A method is proposed that optimises an altitude-based snowfall correction factor (C fsnow ) . This is done through multi-objective calibration of a spatially distributed, multilayer energy and water balance-based snowmelt model (WEB-DHM-S) with observed discharge and remotely sensed snow cover data from the Moderate Resolution Imaging Spectroradiometer (MODIS). The Shufﬂed Complex Evolution–University of Arizona (SCE–UA) automatic search algorithm is used to obtain the optimal value of for minimum cumulative error in discharge and snow cover simulations. Discharge error is quantiﬁed by Nash– Sutcliffe efﬁciency and relative volume deviation, and snow cover error was estimated by pixel-by-pixel analysis. The study region is the heavily snow-fed Yagisawa Basin of the Upper Tone River in northeast Japan. First, the system was applied to one snow season (2002–2003), obtaining an optimised C fsnow of 0.0007 m − 1 . For validation purposes, the optimised C fsnow was implemented to correct snowfall in 2004, 2002 and 2001. Overall, the system was effective, implying improvements in correlation of simulated versus observed discharge and snow cover. The 4 yr mean of basin-average snowfall for the corrected spatial snowfall distribution was 1160 mm (780 mm before correction). Execution of sensitivity runs against other model input and parameters indicated that C fsnow could be affected by uncertainty in shortwave radiation and setting of the threshold air temperature parameter. Our approach is suitable to correct snowfall and estimate its distribution in poorly gauged basins, where elevation dependence of snowfall amount is strong.


Introduction
Solid precipitation (snowfall) is of great importance in mountain snow hydrology, since snow acts as a natural reservoir by storing water in winter and releasing it in spring. Snowmelt discharge from mountain snowpack is an important source of energy for hydropower in the low-flow season and water for agriculture and biodiversity maintenance on local and regional scales. With its intrinsic radiative (high albedo) and thermal (low thermal conductivity) properties, snow can strongly modulate energy and water interactions between the atmosphere and land surface. The considerable spatiotemporal variability of snow distribution at basin scale is important in determining the timing and magnitude of spring snowmelt discharge. Such variability can increase the probability of droughts and snowmelt runoff-induced floods. Hence, accurate prediction of discharge during snowmelt season is imperative to support optimal water resource planning and Published by Copernicus Publications on behalf of the European Geosciences Union. management (Singh and Singh, 2001;Armstrong and Brun, 2008).
Towards better representation and accurate simulation of basin-scale snow processes, many physically based singleor multi-layer energy balance distributed snowmelt models have been developed (e.g. Blöschl et al., 1991;Garen and Marks, 2005;Liston and Elder, 2006;Letsinger and Olyphant, 2007;Shrestha et al., 2012a;Mahat and Tarboton, 2012). Successful parameterization of the physical processes can be achieved by the energy balance snow models reducing calibration efforts and allowing intra-basin transfer of knowledge. However, even though a model is physically based, it may still produce inaccurate results if the forcing is incorrect (e.g. Beven, 2004;Garen and Marks, 2005). Precipitation has the greatest uncertainty among all forcings for distributed hydrological model (e.g. Andréassian et al., 2001;Beven, 2004;Bárdossy and Das, 2008;Moulin et al., 2009), and uncertainty is greater for snowfall than for rainfall. This higher uncertainty in snowfall is caused by effects of multiple factors like wind, topography, blowing and drifting, wetting, and evaporation losses at point scale (Sevruk, 1982;Goodison et al., 1998;Fortin et al., 2008), as well as the wind distribution and orographic dependencies at basin scale (WMO, 1986;Milly and Dunne, 2002;Xia and Xu, 2007;Valery et al., 2010). Inconsistency in basin precipitation and snowmelt discharge is observed because of uncertainty in the snowfall distribution (Milly and Dunne, 2002;Lohmann et al., 2004;Fekete et al., 2004;Tian et al., 2007;Yang et al., 2009;Valery et al., 2010;Bartolini et al., 2011).
Several correction methods (Goodison et al., 1998;Adam and Lettenmaier, 2003;Yang et al., 2005;Fortin et al., 2008) have been developed to overcome systematic errors in snowfall measurements at point scale. Some studies were designed to avoid both systematic and non-systematic (site specific) biases associated with snow gauges (Cherry et al., 2005). The method of Cherry et al. (2005) uses observed snow depth and a physics-based land surface model to solve an inverse problem for snowfall. It reconstructs snowfall by calculating the snowfall that must have occurred to produce the observed snow depth, given the physics of the model. However, all these methods deal with snowfall measurement errors at point scale. The areal distribution of snowfall is still a major problem when extending methods based on a point scale to distributed snowmelt modelling, because of insufficient gauge density across a watershed and methods of interpolating point data (Fassnacht et al., 2003;Valery et al., 2009Valery et al., , 2010. There have been a few studies regarding the correction of snowfall at basin scale. Valery et al. (2009) obtained corrected precipitation (both snowfall and rainfall) as a solution of an inverse problem of the hydrologic cycle at a daily time step, which minimised the difference of observed and simulated discharge volume through a simple water balance formula. Bartolini et al. (2011) followed the same concept to attain monthly precipitation, using a temperature index-based snow process model. However, the use of river discharge data alone is inadequate to correct and estimate the spatial distribution of snowfall, since discharge represents only the integrated response of the catchment water balance. These water balance-based methods therefore have limited applicability to predominantly snow-fed basins, since snow accumulation and ablation are intricately linked not only to water balance but also energy balance.
Along with discharge data, the integration of basin-scale snow properties (e.g. observed snow depth, snow water equivalent, and snow covered area) in a hydrologic model could be one approach to account for spatial snow dynamics. Such an integrated dataset may be referred to as "soft data" for internal model process verification on spatial scales (Seibert and McDonnell, 2002). Basin-scale in situ observations of snow properties are limited and difficult to conduct. Consequently, satellite-derived snow covered area (SCA) may be regarded as a relatively reliable snow product or index for representing large-scale snow variability, and may represent the most effective soft data in hydrologic modelling for quantifying the spatial distribution of snowfall in poorlygauged mountainous river basins.
In recent years, the Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover product has been widely used in multidisciplinary studies, owing to its high spatiotemporal resolution (daily and 8 day product on a 500 m grid) and high accuracy relative to snow depth observations at basin to regional scale (e.g. Klein and Barnett, 2003;Parajka and Blösch, 2006;Pu et al., 2007;Wang et al., 2008;Parajka and Blösch, 2012). MODIS snow data have been exploited as forcing for snowmelt runoff models (e.g. Li and Williams, 2008;Immerzeel et al., 2009;Tahir et al., 2011), as a tool in model evaluation (e.g. Shamir and Konstantine, 2007;Bavera and Michele, 2009), as model input in data assimilation schemes (e.g. Clark et al., 2006;Zaitchik and Rodell, 2009), as integrated soft data in calibrating conceptual models (e.g. Parajka and Blösch, 2008;Sorman et al., 2009;Franz and Karsten, 2013), and for reconstructing the spatial distribution of snow water equivalent using distributed hydrologic models (Molotch and Margulis, 2008).
Unlike prior studies regarding use of the MODIS snow cover product and river discharge in calibration and evaluation of lumped and distributed hydrologic models, this work focuses on a new approach to correction of snowfall in basin scale. This approach uses a comprehensive multilayer, energy balance-based snowmelt model, with both water and energy balance closure in a snow-soil-vegetationatmosphere transfer-based distributed biosphere hydrologic modelling framework (WEB-DHM-S; Shrestha et al., 2010Shrestha et al., , 2012a. The method optimises an elevation-dependent correction factor, using a heuristic algorithm called Shuffled Complex Evolution-University of Arizona (SCE-UA; Duan et al., 1992). First, error between the simulated and observed discharge at basin outlet and error between the simulated and MODIS-derived snow cover pixels in the basin are computed Hydrol. Earth Syst. Sci., 18, 747-761, 2014 www.hydrol-earth-syst-sci.net/18/747/2014/ for model simulation with initial arbitrary value of the snowfall correction factor. The objective function is defined to attain the least cumulative error on discharge and snow cover pixels. Model is re-executed by altering the snowfall correction factor till the minimisation of the objective function is achieved. The corrected snowfall obtained following this approach is regarded as the snowfall that would have likely occurred, given the physics of the model. The study basin is the humid Yagisawa Basin of the Upper Tone River in northeast Japan, where water use in spring is completely dependent on snowmelt discharge since contribution of ground water to spring discharge is zero. Water supply from this basin in spring covers about 45 % of the annual water usage. Thus, it is extremely important to obtain precise water resource forecasts, via input of the correct snowfall amount. This paper is organised as follows. Section 2 briefly discusses the materials and methods related to the hydrologic model, input and evaluation data, and an overview of methodological framework for basin-scale snowfall correction. Section 3 demonstrates how discharge and SCA simulation are improved after application of the snowfall correction. Uncertainty attributable to other model input and model parameters is discussed in this section. Conclusions with remarks on potential application of the methods are presented in Sect. 4.

WEB-DHM-S model description
The model used is the Water and Energy Budget-based Distributed Hydrological Model with improved snow physics (WEB-DHM-S; Shrestha et al., 2010Shrestha et al., , 2012a, which was developed by coupling the three-layer energy balance snow scheme of the Simplified Simple Biosphere 3 model (SSiB3; Xue et al., 2003) and the prognostic albedo scheme of the Biosphere Atmosphere Transfer Scheme (BATS; Yang et al., 1997) with the Water and Energy Budget-based Distributed Hydrological Model (WEB-DHM; Wang et al., 2009a, b). The model runs on an hourly time step and at predefined grid size (500 m in this study). The model consists of the three-layer snow routine, multilayer soil routine, and groundwater flow routine for nine land-use categories (open/forest regions) according to Simple Biosphere Model version 2 (SiB2; Sellers et al., 1996). The snow energy balance algorithm uses specific enthalpy as the prognostic variable, which includes both the internal energy of liquid water or ice and energy of phase change. Exchange of mass and energy fluxes with the atmosphere occur at the surface snow layer only, whereas conductive fluxes dominate energy and mass transport within underlying snow layers. The mass balance for each snow layer is governed by snowfall/rainfall, compaction, snowmelt, runoff, infiltration into the underly-ing snow layer or soil, and evaporation/sublimation at the snow surface.
The basic model process (Wang et al., 2009a) begins with delineation of the basin and sub-basins using the Pfafstetter scheme (Verdin and Verdin, 1999), division of sub-basins into a number of flow intervals based on time lag, prescription of all external parameters (e.g. land use, soil type, hillslope properties, and vegetation parameters), and meteorological forcing including precipitation on each model grid. Water, energy, and CO 2 fluxes are computed on each grid. Each grid maintains its own prognostic snow properties (snow water equivalent, snow depth, snow temperature, snow density, and ice/water content), and/or land surface temperature and soil moisture contents. Then, a grid-hillslope scheme generates slope-driven runoff, which is routed through the river network using the kinematic wave method. Overall model structure is illustrated in Fig. 1. Details of snow processes are given by Shrestha et al. (2010Shrestha et al. ( , 2012a, and other model processes are shown by Wang et al. (2009a, b).
It has been demonstrated that the WEB-DHM-S model is capable for accurate simulation of prognostic variables such as snow depth, snow water equivalent, snow density, snow surface temperature and snowmelt runoff. This was accomplished through rigorous evaluation of the model with comprehensive point snow measurements at Snow Model Intercomparison Project (SnowMIP) sites (Shrestha et al., 2010), the Valdai grassland (Shrestha et al., 2011), and Fraser Experimental Forest (Shrestha et al., 2012b). Moreover, a basinscale evaluation of the model in the Dudhkoshi region of the Nepal Himalaya showed that its simulated spatial distribution of snow cover agreed with MODIS snow cover data to an accuracy of 90 % (Shrestha et al., 2012a). This demonstrates the model capability for capturing spatiotemporal variations in snow cover across the study area.

Study area and data
The Yagisawa River basin (167 km 2 ) lies in a high, steep mountainous region of the Upper Tone Basin, northeast of Tokyo, Japan (Fig. 2a). The Yagisawa dam inlet is considered the basin outlet. Yagisawa dam is an important regulator of snowmelt runoff in spring and flooding in summer. The basin typically supplies 14.27 % of the water to Tone-Ara system to feed water supply to Tokyo metropolitan area. Furthermore, 14.11 m 3 s −1 and 2.918 m 3 s −1 discharge is used for irrigation water and city water respectively in Gunma prefecture. The climate in this region is wet and humid. February is the coldest month, with mean temperature −6 • C, and August is the warmest month, with temperatures averaging 18 • C. Heavy snowfall is common in winter (December through February), owing to a northwest monsoon wind from the Sea of Japan. The snowmelt period is from March through June. Heavy rainfall events in summer (July through October) are commonly associated with typhoons and Mei-yu frontal activity, which produce high flood risks  Wang et al., 2009a); (d) detailed description of vertical 3-layer energy balance snow model, for which T is temperature, e(T ) is vapour pressure at T , R sw and R lw are downward shortwave and longwave radiation, H and λE are sensible and latent heats, and ε, δ and α are emissivity, transmittance and reflectance, respectively. Subscript c refers to canopy, g to soil surface, sn to snow surface, and m to the reference height. Calculation of energy fluxes and aerodynamic resistances (r a , r b , r c , r d and r soil ) follows SiB2 (Sellers et al., 1996). D c and D t are canopy drainage and direct throughfall. H (Z j ) is enthalpy of the snow layer (Z j ) (after Shrestha et al., 2012b). in lower elevation regions. The dataset used is described in the following sections and is summarized in Table 1.

Meteorological data
The atmospheric forcing data necessary to drive the model include air temperature, air pressure, relative humidity, wind speed, and downward shortwave and longwave radiation at hourly time steps. Observed air temperature, wind speed, humidity, pressure and sunshine duration data were available at hourly resolution from meteorological sites of the Automated Meteorological Data Acquisition System (AMeDAS), provided by the Japan Meteorological Agency (JMA; Fig. 2a). Sunshine duration, wind speed, pressure and humidity data were interpolated to each model grid (500 m × 500 m) through the angular distance weighting (ADW) interpolation method (New et al., 2000). Air temperature was interpolated using the detrended ADW. First, air temperature at all AMeDAS sites was converted to a zero-elevation temperature, using a constant lapse rate of 6.5 • C km −1 . Second, the ADW was applied to the detrended data. Third, after the data were interpolated to each model grid, the lapse rate trend was added to each grid, based on its elevation. Then, the sunshine duration, temperature and humidity were used to calculate the downward shortwave radiation at each grid, using a hybrid model developed by Yang et al. (2001). Longwave radiation at each grid was then estimated from temperature, relative humidity, pressure and shortwave radiation,  using a relationship between shortwave and longwave radiation (Crawford and Duchon, 1999).

DEM, land use and soil data
Digital elevation model (DEM) data at 50 m resolution (Fig. 2b) and land-use data at 100 m resolution were obtained from the Japan Geographical Survey Institute. Basin elevation ranges from about 740 m to 2140 m, with mean 1285 m. Grid slopes vary from 3 • to 37 • , with mean 24 • . The land-use data were reclassified according to SiB2 categories (Sellers et al., 1996) with broadleaf deciduous trees being the dominant type (about 96 % of the basin; Fig. 2c). Static vegetation parameters, including morphological, optical and physiological properties for those SiB2 categorized vegetations, were defined following Sellers et al. (1996). We used dynamic vegetation parameters such as leaf area index (LAI) and fraction of photosynthetically active radiation (FPAR) absorbed by the green vegetation canopy as model inputs. These data are 8 day composites of MOD15A2 version 5.0 products, obtained from the MODIS Terra satellite at 1 km spatial resolution. The soil map was processed from a 1 : 200 000 scale Gunma Prefecture geologic map. The dominant soil type is forest soil, which covers about 60 % of the basin. Black, high-permeability, and red soils cover approximately 15, 15, and 10 %, respectively. Soil static parameters include saturated soil moisture content, residual soil moisture content, saturated hydraulic conductivity for soil surface and groundwater, and van Genuchten parameters (α and n; van Genuchten, 1980). Values of hydraulic conductivities were given by Yang et al. (2004), and other soil parameters were obtained from Food Agriculture Organisation (FAO, 2003).

Precipitation and discharge data
Precipitation gauge data were from meteorological sites of the AMeDAS. AMeDAS gauges in this region are the overflow tipping bucket type (JMA-RT4), with a heated reservoir and wind shield. Observed hourly discharge data at the dam inlet were obtained from the Ministry of Land, Infrastructure, Transport and Tourism, Japan.

MODIS snow cover data
The remotely sensed snow cover data used are from the 8 day maximum snow extent dataset (MOD10A2) of MODIS, aboard the Terra satellite. Tile h29V05 covered the entire study area. MOD10A2 depicts the maximum snow cover extent over an 8 day period at 500 m resolution, which is derived from daily snow cover products (MOD10A1) over such periods. MODIS 8 day product represents the maximum extent of snow cover over eight days, which effectively provides a temporal filter of MODIS daily data minimising cloud coverage. For forested regions, the normalized difference vegetation index (NDVI) and normalized difference snow index (NDSI) were jointly used to discriminate snow-free and snow covered forests, using the algorithm of Klein et al. (1998). All MODIS datasets were acquired from the NASA Earth Observing System Data and Information System (EOSDIS) (http://reverb.echo.nasa.gov), and processed by the MODIS Reprojection Tool (MRT, 2011).

Methodology for snowfall correction
The overall methodological framework is depicted in Fig. 3. All time-variant input data (downward longwave and shortwave radiation, wind speed, air temperature, relative humidity, LAI and FPAR at 1 h temporal resolution) and static data (DEM, soil, land use) were prepared on a 500 m grid. First, interpolated gridded air temperatures (T grid ) were used to define precipitation phase. We used 0 • C as the static threshold air temperature (T th ), below which all precipitation was assumed to be snowfall and above which it was assumed to be rainfall. The gauge precipitation was then interpolated to the grid using the ADW (New et al., 2000;Hofstra and New, 2009) interpolation method, taking into account the elevation dependent correction factor. Precipitation at each grid was estimated as the weighted sum of neighbouring precipitation gauges, with the weighted elevation-dependent correction factor considering the elevations of all neighbouring gauges used in ADW interpolation. This is described as follows.
where P grid (z) is corrected precipitation (m) at elevation z (m), P gauge (z i ) is observed precipitation for gauge i at elevation z i (m), W i is the angular distance weight factor for gauge i, ng is the total number of nearest-neighbour gauges contributing to the grid point during interpolation, and C fsnow /C frain is a calibration parameter (m −1 ) for orographic correction of snowfall/rainfall. Following New et al. (2000), for a grid point value, the weight factor (W i ) for gauge i, out of a total of ng contributing gauge stations, is where the position of the ith gauge is defined in terms of its distance x i and its angle to north θ i , relative to the specified grid point; w i is the distance weight. The correction factor is defined according to the precipitation phase (C frain for rainfall and C fsnow for snowfall). Earlier studies in the region show that rainfall-induced discharge matched well the observed flood peaks (Yang et al., 2004;Wang et al., 2009b). The correction factor for rainfall is not sensitive in the region, and thus we assumed C frain to be zero. However, C frain can be calibrated for basins where rainfall distribution is sensitive in the snow season.
For model execution, initial conditions for soil moisture and ground water storage were attained by running the model 100 number of time until hydrologic equilibrium was reached. The equilibrium condition is defined by setting the value of relative change (between two runs) in soil moisture and ground water storage to 0.1 %. Values for model parameters were taken from previous studies in the Upper Tone region, as indicated in Table 2. After simulation of model processes with initial arbitrary value of C fsnow , error between observed and simulated discharge at the basin outlet and error between MODIS-derived and simulated snow cover pixels were computed. The C fsnow was then optimised using the SCE-UA automatic search algorithm, which minimises the multi-criterion objective function (Z Err ), including the weighted components of discharge error (Q Err ) and snow pixel error (S Err ). The objective function is expressed by where α is the weight. The value of α should be given in such a way that both error components give the equal weights to the total error so that none of the error components overrule each other. The α is defined by For computation of discharge error, Nash-Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) and relative volume error (RVE) were used. Discharge error (Q Err ) is expressed by where Q oi is observed discharge at hour i, Q si is simulated discharge, and Q o is the average of observed discharge over the simulation (calibration/verification) period of comparison (N hours). NSE determines the relative magnitude of residual variance relative to measured data variance. RVE measures the average tendency of simulated values to be larger or smaller than observed ones. Positive values indicate overestimation bias and negative values indicate underestimation bias. NSE equal to 1 and RVE equal to zero correspond to perfect matching between simulated and observed discharge. S Err is expressed by combining the model overestimation error (M OE_AVG ) and model underestimation error (M UE_AVG ), averaged over the period of comparison (N d days).
These errors are calculated by performing a grid-to-grid analysis, to examine whether the MODIS and model simulation agree on whether the grid point is covered by snow. The model simulates the amount of snow (in snow depth or water equivalent) on each grid. The grid is considered snowcovered for snow depth greater than a threshold value. The MODIS snow cover data show only if the grid is covered by snow or land, or is classified as missing information (mostly caused by clouds). A 2 × 2 confusion matrix (

Results and discussion
The 4 yr hydrometeorological data were prepared for model simulation for October 2000-September 2004. As discussed earlier, error for the multi-criteria objective function (combination of discharge and snow pixel errors) was minimised using the SCE-UA search algorithm, for which an optimised value of C fsnow was obtained during a 1 yr calibration period (2002)(2003). as the calibration year because of the longer snow cover period and higher snowfall amount. The optimised C fsnow was 0.0007 m −1 , which indicated that the average snowfall would likely be 1.7 times that for an elevation difference of 1 km. The optimised or calibrated C fsnow was used to perform altitude correction of snowfall for 3 yr (2000-2001, 2001-2002 and 2003-2004), for validation of the modelling approach. Model performance indicators for discharge (NSE and RVE) and snow cover (M OE_AVG and M UE_AVG ) were determined for the validation period. These results were presented along with those for simulations with C fsnow = 0, to ascertain improvements in simulation of discharge and SCA. Further, model sensitivity owing to uncertainty in other model input and model parameters were analysed.

Discharge simulation
In the study region, snowfall normally accumulates from mid November through mid-March, and snowpack begins melting at the end of March. This produces remarkable diurnal variation of snowmelt discharge in the river. Snowmeltdriven discharge continues until mid June. Simulated vs. observed hourly discharges at the basin outlet are shown in Fig. 4a and b for March through June in calibration year 2003. These were for zero (no snowfall correction) and optimised values of C fsnow . A slow increase in discharge was seen for 1-11 April, after which there was a strong and continuous increase of observed discharge. Cold days on 23-26 April led to small discharge rates. Since snowpack amount had declined substantially by the beginning of June, the influence of strong melt events decreased. Correspondingly, the daily mean discharge rate showed less variability and diurnal fluctuations decreased in amplitude. Simulated hydrographs matched well the diurnal variation of observed discharge from the beginning of melt season (late March) to the beginning of May, for C fsnow = 0. On some dates, the model was unable to capture the diurnal variation of discharge well (e.g. 20-22 April), which may be due to less available energy estimation in the model. In May, the discharge was largely underestimated, because of less snowfall in winter. The volume deviation for discharge underestimation was 39 % (RVE = −0.39), with poor NSE (−0.03). After optimisation for snowfall correction, discharge simulation in late melt season (beginning of May to  (2003) and validation (2001, 2002 and 2004) periods, for optimised C fsnow and C fsnow = 0. mid June) greatly improved. Overall, the observed discharge curve was reasonably reproduced by the model, with NSE of 0.66 and RVE of −0.06 (Table 4).
For validation, the model was run using optimised C fsnow for 2001, 2002 and 2004. Observed and simulated hourly discharges for those years are shown for this optimised (Fig. 4d, f and h) and zero C fsnow (Fig. 4c, e and g). For zero C fsnow , the volume deviation was about 50 %, with poor NSE for all 3 yr. For the optimised value, the diurnal variation of dis-  Table 4). The simulated hydrograph for 2001 showed remarkable underestimation of discharge (about 15 %). The large discrepancies were observed from 3rd week of May to mid June, mainly due to underestimation of snowfall in winter. This implies that the lack of snowfall input was not much compensated by the optimised C fsnow . However, the Hydrol. Earth Syst. Sci., 18, 747-761, 2014 www.hydrol-earth-syst-sci.net/18/747/2014/  simulation of different discharge regimes in the 4 yr (onset of strong snowmelt discharge) was well replicated by the model. To illustrate the scatter of observed and simulated discharges, scatterplots of simulated versus observed discharge in the 4 yr are shown in Fig. 5a and b, for simulations with and without snowfall correction. These plots indicate that the discharge underestimation was greatly improved, and the correlation coefficient increased from 0.32 to 0.85 after snowfall correction.

Snow covered area simulation
Maps of the spatial distribution of snow cover derived from MODIS and that model at selected dates are presented in Fig. 6 (a: calibration period; b-d: validation period), for simulations with and without snowfall correction. These represent snow accumulation and snowmelt phases of the snowpack. A model-output pixel was considered snow-covered when snow depth exceeded or equaled 4 cm (e.g. Klein and Barnett, 2003;Wang et al., 2008). Model output was represented as the maximum snow extent over 8 day periods, corresponding to the dates of the MODIS dataset. The method used here did not assimilate snow covered area pixel by pixel; instead, it minimised the error between MODIS and model simulations during multi-objective optimisation of C fsnow . Thus, the model might not necessarily reproduce the same SCA as MODIS after correction.
During the calibration period (2002)(2003), simulated SCAs were underestimated at high elevations during early accumulation and late melt seasons (after 1 May), for model simulation with zero C fsnow (Fig. 6a). For optimised C fsnow , SCA was overestimated in late melt season. Table 5 summarizes a pixel-by-pixel comparison of simulated and MODIS SCAs, showing the ratio of pixels at which the simulated and MODIS snow cover agreed and disagreed (underprediction is when the model pixel was snow-free but MODIS was snow covered, and overprediction was vice versa), averaged over the comparison period. Average model underestimation error (M UE_AVG ) improved from 0.15 for zero C fsnow to 0.04 for optimised C fsnow . Average M OE_AVG showed opposite characteristics with regard to M UE_AVG , as expected. In the validation period after snowfall correction, M OE_AVG was 0.098, 0.135 and 0.10 in 2004and 2001. M UE_AVG was from 0.014 to 0.021. Snow pixel underestimation error on 17 May was the highest (0.12) in 2001 among all years which supported the clarification of large RVE in discharge simulation due to underestimation of winter snowfall. For simulations without correction, overestimation error was nearly zero because of early snowmelt, with underestimation error from 0. 06 (2004 and 2002) to 0.10 (2001).
It is interesting that the model generally overestimated snow pixels for corrected snowfall; i.e. total S Err increased with respect to error for C fsnow = 0. The optimisation of C fsnow against SCA could only reduce the error, but it would definitely produce large error in the discharge components. Similarly, the optimisation against discharge would only reduce discharge error, but it would largely overestimate SCA in melt season. Thus, the multi-criteria objective function was needed to reach a tradeoff between the two error components. Although discharge was well simulated by use of the optimal correction factor, there was considerable deterioration of model performance regarding overestimation of snow pixels during melt season. There may be two basic contributions to this error. The first is uncertainty of the MODIS dataset algorithm in mapping snow cover in forested regions during melt season. Hall et al. (2001) stated that snow mapping performance error of the MODIS daily product for such regions was the largest (≈ 15 %) of all land cover types.  (2004) reported that MODIS mapping accuracies were poorest in evergreen forests, with an error rate of 20 % during snowmelt. In Austria, Parajka and Blöschl (2006) reported a mean misclassification error for the shrub class around 10 %, and for pastures and forest around 6 %. Parajka et al. (2012) addressed MODIS snow-cover mapping accuracy in a forested region, finding that most of the mapping errors were at the end of snowmelt season for patchy and shallow depths (mean snow depth typically less than 10-15 cm). The model outputs snow on the ground surface, whereas MODIS shows a satellite view of snow in the visible spectrum, as a mosaic of the vegetation (NDVI) and snow (NDSI) algorithms. The diurnal variation of observed discharge clearly showed that there must have been snow in high-elevation areas during late May, which MODIS was unable to capture. The second contribution to the error may be attributed to noise during interpolation of precipitation, for model grids distant from gauges and for very weak correlation of snowfall between gauges. A comprehensive study should be made in the future, to ascertain the reasons for misclassification of snow pixels by the current MODIS snow mapping algorithm. In general, the model was able to simulate seasonal and interannual variability of snow cover relative to MODIS, including light snow accumulation in mid-November, full snow coverage from December through April, and persistence of snow at high elevations in May. Figure 7 portrays the spatial distribution of total snowfall before and after correction, for the 4 yr of the study period. Total snowfall was calculated as the sum of snow accumulation from October through May. The results for zero C fsnow clearly show a region of snowfall greater than 800 mm in the northwest part of the basin. The year 2003 had heavy snowfall, so the 800 mm contour moved toward lower elevations. Basin-average total snowfall was 756, 786, 850 and 727 mm over -2004). Monthly analysis of total snowfall illustrated that the basin had three-peak total snowfall in 2003 (November, January and March with the highest in January) and one-peak snowfall in 2001 (February), 2002 (February) and 2004 (January).

Spatial distribution of snowfall
For optimised C fsnow , the corrected total snowfall varied from ≈ 800 mm at low elevations (near the basin outlet) to 2000 mm at high elevations of the upper region of basin. After correction, the 800 mm contour for uncorrected snowfall became about 1200 mm. Basin-average corrected total snowfalls for 2001-2004 were 1124, 1169, 1263 and 1077 mm. These values are about 49 % higher than uncorrected ones (Table 6). Through our efficient calibration and validation approaches, the spatial distributions of snowfall estimated herein can be taken as reference in hydrologic modelling and in the correction of radar data and reanalysis datasets in the region.

Sensitivity to uncertainties in model input
Despite the precipitation correction, there may be uncertainties in other model input, which arise from observations, observation-based meteorological models, and remotesensing products. To obtain insight regarding the response of discharge and SCA simulations to these uncertainties, we performed sensitivity analyses of other inputs. This analysis was done for 2003, which was referred to as the control run. Several simulations were made for a prescribed range of variations (±10 %) of model input (shortwave radiation, wind speed, air temperature, and LAI/FPAR). Longwave radiation is a function of air temperature, so its variation was  (2003) and validation (2004, 2002 and 2001) periods, for simulations with and without correction for snowfall. not considered. The discharge and snow pixel error components were analysed for each model run.
It was speculated that decreased shortwave radiation increased the overestimation error of snow pixels to 0.1853, and reduced underestimation error to 0.0108, giving an NSE of 0.53 and RVE of −0.03. For increased shortwave radiation, discharge was slightly high in early May and was slightly underestimated in mid June. M OE_AVG and M UE_AVG were 0.07 and 0.05, with NSE of 0.57 and RVE of −0.11 (Table 7). With increase of air temperature, the ratio of snowfall  to total precipitation normally decreases, but if that temperature is below freezing during the entire snowfall period, this ratio will have much less sensitivity. Since incoming longwave radiation was treated as a function of air temperature, the higher that temperature, the greater the longwave radiation and vice versa. However, variation of air temperature results in very slight change in longwave radiation. For instance, a 10 % change of temperature in the study area would cause 0.5 to 2 % variation in longwave radiation, depending on season. A small increase in discharge in early April was simulated for a 10 % increase in temperature, causing NSE to decline to 0.6315 and an RVE of −0.10. M OE_AVG decreased to 0.0966, but M UE_AVG increased to 0.0357. The decrease in air temperature caused a late snowmelt, owing to a small decrease of discharge in early April. This caused an increase of M OE_AVG and decrease of M UE_AVG (Table 7). Variation in dynamic vegetation parameters such as LAI may affect the water/snow holding capacity of the canopy at point scale. Since WEB-DHM-S inherited SiB2 (Sellers et al.,  Table 7. Change in wind speed affects the computation of turbulent fluxes. As expected, this change did not greatly affect discharge, and thus showed the least sensitivity. There was greater sensitivity to a ±10 % change in shortwave radiation among all inputs, since it directly controls available energy.

Sensitivity to uncertainties in model parameters
A large number of uncertainties may exist in model parameters (e.g. snow albedo, roughness length for snow surface, threshold temperature for snow/rain, and morphological parameters of the canopy), initial conditions, and soil properties (Xue et al., 1997). Among these parameters, threshold temperature is critically important for discharge generation, since it determines the proportions of rainfall and snowfall in total precipitation. Its effect is dominant in lower and more humid elevations, but diminishes gradually with increasing elevation. Simulations were carried out for threshold temperatures (T th ) of 0.5, 1.0, 1.5 and 2 • C. With increasing T th , NSE and RVE decreased from 0.66 to 0.53 and −0.06 to −0.03 (see Table 8). Model overestimation increased from 0.1187 to 0.1818, whereas underestimation error showed less variability. It is evident that with the increase of T th , C fsnow decreases and overestimation error can be reduced. We executed the SCE-UA search algorithm to optimise C fsnow for T th of 2 • C. It was found that the minimum objective function value (Z Err ) was higher than that of the control run. From the sensitivity runs, T th was a key parameter for simulating discharge and SCA; however, the default parameter (T th of 0 • C) showed the best result. We also tested the sensitivity of model simulation to a change of fresh snow albedo, which influences net radiation flux. The effect of this albedo was examined by changing the  (Table 8). For decreased albedo, discharge error slightly worsened, but snow pixel error improved. Nonetheless, in both cases, total error exceeded that of the control run. Overall, these results show that the optimised C fsnow of the control run was the best value for minimisation of errors in discharge and snow cover.

Conclusions
This study focused on correction and estimation of the spatial distribution of snowfall. We used the spatially distributed multilayer water and energy balance-based snowmelt model (WEB-DHM-S) and remotely sensed snow cover data from MODIS aboard the Terra satellite, following multi-objective optimisation of altitude-based snowfall correction factor C fsnow within the framework of the SCE-UA automatic search algorithm. The minimum objective function was achieved by minimising the difference between observed and simulated discharge at the basin outlet and the difference between MODIS-derived and simulated snow-covered pixels.
The system was applied to Yagisawa Basin of the Upper Tone River in Japan.
The optimisation of C fsnow was carried out for the snow season 2002-2003. For the minimum errors in discharge and SCA simulation, optimum C fsnow was 0.0007 m −1 . To assess improvement in model simulation over that with no snowfall correction, the model was also run for C fsnow = 0. NSE improved from −0.03 to 0.66, and RVE from −0.39 to −0.06. The snow covered pixel underestimation error greatly decreased, but overestimation error increased for C fsnow = 0.0007. The optimised C fsnow was validated by model simulations for 2004, 2002 and 2001, for which discharge and SCA were well simulated. The spatial distribution of total snowfall was estimated in each year, and varied from 1077 mm to 1260 mm. The 4 yr average of basin-scale total snowfall for optimised C fsnow was 1160 mm, which was about 1.49 times that for zero C fsnow .
Hydrol. Earth Syst. Sci., 18, 747-761, 2014 www.hydrol-earth-syst-sci.net/18/747/2014/ Although C fsnow was optimised for the least cumulative error for discharge and snow-covered pixels, snow pixels were generally overestimated by the former in late melt season (May). This overestimation may be attributable to a known weakness of MODIS in mapping snow over forested regions. It may be argued that forest is one of the most complex land uses for the snow-mapping algorithm. Previous studies have shown that MODIS snow-mapping performance error for forest areas was the greatest (around 10-20 %) of all land cover types, especially when the snowpack was shallow and patchy. The diurnal variation of observed discharge clearly showed that there must have been snow in high elevation areas during late May, which MODIS was unable to capture. However, detailed analysis should be done in the future to determine the reasons for the misclassification of snow pixels by the current MODIS snow-mapping algorithm.
Model sensitivity runs were executed to observe the impact of uncertainty in other model input and model parameters. It revealed the model's high sensitivity against the shortwave radiation input data (relative to air temperature, wind speed and LAI). In addition, the uncertainty in threshold air temperature parameter showed the considerable variation in model simulations. An attention should be given in the preparation of this sensitive input and in the selection of T th value. However, the default model input and parameters yielded the best value for the optimised C fsnow which demonstrated the robustness of our results.
The method used here is simple and robust, and can be applied to any snow-fed river basin to obtain a reliable estimate of snowfall correction factor via the energy balance-based distributed snow model. Research using such a hydrologic modelling approach to different basins is needed, to further develop or validate the approach for climate and land-surface hydrologic research on snow-dominated, mountainous river basins. A future step would be application of the approach to correction of snowfall amount in reanalysis products and atmospheric model outputs, which would contribute to bias correction of climate model projections. Finally, a satellite snow-cover data assimilation system could be a future prospect for optimisation of the snowfall correction factor.