Spatial snowfall distribution in mountainous areas estimated with a snow model and satellite remote sensing

: We describe an approach for reconstructing snowfall that combines satellite observations of the snow disappearance date (SDD) with a snow model for two mountainous areas in Japan having distinct snow climatology. This approach allows assessment of the distribution of snow and topographical effects on snowfall within a catchment. We also evaluated how the reconstructed snowfall affects the catchment snow hydrology. Validation at observation sites demonstrated that a combination of the snow model and snowfall reconstructions successfully estimated the seasonal changes of snow water equivalent (SWE). In Japan, the dependence of snowfall on elevation is stronger in mountainous areas along the coast on the Sea of Japan side of Japan’s central mountain spine (where snowfall triples with every 1000 m increase in elevation; C sf = 0.002 m − 1 ) compared with inland locations of the same region (where snowfall doubles with every 1000 m increase; C sf = 0.001 m − 1 ). Moreover, the reconstructed snowfall improved the estimation of maximum catchment SWE. Maximum total SWE, estimated with reconstructed snowfall, was 3.8 × 10 8 m 3 in Kurobe catchment (along the coast on the Sea of Japan side of Japan’s central mountain spine), while that, estimated by convectional method with the spatially-constant C sf = 0.001 m − 1 , was 2.0 × 10 8 m 3 . As a result, estimations of the snow disappearance date and of the catchment snowmelt were also improved. These results suggest that it is useful to estimate the spatial snow distribution, especially where steep topography causes large gradients of snowfall amount with respect to elevation.


INTRODUCTION
Snowfall is a key component of the hydrological cycle in snowy climates. As snowfall is the first input in a hydrological analysis in such regions, spatially accurate snowfall data are needed for both water balance analysis and water resources assessment with hydrological models (Tian et al., 2007;Asaoka et al., 2010). In Japan, Motoya et al. (2001) and Shinohara et al. (2009) evaluated water balance over heavy snow catchments with hydrological models; these analyses were mostly driven by spatial datasets interpolated from routine meteorological station data. A serious problem with interpolated meteorological data is the scarcity of rain gauges and snow survey data in mountainous areas, where the complex topography introduces large spatial variability in the snowfall and snow water equivalent (SWE) (Asaoka et al., 2002).
Several factors can influence the estimation of spatial snowfall: whether precipitation is classified as snow, undercatchment of snowfall in gauges due to the wind effect, and spatially complicated topography. Yamazaki (2001) suggested a formula for the ratio of snowfall to total precipitation as a function of wet-bulb temperature, as analyzed by Sugaya (2005). Yokoyama et al. (2003) proposed a formula for the gauge undercatchment as a function of wind speed. Hirabayashi et al. (2008) constructed a global snowfall dataset by separating precipitation types and correcting for gauge undercatchment. To handle the heterogeneous distribution of snowfall, researchers have often used a simple function of elevation for interpolating snowfall from routine meteorological observations, under the assumption that the relationship of snowfall to elevation is constant over a large area.
Satellite-derived snow cover area (SCA) is effective for evaluating the uncertainty in the spatial snow distribution. Shimamura et al. (2003) estimated SWE in a catchment by combining satellite-derived SCA with in-situ data. Molotch et al. (2004) and Kazama et al. (2008) successfully combined the satellite-derived SCA data with distributed snowmelt models to estimate the spatial SWE during the snowmelt season. In most research, an uncertain parameter is assumed to be constant over the catchment or other large area. In reality, the complex terrestrial topography of Japan greatly affects the spatial snowfall or snowmelt. Thus, our research attempts to obtain the grid-based topographic parameter, the elevation dependence of snowfall, for reconstruction of the snowfall distribution over mountainous areas.
Our purpose is to reconstruct the snowfall distribution in Japanese catchments in mountainous terrain by combining satellite-derived SCA data with a snow model. We used SPOT/VEGETATION (SPT/VGT) data for determining the snow disappearance date (SDD). The elevation dependence of snowfall is analyzed in two geographically distinct mountainous regimes in Japan: 1) along the coast on the Sea of Japan side of Japan's central mountain spine; and 2) within more inland areas on the same side of the central mountain spine. We compared our reconstruction approach with the conventional methods for snowfall estimation that assume that the elevation dependence of snowfall is spatially constant when interpolating meteorological station data. Finally, we evaluated the effects of our reconstructed snowfall on catchment snow hydrology.

DATA AND METHODS
The snowfall reconstruction approach was used in the two mountainous areas (Supplement Figure S1) which are represented by the Kurobe Dam catchment (185 km 2 ) and the Tadami Dam catchment (856 km 2 ). Kurobe Dam (36°33'59''N 137°39'44''E) is located close to the coast on the Sea of Japan side of Japan's central mountain spine and its catchment ranges in elevation from about 1500 m at the dam site to 3000 m at the highest peak. Tadami Dam (37°20'9''N 139°18'1''E) is located more inland than Kurobe Dam on the same side of the central mountain spine; the elevation of its catchment ranges from 400 to 2300 m. We used meteorological data from the Japan Meteorological Agency's AMeDAS (Automated Meteorological Data Acquisition System) stations, which automatically record temperature, precipitation, sunshine duration and wind velocity, as well as snow depth at some sites. Moreover, we used observations of SWE (Yamaguchi et al., 2007) automatically recorded at four sites (St. A, B, C and D in Supplement 1) outside of the catchments for validation of the reconstruction approach.
The reconstruction approach estimates cumulative snowmelt during the season of snow cover on the assumption that it equals total snowfall. The duration of snow cover in each grid was estimated from information on SCA derived from the SPT/VGT. Next, total snowfall was apportioned to daily snowfall values with a snow model. SPT/VGT images are composites of 10 days of Normalized Difference Vegetation Index data: for each pixel of these images, a maximum value is found from 10 days of daily SPT/VGT images (Holben, 1986;Duchemin et al., 2002). These data have a spatial resolution of approximately 1 km at satellite nadir. Snow cover detection is accomplished using the S3 index (Saito and Yamazaki, 1999), which is based on three bands of electromagnetic radiation: red (0.61-0.68 μm), near-infrared (0.78-0.89 μm) and short-wave infrared (1.58-1.75 μm). The pixels with S3 greater than 0.05 are mapped as snow covered (Shimamura et al., 2006). We performed the detection for 21 SPT/VGT scenes from January to July 2000, and extracted a snow disappearance date (SDD) in each grid of satellite images at 10 day intervals (Figures 1a  and 1c). Some pixels may be covered by cloud even if SPT/ VGT images are 10-day composites. SPT/VGT images have information about cloud cover. Therefore, cloud-cover pixels are excluded from detecting snow SDD. The extracted SDD in the SPT/GVT images were resampled to the grid of a digital elevation model with 1 km resolution using a nearestneighbor method. The snow model estimates the beginning of the season of snow cover, in a manner described later in this chapter. Snowmelt rate was estimated with a degreeday method (Rango and Martinec, 1995) as follows: where M is daily snowmelt (mm), k is a degree-day factor (mm day −1 °C −1 ) and T a is daily mean air temperature (°C). A base temperature below which snowmelt does not generate, is defined as 0°C. The degree-day factors were optimized from meteorological data at the AMeDAS station at Hakuba for the Kurobe area, and the station at Tadami for the Tadami area, by assimilating cumulative snowmelt with cumulative snowfall during the snow cover season. We obtained snowmelt rate factors of approximately 7.0 mm day −1 °C −1 at both sites; they were assumed to be spatially constant in this study.
Daily precipitation was interpolated from precipitation at AMeDAS sites over each grid as follows: where Pr is daily precipitation at the elevation z (m), C sf is the coefficient of the vertical gradient of snowfall with respect to height z, and elv is the elevation. The subscript "obs" refers to AMeDAS observation sites. This approach does not need to use the coefficient for rainfall because the result of snowfall reconstruction does not depend on rainfall.
The ratio s of snowfall to total precipitation was obtained from the wet-bulb temperature T w as follows (Yamazaki, 2001):  (2) T w is obtained from the daily temperature and water vapor pressure. Daily temperature at any point was calculated with a lapse rate of 6.5°C/km and we used daily temperature data at the AMeDAS sites. Water vapor pressure is interpolated from the nearest meteorological observatory or weather station of the Japan Meteorological Agency. We assumed that the gauges catch 0.8 of the snowfall because of windinduced undercatchment of snow (Kominami et al., 2005), and all (1.0) of the rainfall. The coefficient C sf was obtained over each grid with the assumption that cumulative snowfall equals cumulative snowmelt during the snow cover season as follows: ( 3) where s i is the ratio of snowfall to total precipitation and i is time interval and i = 1 represents the first date of snow cover. Pr was calculated with the three nearest AMeDAS satations by equation (1) and the inverse distance weighted interpolation, and we assumed that C sf is constant for the three stations.
Finally, in order to validate this approach and to evaluate the impact of reconstructed snowfall on catchment hydrology, snow water equivalent was estimated as follows: where SF is daily snowfall and can be defined as sPr. The time interval is one day and we assumed that all grids had no snow cover at the beginning of model calculations.

Validation of the snowfall reconstruction approach with SWE
The reconstructed snowfall was validated by comparing estimated SWE with observed SWE data at four sites ( Figure  2). The degree-day factor for each site was optimized at the closest AMeDAS station, and two types of snowfall data were used for SWE estimation: 1) reconstructed snowfall (thin line in Figure 2) and 2) precipitation data at the closest AMeDAS station (i.e. C sf = 0 m −1 ), after separation of precipitation type at the validation points and corrections for undercatchment (dashed line in Figure 2). The snow model, when coupled to the snowfall reconstruction approach, successfully estimated the seasonal change of SWE and the duration of the snow cover season. There were short-term errors near the peak values of SWE; the relative errors at the peak season were 23% at St.D and 10-15% at the other sites. The snowfall forcing data were based on interpolation from AMeDAS sites in lowland areas. Thus, interpolated snowfall was underestimated wherever snowfall occurred in the highlands but did not occur in the lowlands. The C sf values obtained by the reconstruction approach were 0.001, 0.001, 0.002, and 0.0015 m −1 at site A, B, C, and D, respectively. These results lead us to suspect optimization of C sf over each grid and constructed snowfall input into a snow model would increase accuracy in SWE and snowmelt estimation. The other estimation approach, which used the snow model and AMeDAS precipitation data, underestimated SWE throughout the snow cover season. As a result, less snowfall was input into the model, which caused estimates of SDD to be too early, especially at St.C (by 25 days) and St.D (by 13 days). Spatial snowfall distribution during the snow cover season Supplement Figure S2 represents the cumulative snowfall during the snow cover season in two study areas. The spatial average of optimized C sf was approximately 0.002 m −1 for the Kurobe area, where snowfall accumulation triples with every 1000 m rise in elevation, and was 0.001 m −1 for the Tadami area, where snowfall doubles with every 1000 m rise. Also note that C sf = 0.002 and 0.001 m −1 for St.C and St.B, respectively, for the validation represented in the previous section. Kondo et al. (1995) had obtained C sf of approximately 0.001 m −1 from a water budget in a mountainous catchment on the Sea of Japan side of the central mountain spine of the northeast part of Japan. Therefore, Motoya et al. (2001), Asaoka et al. (2007) and Kazama et al. (2008) applied a C sf value of 0.001 m −1 to estimate SWE distribution in Japan, and proved that this value was useful for estimating the average distribution of SWE in the central mountain spine of Northern Japan. Our study also obtained a spatial average of 0.001 m −1 for the optimized C sf in the inland area west of the central mountain spine; i.e. the Tadami area and St.B (Okutadami-Maruyama). The consistency of our optimized C sf values with values from previous research suggests that our reconstruction approach is valid, but it does not prove that a C sf value of 0.001 m −1 is valid for smaller scales or for regions outside of the central mountain spine. In some mountainous areas adjacent to the Sea of Japan, annual precipitation ranges from 4000 to 12000 mm (Wada, 2004;Tsuchiya, 2005), which is much greater than the precipitation estimated with a C sf value of 0.001 m −1 .
The steep topography in the Kurobe area and around St.C (Myokoo-Sasagamine) leads to pronounced effects on precipitation when the wind impinges upon the steep mountain topography with saturated water vapor from the Sea of Japan. The Kurobe area is closer to the Sea of Japan and has steeper topographical gradients than does Tadami; the ground around Kurobe rises from 0 to 3000 m over a horizontal distance of 20 km. As a cold and dry winter monsoon blows from the Asian continent during the winter, heat and vapor are imparted from the relatively warm sea surface to the lowest atmospheric layer over the Sea of Japan. This causes the stability of the lower atmosphere to decrease, which leads to formation of cumulonimbus clouds, which when lifted over steep mountain slopes release huge amounts of precipitation. In Tadami, by contrast, lower elevaton at the mountain peaks and less steep mountain topography result in a more moderate dependence of snowfall rate on elevation. Another reason is that the winter wind, originally saturated with water vapor on arriving from the Sea of Japan, arrives at Tadami after having relinquished a large fraction of its water vapor as precipitation upon the intervening mountain ranges.
Snowfall at ground level depended strongly on the elevation of the ground (Figure 3) in the Kurobe Dam catchment. Both catchments received more snowfall from 800 to 1400 m in elevation than the regression line would indicate. Most snow-producing clouds, generated over the Sea of Japan and blown by the wind against the mountain slopes, either drift upward or diverge around the mountain range in a horizontal direction. The height of cloud base is around 1 km above sea level (Nakai and Endoh, 1995). Hence, there is a steep dependence of snowfall on elevation in this altitude range. Moreover, winds cause snow particles to drift on the lower mountain slopes. Optimization of the C sf coefficient in each grid would reflect these phenomena, except that this approach would not detect where the land surface is exposed on mountain peaks or ridges due to strong wind (Shimamura et al., 2005), as the spatial resolution of SPT/VGT is too coarse to resolve such small-scale phenomena. Also note that C sf indicates the relationship of snowfall between each grid and AMeDAS sites, therefore, the snowfall gradient with respect to elevation in Figure 3 is smaller than the spatial average of optimized C sf by the reconstruction approach.

Impact of snowfall on catchment snow hydrology
To evaluate the effect of the snowfall reconstruction approach, we estimated the seasonal changes of catchment SWE and snowmelt in the study areas by using two methods: 1) with the snow model and reconstructed snowfall data (using an optimized C sf in each grid) and 2) with the snow model and interpolated snowfall data, under the assumption that C sf was constant at 0.001 m −1 (Figure 4). The two methods yielded the same seasonal variation of SWE in the Tadami catchment, but a large difference in SWE in the Kurobe catchment: the maximum total SWE was 3.8 × 10 8 m 3 with the reconstruction approach but only 2.0 × 10 8 m 3 with the spatially constant C sf . The mean optimized value of C sf in the Kurobe area was approximately 0.002 m −1 , and spatially optimized values of C sf tended to increase with elevation. In the Tadami area, on the other hand, the mean optimized C sf was approximately 0.001 m −1 , and this value did not tend to increase with elevation, so that both The differences in the estimated SWE between the two approaches had a large impact on snowmelt after the time of the maximum SWE in the Kurobe catchment. Estimated daily snowmelt in the catchment was about 300% greater using the reconstruction method than it was using the spatially constant C sf , when averaged from DOY (day of the year) 140 (mid-May) to DOY 170 (mid-June). With remote sensing data, we detected that snow started to disappear after DOY 150 (end of May) in the year of this study (Figure 1a, while the method that uses a spatially constant C sf estimated that snow would disappear after DOY 130 (early May) (Figure 1b). Therefore, the two methods differed most of all in their estimates of snowmelt after DOY 130. Combination of a snow model with a satellitederived SDD can avoid an estimate for disappearance of snow that would be too early. Use of the reconstructed snowfall also resulted in successful estimates of the seasonal SWE in the catchment, and the amount of snowmelt after the time of peak SWE. In contrast, the two methods were in good agreement for the Tadami catchment because the SDD estimated by the snow model corresponded well with the satellite-derived SDD (Figures 1c and 1d).

CONCLUSION
Remotely-sensed snow information was combined with a snow model to reconstruct the snowfall in mountainous areas. This approach improved the seasonal estimation of SWE, especially its peak seasonal value and the snow disappearance date at observation sites in heavy snow areas. The conventional approach for estimating snowfall assumed the elevation dependency of snowfall was constant over a catchment or large area. Our approach allowed grid-based assessment of snowfall, and demonstrated that the elevation dependency of snowfall is distinct depending on the catchment or geographical condition. The computed spatial snowfall distribution was a steeper function of elevation in the mountainous area along the coast on the Sea of Japan side of the central mountain spine, where snowfall tripled with every kilometer increase in elevation (C sf = 0.002 m −1 ), while it doubled with every kilometer increase (C sf = 0.001 m −1 ) on the more inland side of the central mountain spine. Moreover, we found that using reconstructed snowfall in catchment snow hydrology prevented underestimation of the maximum SWE or snow disappearance date (SDD) that was too early. Especially, this approach estimated 3.8 × 10 8 m 3 maximum of total SWE in Kurobe catchment (along the coast, the Sea of Japan side of the central mountain spine), while the conventional approach, with the spatially-constant coefficient of the vertical gradient of snowfall with respect to height (C sf = 0.001 m −1 ), estimated 2.0 × 10 8 m 3 . The beginning of SDD by the conventional method was approximately 20 days earlier than that of satellite-derived SDD. It followed that this approach reduced the present large underestimation of catchment snowmelt and duration of snowmelt season. Such results confirm that this approach is a more effective way to estimate snow hydrology, particularly in the heavy snow areas affected by the steep topography and high moisture content in air arriving from the Sea of Japan: that is, along the coast on the Sea of Japan side of Japan's central mountain spine. This method of snowfall estimation also has the potential to provide more accurate and useful information to water resources managers in heavy snow areas. (GREF) of the Ministry of the Environment of Japan for financial support. We are also grateful to Dr. Satoru Yamaguchi of the National Research Institute of Earth Science and Disaster Prevention for providing snow observation data. Two anonymous reviewers and the associate editor are gratefully acknowledged for their insightful and constructive remarks.