Assimilation of snow cover and snow depth into a snow model to estimate snow water equivalent and snowmelt runoff in a Himalayan catchment

<jats:p>Abstract. Snow is an important component of water storage in the Himalayas. Previous snowmelt studies in the Himalayas have predominantly relied on remotely sensed snow cover. However, snow cover data provide no direct information on the actual amount of water stored in a snowpack, i.e., the snow water equivalent (SWE). Therefore, in this study remotely sensed snow cover was combined with in situ observations and a modified version of the seNorge snow model to estimate (climate sensitivity of) SWE and snowmelt runoff in the Langtang catchment in Nepal. Snow cover data from Landsat 8 and the MOD10A2 snow cover product were validated with in situ snow cover observations provided by surface temperature and snow depth measurements resulting in classification accuracies of 85.7 and 83.1 % respectively. Optimal model parameter values were obtained through data assimilation of MOD10A2 snow maps and snow depth measurements using an ensemble Kalman filter (EnKF). Independent validations of simulated snow depth and snow cover with observations show improvement after data assimilation compared to simulations without data assimilation. The approach of modeling snow depth in a Kalman filter framework allows for data-constrained estimation of snow depth rather than snow cover alone, and this has great potential for future studies in complex terrain, especially in the Himalayas. Climate sensitivity tests with the optimized snow model revealed that snowmelt runoff increases in winter and the early melt season (December to May) and decreases during the late melt season (June to September) as a result of the earlier onset of snowmelt due to increasing temperature. At high elevation a decrease in SWE due to higher air temperature is (partly) compensated by an increase in precipitation, which emphasizes the need for accurate predictions on the changes in the spatial distribution of precipitation along with changes in temperature.                     </jats:p>

section_content: In the Himalayas a part of the precipitation is stored as snow and ice at high elevations. This water storage is affected by climate change resulting in changes in river discharge in downstream areas (Barnett et al., 2005; Bookhagen and Burbank, 2010; Immerzeel et al., 2009 Immerzeel et al.,, 2010)). The Himalayas and adjacent Tibetan Plateau are important water towers, and water generated here supports the water demands of more than 1. 4 billion people through large rivers such as the Indus, Ganges, Brahmaputra, Yangtze and Yellow River (Immerzeel et al., 2010). So far, the main focus has been on the effect of climate change on the glaciers and the resulting runoff. However, snow is an important short-term water reservoir in Published by Copernicus Publications on behalf of the European Geosciences Union. 

the Himalayas, which is released seasonally, contributing to river discharge (Bookhagen and Burbank, 2010; Immerzeel et al., 2009). The contribution of snowmelt to total runoff is highest in the western part of the Himalayas and lowest in the eastern and central Himalayas (Bookhagen and Burbank, 2010; Lutz et al., 2014). 

Although Himalayan snow storage is important for the water supply in large parts of Asia, in situ observations of snow depth are sparse throughout the region. Many studies benefit from the continuous snow cover data retrieved from satellite imagery to estimate snow cover dynamics or contribution of snowmelt to river discharge (Bookhagen and Burbank, 2010; Gurung et al., 2011; Immerzeel et al., 2009; Maskey et al., 2011; Wulf et al., 2016). Studies about snowmelt in the Himalayas have predominantly relied on remotely sensed snow cover and a modeled melt flux estimating melt runoff resulting from this snow cover (e. g., Bookhagen and Burbank, 2010; Immerzeel et al., 2009; Tahir et al., 2011; Wulf et al., 2016). However, this approach provides no or limited information on snow water equivalent (SWE), which is an important hydrologic measure as it indicates the actual amount of water stored in a snowpack. SWE can be reconstructed based on the integration of a simulated melt flux over the time period of remotely sensed observed snow cover. However, this method provides only information on the peak SWE value and introduces errors when snowfall occurs during the melt season (Durand et al., 2008; Molotch, 2009; Molotch and Margulis, 2008). Currently there is only limited reliable information available on SWE for the Himalayas (Lutz et al., 2015; Putkonen, 2004). SWE can be retrieved with passive microwave remote sensing, but the results are highly uncertain, especially for mountainous terrain and wet snow (Dong et al., 2005). In addition, the spatial resolution is coarse and therefore inappropriate for catchment scale studies in the Himalayas. Estimating both the spatial and temporal distribution of SWE and snowmelt is important for flood forecasting, hydropower and irrigation in downstream areas. 

Selection of a suitable snow model is critical to correctly represent snow cover and SWE. Snow models of different complexity exist and can be roughly divided into physically based and temperature-index models. Several studies have compared snow models of different complexity and their performance. Physically based models typically outperform temperature-index models in snowpack runoff simulations on a sub-daily timescale (Avanzi et al., 2016; Magnusson et al., 2011; Warscher et al., 2013). However, physically based and temperature-index models have a similar ability to simulate daily snowpack runoff (Avanzi et al., 2016; Magnusson et al., 2015). Avanzi et al. (2016) showed that the use of a temperature-index model does not result in a significant loss of performance in the simulation of SWE and snow depth with respect to a physically based model. Even though physically based models outperform temperature-index models in some cases, temperature-index models are often preferred, as data requirements and computational demands are lower. 

Especially in the Himalaya, data availability constrains the choice of a snow model. 

Assimilation of remotely sensed snow cover and groundbased snow measurements has been proved to be an effective method to improve hydrological and snow model simulations (Andreadis and Lettenmaier, 2006; Clark et al., 2006; Leisenring and Moradkhani, 2011; Liu et al., 2013; Nagler et al., 2008; Saloranta, 2016). Although different data assimilation techniques exist, Kalman filter techniques are often selected, due to their relatively low computation demand. They estimate the most likely solution using an optimal combination of observations and model simulations. Especially in catchments with strong seasonal snow cover, assimilation of remotely sensed snow cover is expected to be most useful as a result of fast changing conditions in the melting season (Clark et al., 2006). 

The aim of this study is to estimate SWE and snowmelt runoff in a Himalayan catchment by assimilating remotely sensed snow cover and in situ snow depth observations into a modified version of the seNorge snow model (Saloranta, 2012 (Saloranta,, 2014 (Saloranta,, 2016)). Climate sensitivity tests are subsequently performed to investigate the change of SWE and snowmelt runoff as result of changing air temperature and precipitation. The approach of modeling snow depth allows us to validate the quantity of simulated snow rather than snow cover alone and is a new approach in Himalayan snow research.
section_name: Introduction

section_name: Methods and data

section_content: The study area is the Langtang catchment, which is located in the central Himalayas approximately 100 km north of Kathmandu (Fig. 1 ). The catchment has a surface area of approximately 580 km 2 from the outlet near Syabru Besi upwards. The elevation ranges from 1406 m above sea level (a. s. l. ) at the catchment outlet to 7234 m a. s. l. for Langtang Lirung, which is the highest peak in the catchment. The climate is monsoon dominated and 68-89 % of the annual precipitation falls during the monsoon (Immerzeel et al., 2014). Spatial patterns in precipitation are seasonally contrasting, and there is a strong interaction between the orography and precipitation patterns. At the synoptic scale, monsoon precipitation decreases from south to north, but at smaller scales local orographic effects associated with the aspect of the main valley ridges (Barros et al., 2004) determine the precipitation distribution. Numerical weather models suggest that monsoon precipitation mainly accumulates at the southwestern slopes near the catchment outlet at low elevation, while winter precipitation mainly accumulates along highelevation southern-eastern slopes (Collier and Immerzeel, 2015). Winter westerly events can also provide significant snowfall. Snow cover has strong seasonality with extensive, but sometimes erratic, winter snow cover and retreat of the snowline to higher elevations during spring and summer. For the upper part of the catchment (upstream of Kyangjin) it has been estimated that snowmelt contributes up to 40 % of total runoff (Ragettli et al., 2015).
section_name: Study area

section_content: Remotely sensed snow cover, in situ observations and a modified version of the seNorge snow model were combined to estimate SWE and snowmelt runoff dynamics. The remotely sensed snow cover (Landsat 8 and MOD10A2 snow maps) was first validated with in situ snow cover observations provided by surface temperature and snow depth measurements. 

The snow model was used to simulate daily values of SWE and runoff and was forced by daily in situ meteorological observations of precipitation, temperature and incoming shortwave radiation. MOD10A2 snow cover and snow depth measurements were assimilated to obtain optimal model parameter values using an ensemble Kalman filter (EnKF; Evensen, 2003). The optimized parameters were used for a simulation without assimilation of the observations (open loop). Finally, the model outcome was validated with observed snow depth and Landsat 8 snow cover.
section_name: Calibration and validation strategy

section_name: Datasets

section_name: Remotely sensed snow cover

section_content: MOD10A2 is a Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover product available at http: //reverb. echo. nasa. gov/. The online sub-setting and reprojection utility was used to clip and project imagery for the Langtang catchment. MOD10A2 provides the 8-day maximum snow extent with a spatial resolution of ∼ 500 m. If there is one snow observation within the 8-day period, then the pixel is classified as snow. The 8-day maximum extent offered a good compromise between the temporal resolution and the interference of cloud cover. The snow mapping algorithm used is based on the normalized difference snow index (NDSI; Hall et al., 1995). The NDSI is a ratio of reflection in shortwave infrared (SWIR) and green light (GREEN) and takes advantage of the properties of snow -i. e., snow strongly reflects visible light and strongly absorbs SWIR -Eq. ( 1 ):

The NDSI is calculated with MODIS spectral bands 4 (0. 545-0. 565 µm) and 6 (1. 628-1. 652 µm). Pixels are classified as snow when the NDSI ≥ 0. 4. Water and dark targets typically have high NDSI values, and, to prevent pixels from being incorrectly classified as snow, the reflection should exceed 10 and 11 % for spectral bands 2 (0. 841-0. 876 µm) and 4 respectively for a pixel to be classified as snow (Hall et al., 1995). A full description of the snow mapping algorithm is given by Hall et al. (2002).
section_name: MOD10A2

section_content: Landsat 8 imagery from 15 April 2013 to 5 November 2014 was downloaded from http://earthexplorer. usgs. gov/. Cloudfree scenes (10 out of 34), based on visual inspection, were used to derive daily snow maps with high spatial resolution (30 m). For each image digital numbers were converted to top of atmosphere reflectance. For Landsat 8 the NDSI was calculated with Eq. (1) with spectral bands 3 (0. 53-0. 59 µm) and 6 (1. 57-1. 65 µm). The chosen threshold value was equal to that used for the MOD10A2 snow cover product. The NDSI has proven to be a successful snow mapping algorithm for various sensors with a threshold value around 0. 4 (Dankers and De Jong, 2004). Although the spectral bands have slightly different band widths and spectral positions, a threshold value of 0. 4 gave satisfactory results when compared with in situ snow observations. In addition, the reflection in near-infrared light should exceed 11 % to prevent water from being incorrectly classified as snow (Dankers and De Jong, 2004). Therefore, a pixel is classified as snow when the NDSI value ≥ 0. 4 and the reflectance in near-infrared light > 11 %.
section_name: Landsat 8

section_content: Different types of snow and meteorological observations were available for the study period (January 2013-September 2014; Table 1, Fig. 1 ). Two transects of surface temperature measurements on a north-and south-facing slope provided information on snow cover. The 13 temperature sensors (Hobo Tidbits) were positioned on the surface and covered by a small cairn and recorded surface temperature with 10 min sampling intervals. sonic ranging sensors at four locations at 15 min intervals. Hourly measurements of snow depth were also made at the Kyangjin and Yala base camp automatic weather stations (AWS K and AWS Y; Fig. 1 ). Hourly means (or totals) of air temperature, liquid and solid precipitation, and incoming shortwave radiation were also recorded at AWS Kyangjin (Shea et al., 2015). Air temperature data were also acquired at several locations with 10 and 15 min recording intervals.
section_name: In situ observations

section_content: The snow model was forced with daily average and maximum air temperature, cumulative precipitation and average incoming shortwave radiation for the time period January 2013-September 2014. Hourly measurements of air temperature, precipitation and incoming shortwave radiation at AWS Kyangjin (Shea et al., 2015) were therefore aggregated to daily values. This study period was chosen based on availability of forcing data and observations. Daily temperature lapse rates were interpolated from the air temperature measurements throughout the catchment and used to extrapolate (average and maximum) daily air temperature observed at AWS Kyangjin (Fig. 1 ). The derived temperature lapse rates agree with the values found by Immerzeel et al. (2014). The daily observed precipitation and temperature lapse rates were corrected in the modified seNorge snow model with the correction factors P and T lapse respectively to account for the uncertainty related to undercatch and the derived temperature lapse rates (Table 2 ). Although temperature has a strong relation with altitude and can be accurately derived from multiple weather stations at different altitudes, small differences in the temperature lapse rate (e. g., 0. 001 • C m -1 ) can result in temperature differences of up to several degrees at high altitude in Langtang due to the extreme topography (Immerzeel et al., 2014). Hence, there is a need to consider a potential correction on the temperature lapse rate. A correction is also applied to the daily observed precipitation as precipitation measurements are typically biased due to wind-induced undercatch, especially for solid precipitation (Wolff et al., 2015). Collier and Immerzeel (2015) modeled the spatial distribution of precipitation in Langtang using an interactively coupled atmosphere and glacier mass balance model (Collier et al., 2013). Their study revealed seasonally contrasting spatial patterns of precipitation within the catchment. Monthly modeled precipitation fields from this study were therefore normalized and used to distribute the observed precipitation at AWS Kyangjin. Similarly, a radiation model (van Dam, 2001; Feiken, 2014) was used to extrapolate observed incoming shortwave radiation. The radiation model takes into account the aspect, slope, elevation and shading due to surrounding topography. 

The model initial conditions for January 2013 (i. e., SWE and snow depth) were set by simulating year 2013 three times.
section_name: Model forcing

section_content: The seNorge snow model (Saloranta, 2012 (Saloranta,, 2014 (Saloranta,, 2016) ) is a temperature-index model which requires only data of air temperature and precipitation. In addition, the seNorge snow model includes a compaction module that can be used to assimilate and validate snow depth rather than snow cover only. The low data requirements and the compaction module make the seNorge snow model suitable for application in this study. 

The seNorge snow model was rewritten from its original code into the environmental modeling software PCRaster Python (Karssenberg et al., 2010) to allow spatiotemporal modeling of the SWE and runoff within the catchment. The snow is modeled as a single homogeneous layer with a spatial resolution of 100 m and a daily time step. The seNorge model was further improved by implementing a different melt algorithm, albedo decay and avalanching. These novel model components are described hereafter, and the model parameters used are given in Table 2.
section_name: Modified seNorge model

section_content: Precipitation in the model is partitioned as rain or snow based on an air temperature threshold thr snow ( • C). The snowpack E. E. Stigter et al. : Assimilation of snow cover and snow depth into a snow model consists of a solid component and possibly a liquid component. Meltwater and rain can be stored within the snowpack until its water holding capacity is exceeded and has the possibility to refreeze within the snowpack. The original melt algorithm of the seNorge snow model is substituted by the enhanced temperature-index approach (Pellicciotti et al., 2005 (Pellicciotti et al.,, 2008)). When air temperature (T ; • C) exceeds the temperature threshold for melt onset (T T ; • C), the potential melt (M pot ; mm d -1 ) is calculated for each pixel by Eq. ( 2 ):

where

) is a temperature melt factor, α (-) is the albedo of the snow cover and R inc (W m -2 ) is the incoming shortwave radiation. In case that the threshold temperature is negative, the potential melt can become negative when the radiation melt component is not positive enough to compensate for the negative temperature melt component. When the potential melt is negative it is set to zero to prevent negative values. 

The simulated runoff in the seNorge snow model is the total runoff, i. e., the sum of snowmelt and rain. As the focus of this study is on snowmelt runoff it is necessary to split the runoff in snowmelt and rain runoff. Meltwater and rain fill up the snowpack until its water holding capacity is exceeded. The surplus is defined as snowmelt and rain runoff respectively. If both rain and snowmelt occur it is assumed that rain saturates the snowpack first. Rain falling on snow-free portions of the basin is included in the rain runoff totals.
section_name: Water balance and snowmelt

section_content: Decay of the albedo of snow is calculated with the algorithm developed by Brock et al. (2000) in which the albedo is a function of cumulative maximum daily air temperature T max ( • C). When T max is above 0 • C the air temperature is summed as long as snow is present and no new snow has fallen. When T max is below 0 • C the albedo remains constant. Albedo decay is calculated differently for deep snow (SWE ≥ 5 mm) and shallow snow (SWE < 5 mm). The albedo decay for deep snow is a logarithmic decay, whereas the decay for shallow snow is exponential. This results in a gradual decrease of the albedo for several weeks, which agrees with reality (Brock et al., 2000). When new snow falls the albedo is set to its initial value. In Langtang the observed albedo of fresh snow is 0. 84 and the observed minimum precipitation rate to reset the snow albedo is 1 mm d -1 (Ragettli et al., 2015).
section_name: Albedo decay

section_content: After snowfall events, avalanching occurs regularly on steep slopes in the catchment. Therefore, snow transport due to avalanching is considered to be an important process for redistribution of snow in the Langtang catchment (Ragettli et al., 2015). Snow avalanching is implemented in the model us-ing the SnowSlide algorithm (Bernhardt and Schulz, 2010). For each cell a maximum snow holding depth SWE max (m), depending on slope S ( • ), is calculated using an exponential regression function following Eq. ( 3 ):

where SS 1 and SS 2 are empirical coefficients. If SWE exceeds SWE max and the slope exceeds the minimum slope S min for avalanching to occur, then snow is transported to the adjacent downstream cell. Snow can be transported through multiple cells within one time step. 

As the snowpack is divided into an ice and liquid component, both the ice and liquid components should be transported downwards. Avalanches in the Langtang catchment mainly occur at high elevations where temperatures are low and (almost) no liquid water is present in the snowpack. It is therefore assumed that avalanches are dry avalanches and that no liquid water is present in the avalanching snow. When there is, in rare circumstances, liquid water present in avalanching snow, the liquid water is converted to the ice component to ensure water balance closure.
section_name: Avalanching

section_content: The compaction module is described in detail in Saloranta (2014 Saloranta (, 2016)). In this module SWE is converted into snow depth. Change in snow depth occurs due to melt, new snow and viscous compaction. The change in snow depth due to new snow is adapted such that an increase in snow depth can occur due to both snowfall and deposition of avalanching snow. The increase in snow depth due to deposition of avalanching snow is calculated using a constant snow density for dry avalanches (200 kg m -3 ; Hopfinger, 1983).
section_name: Compaction and density

section_name: Data assimilation

section_content: In order to assess which model parameters to calibrate, a local sensitivity analysis was performed by varying the value of one parameter at a time while holding the values of other parameters fixed. This gives useful first order estimates for parameter sensitivity, although it cannot account for parameter interactions. Plausible parameter values were based on the literature (Table 2 ). The model was run in Monte Carlo (MC) mode with 100 realizations for each parameter. The values for the parameters were randomly chosen from a uniform distribution with defined minimum and maximum values for the parameters. The snow extent and snow depth were averaged over the study period and study area for the sensitivity analysis. The sensitivity of the modeled mean snow extent and mean snow depth were compared to the changes in parameter values. A pixel is determined to be snow covered in the model when the simulated SWE exceeds 1 mm. All the parameters were varied independently per run, except for the melt factors (F T and F SR ), as these are known to be dependent on each other (Ragettli et al., 2015). Therefore, F T and F SR were varied simultaneously in the sensitivity analysis using a linear relation between these melt factors.
section_name: Sensitivity analysis

section_content: Using the ensemble Kalman filter (Evensen, 1994), data assimilation of snow extent and snow depth observations was used to calibrate model parameters using the framework developed by Wanders et al. (2013). Both the EnKF and particle filter (PF) have been used in several studies to assimilate snow observations into snow models (e. g., Charrois et al., 2016; Leisenring and Moradkhani, 2011; Liu et al., 2013; Magnusson et al., 2016). The EnKF and PF are similar in their approach (estimate the model uncertainty from the particle or ensemble spread). The EnKF can only be used for assimilation of continuous values and not for binary values (i. e., snow cover present or not). Therefore, it is necessary to assimilate snow extent (continuous values) into the model, which results in a partial loss of spatial information of snow cover. However, the EnKF has a higher efficiency when it deals with Gaussian data and related errors. The computational demand required for a PF exceeds the EnKF's computer requirements, due to the need to cover the entire (non-Gaussian) distribution. When the number of particles becomes too low, there is an additional risk of particle collapse, especially when one wants to take into account all the grid cells in the simulation with or without snow. This would require a total particle number exceeding the total number of grid cells in the domain, in combination with all the possible parameter combinations to avoid collapse of the filter. For a single site or small sites a PF would be a good alternative (e. g., Charrois et al., 2016; Magnusson et al., 2016), but, limited by the current available computational power, this is only feasible with an EnKF implementation. As we deal with continuous values, it is computationally efficient and allows for dual-state parameter estimations. The lower number of ensemble members compared to a PF allowed us to run multiple simulations over longer time periods, providing a better estimate of the potential of the EnKF improvements. 

An advantage of the EnKF calibration framework is that it allows for the obtaining of an uncertainty estimate for the calibrated parameters. The EnKF obtains the simulation uncertainty by using an MC framework, where the spread in the ensemble members represents the combined uncertainty of parameters and input data. Unfortunately, the EnKF does not allow us to reduce and estimate the model structure uncertainty, since it relies on the assumption that the ensemble members are normally distributed. This assumption is no longer valid if multiple model schematizations are used. Therefore, it is assumed that the model is capable of accurately simulating the processes, when provided with the correct parameters. Besides the parameter and model uncertainty, there is uncertainty in the observations which are as-similated. The EnKF finds the optimal solution for the model states and parameters, based on the observations and modeled predicted values and their respective uncertainties. With sufficient observations the parameters will convert to a stable solution with an uncertainty estimate that is dependent on the observations error and the ability of the model to simulate the observations. It was found that 50 ensemble members are sufficient to obtain stable parameter solutions and correctly represent the parameter uncertainty. 

The EnKF was applied for each time step that observations were available. The MOD10A2 snow extent was divided into six elevation zones. The snow extent per elevation zone was derived from the MOD10A2 snow cover and used for assimilation to include more information on spatial distribution of snow. The elevation zone breakpoints are at 3500, 4000, 4500, 5000 and 5500 m a. s. l. Snow maps with more than 30 % cloud cover and with obvious misclassification of snow were exempted from assimilation (3 snow maps out of 88). Only for cloud-free pixels, comparisons were made between modeled and observed snow extent. Two snow depth observation locations (Pluvio Langshisha and AWS Kyangjin; Fig. 1 ) were also assimilated. 

The EnKF framework allows for the inclusion of an uncertainty in the assimilated observations. Point snow depth measurements have high uncertainties that are related to limited representativeness of point snow depth observations in complex terrain due to local influence of snow drift (Grünewald and Lehning, 2015). For the snow depth measurements a variance of 25 cm was chosen to represent the uncertainty of point snow depth measurements. The MOD10A2 snow extent was assigned an uncertainty based on the classification accuracy (fraction of correctly classified pixels) determined with the in situ snow observations (Sect. 3. 1. 2). The uncertainty is dependent on the snow extent (SE; m 2 ), i. e., an increase in uncertainty for an increase in snow extent. To prevent the uncertainty from becoming zero when there is no snow cover, the minimum variance for each zone was restricted to the average snow extent SE zone (m 2 ) × the accuracy (-). Therefore, the variance σ 2 per elevation zone is defined following Eq. ( 4 ):

The four most sensitive parameters (T T, T lapse, P and C 6 ) resulting from the sensitivity analysis were optimized based on the assimilation of snow depth and MOD10A2 snow extent. The first three parameters (T T, T lapse and P ) influence both snow depth and snow extent and were optimized by assimilating MOD10A2 snow extent. The fourth parameter (C 6 ) is an empirical coefficient relating viscosity to snow density and only influences snow depth. C 6 was optimized by assimilating snow depth observations and taking into account the full uncertainty in the previously determined parameters. 

The two-step approach was chosen to restrict the degrees of freedom and to prevent unrealistic parameter estimates. 3 ). 

Figure 2 shows the monthly cumulative precipitation and the average daily maximum temperature per month measured at AWS Kyangjin for the study period. These data are also available for the time period 1988-2009 and are used to characterize the climatology of the catchment. Comparison of the measurements of the 1988-2009 period and the study period shows that the maximum temperature is similar for both time periods, whereas more variability exists in the cumulative precipitation. Especially in October, a large difference exists in cumulative precipitation, which is caused by a large precipitation event of approximately 100 mm during the study period.
section_name: Parameter calibration

section_name: Results and discussion

section_name: Validation of snow maps with in situ observations

section_content: Surface temperature is an indirect measure of presence of snow. Figure 3 shows observed surface temperature for two locations. Snow cover is distinguishable based on the low diurnal variability in surface temperature when snow is present due to the isolating effect of snow (Lundquist and Lott, 2008 ). An optimal threshold for distinguishing between snow and no snow was determined to be a 2 • C difference between daily minimum temperature and maximum temperature. The use of a larger temperature interval as threshold value was Figure 3. Observed surface temperature with 10 min interval at two locations (Table 1 ). The blue vertical lines indicate the start and end of the snow cover. 

explored; however, as diurnal temperature variability is small during monsoon (Immerzeel et al., 2014) setting, the diurnal cycle temperature threshold above 2 • C may result in incorrect monsoon snow observations.
section_name: In situ snow observations

section_content: Both observed surface temperature and snow depth measurements were converted to daily and 8-day maximum binary snow cover values to validate Landsat 8 and MOD10A2 snow cover respectively. We find that the classification accuracy of MOD10A2 and Landsat 8 snow maps based on all in situ snow observations is 83. 1 and 85. 7 % respectively. The classification accuracy is defined as the number of correctly classified pixels divided by the total number of pixels. Table 4 shows the confusion matrices. Misclassification can be a result of variability of snow conditions within a pixel and classification of ice clouds or high cirrus clouds as snow (Parajka and Blöschl, 2006). Large viewing angles, and consequently larger observation areas, may also result in misclassification (Dozier et al., 2008). MOD10A2 has a lower spatial resolution than Landsat 8 which likely causes the slightly lower accuracy for the MOD10A2 snow cover product (Hall et al., 2002). Visual inspection of MOD10A2 snow maps also revealed that some clouds are erroneously mapped as snow cover. 

The accuracy of MODIS daily snow cover products are reported to be 95 % for mountainous Austria (Parajka and Blöschl, 2006) and 94. 2 % for the upper Rio Grande basin (Klein and Barnett, 2003). The lower accuracy presented in this study is likely a result of the simplification of the 8-day composite product and more extreme relief and consequently larger spatial variability in snow cover. Besides classification errors, uncertainty in the in situ snow observations should be considered as well. For the in situ snow cover observations provided by surface temperature, there are relatively many observations for which snow is not observed in situ, while the MOD10A2 and Landsat 8 snow maps indicate that snow should be present (Table 5 ). This may be caused by the fact that a thin snow layer may not result in sufficient isolation to reduce the diurnal temperature fluctuations for observation as snow (Lundquist and Lott, 2008). This observation bias in the temperature-sensed snow cover data would indicate that MOD10A2 and Landsat 8 snow maps possibly have even higher accuracies than presented here based on this validation approach.
section_name: Remotely sensed snow cover

section_content: The results of the sensitivity of mean snow extent and mean snow depth to parameter variability are shown in Table 2. The sensitivity analysis shows that the threshold temperature for melt onset (T T ), precipitation bias (P ), temperature lapse rate bias (T lapse ) and the coefficient for conversion for viscosity (C 6 ) are the most sensitive parameters. For the snow com- paction parameters, snow depth is most sensitive for changes in C 6, which is in agreement with Saloranta (2014). The melt parameters F SR and F T influence melt directly but show small sensitivity, as these parameters are dependent on each other. A higher value for F T coincides with a lower value for F SR where the value of both parameters is climate zone dependent (Ragettli et al., 2015). 

Only the four most sensitive parameters were chosen to be calibrated by the EnKF to limit the degrees of freedom and to prevent the absence of convergence in the solutions for the parameters. Table 6 shows the prior and posterior parameter distribution resulting from the assimilation of snow extent per zone and snow depth. The parameter values for T lapse, P and C 6 show a narrow posterior distribution (i. e., small standard deviation) indicating that parameter uncertainty is small. T lapse and P represent measurement uncertainties of the model inputs. After calibration the modeled precipitation is increased and the temperature lapse rate is slightly steeper (more negative) than derived. The calibrated value of T T shows a large standard deviation indicating absence of convergence in parameter solutions. This can be either a result of insufficient data to determine the parameter value or insensitivity of the model to the parameter value. A negative value for T T is plausible as melt can occur with air temperatures below 0 • C when incoming shortwave radiation is sufficient. Especially at low latitudes and high elevation, solar radiation is an important cause of snowmelt (Bookhagen and Burbank, 2010). T T is reported to be as negative as -6 • C for Pyramid Station, Nepalese Himalayas (Pellicciotti et al., 2012). Here T T lies in a range which is even more negative than -6 • C. This is likely to be partly a result of the model structure. When T T is negative the melt algorithm (Eq. 2) can give negative values. The temperature term in Eq. (2) becomes negative in case the air temperature is below zero degrees but higher than T T. The reason for negative melt to occur in a few rare cases is a limitation of the EnKF calibration in combination with the enhanced temperature-index method. The EnKF does not allow us to constrain parameter ranges and this results in a relative low T T, which may occasionally lead to negative melt when incoming shortwave radiation is low and the air temperature is above T T. In those cases when negative melt occurs, it is capped to zero, and as a results the model is relatively insensitive for low temperatures close to the T T and the EnKF does not converge into a parameter solution.
section_name: Model calibration

section_name: Model validation

section_content: Both the modeled and MOD10A2 snow extent show strong seasonality of snow cover in the catchment (Fig. 4 ). After calibration, modeled snow extent shows notable improvement in elevation zone 3500-4000 m a. s. l. during the melt season in 2014. After calibration the threshold temperature for melt onset is lower, resulting in more and earlier onset of snowmelt. Consequently there is a decreased snow extent. The zones in the lower areas are expected to show most improvement, as this is the area where snow cover is ephemeral, and considerable improvements of the modeled snow extent in elevation zone 3500-4000 m a. s. l. are indeed observed (Fig. 4 ). 

The root mean square error (RMSE) decreased from 14. 2 to 11. 2 km 2 after calibration. The simulated snow extent agrees well with MOD10A2-observed snow cover for the higher elevation zones (> 4500 m a. s. l. 

). An exception is the snow extent in summer 2013 in the elevation zone 5000-5500 m a. s. l. The snow model underestimates the snow extent compared to the MOD10A2 snow extent. This discrepancy is possibly the result of (i) overestimation of simulated melt, (ii) an actual snow event that is simulated as rain by the model due to toohigh air temperature or (iii) erroneous mapping of clouds as snow in the MOD10A2 snow cover. 

The model classification accuracy of snow cover after calibration is 85. 9 % based on pixel comparison between modeled 8-day maximum snow extent and MOD10A2 snow extent. The classification accuracy is the average classification accuracy over all members. There is only a slight increase of 0. 2 % in accuracy after calibration; however, the performance was already high (85. 7 %) before calibration. The classification accuracy is lower on steep slopes where avalanching is common, and as the snow extent in avalanching zones is highly dynamic, this is not well captured in the model. Calibration of parameters that influence avalanching might overcome this discrepancy to some degree; however, a more advanced approach to avalanche modeling may be required. In addition, the spatial resolution of the remotely sensed snow cover is likely to be insufficient to detect the avalanche dynamics. Other potential explanations for lower classification accuracies are uncertainties related to the simulated precipitation phase (rain or snow) and the simulated spatial distribution of precipitation based on Collier and Immerzeel (2015). 

Landsat 8-derived snow extent is lower in winter than the modeled snow extent and the MOD10A2 snow extent (Fig. 4 ). Distinct differences between the Landsat 8 instantaneous snow cover observations and the MOD10A2 8-day maximum snow cover extents (Fig. 4 ) can be attributed to (i) the sensitivity of the Landsat 8 snow cover maps to misclassified snow pixels in the shaded area, (ii) the much higher spatial resolution of Landsat 8 (Hall et al., 2002), and (iii) the difference between an instantaneous image and an 8-day composite. 

The model classification accuracy, based on pixel comparison with Landsat 8 snow maps, increased from 74. 7 to 78. 2 % after calibration. In Table 7 individual model classification accuracy is given based on comparison with each Landsat 8 snow map. Relative low accuracies occur in winter (especially on 20 December 2013 and 5 January 2014), and the model overestimates snow cover compared to the Landsat 8 snow maps (Fig. 4 ). The overestimation of snow cover by the model on 20 December 2013 is particularly large, and it can be explained by a small snow event (2. 3 mm measured at Kyangjin) a few days before the acquisition. With below zero temperatures the model simulates a large snow cover extent, but based on a very small amount. Snow redistribution by wind, a patchy snow cover and/or sublimation may also explain the mismatch with the Landsat 8 snow cover in this particular case.
section_name: Snow cover

section_content: The observed and modeled snow depths at four locations are shown in Fig. 5. The simulated snow depth is given for the model simulations (i) without calibration, (ii) after calibration of snow extent, and (iii) after calibration of both snow extent and snow depth. After calibration with snow extent there is an increase in snow depth for Yala Pluvio and Yala BC for the entire snow season as result of increased simulated precipitation. For Langshisha and Kyangjin the snow depth mainly decreased after calibration with snow extent. These stations are at a lower elevation, and, since the threshold temperature for melt onset is lowered after calibration, this leads to reduced snow depth. At all locations the modeled snow depth decreased after calibration with both snow extent and snow depth due to lowering of the parameter relating snow density to snow depth. After calibration with both snow extent and snow depth, comparison of modeled and observed snow depth at Langshisha shows good agreement. Especially after calibration, the timing of the melt onset during spring is improved. For Yala Pluvio and Yala BC the agreement between modeled and observed snow depth is also good, though improvement of the timing of melt onset is limited. For Kyangjin the modeled snow depth does not agree as well with observed snow depth in spring 2013, but it improves in 2014. In spring the snow cover duration of snow events decreases after calibration and improves the fit with the observed snow depth. Yala Pluvio and Yala BC are the only locations that serve as an independent validation of snow depth, as these stations are not used for the assimilation. The simulated melt onset in spring is later compared to what is observed. The diurnal variability of air temperature is high during the pre-monsoon season (March to mid-June; Immerzeel et al., 2014). Though daily average air temperatures are below zero, positive temperatures and snowmelt can occur in the afternoon above 5000 m a. s. l. (Shea et al., 2015; Ragettli et al., 2015). This can explain the difference between simulated and observed melt onset. Using an hourly time step might therefore improve the simulation of snowmelt in spring (Ragettli et al., 2015). While the timing of snowpack depletion at Yala Pluvio and Yala BC are offset from the observations, the modeled quantity of snow is in the same order of magnitude for both modeled and observed time series. Hence, there is no substantial overestimation or underestimation of snow depth. The RMSE between simulated and observed snow depth decreases after calibration with both snow extent and snow depth compared to the uncalibrated simulation of snow depth and after calibration of only snow extent. This shows the benefit of assimilating both snow extent and snow depth into the snow model to obtain optimal parameter values. 

While this study shows an approach in using snow depth observations for assimilation and validation, only four locations with snow depth observations were available. The number of available snow depth observations and the choice of different stations for assimilation might influence the results. Four snow depth observations are insufficient for systematic assimilation and independent validation. However, our approach is useful and is recommended for future studies in the Himalayas, in particular when more point observations of snow depth are available.
section_name: Snow depth

section_content: The cumulative basin-wide mean snowfall is 1222 mm for the simulation period. Nearly one-third (31. 4 %) of the snowfall is transported to lower elevations due to avalanching, and 16. 2 % of the snowfall is transported to elevations lower than 5000 m a. s. l. Transport of snow to lower elevations contributes to snowmelt runoff and has been estimated to be 4. 5 % of the total water input for the upper part of the Langtang catchment (Ragettli et al., 2015). 

The simulation of the SWE for the study period shows a pattern of increasing SWE with increasing elevation ( and 7). At higher elevation, air temperature is lower with more snow accumulation than melt, resulting in a higher gain in SWE over time. The glaciers Langtang and Langshisha are positioned at approximately the same elevation (Ragettli et al., 2015), though the SWE is considerably higher at the Langshisha glacier (Fig. 6 ) due to the precipitation distribution approach we use. Also, some areas at higher elevation 3 ). 

Figure 8. Change in SWE averaged over the simulation period and all members for each climate sensitivity test (Table 3 ). 

show less SWE than surrounding areas at the same elevation. These areas represent the steep slopes in the catchment where avalanching occurs regularly. The transported snow accumulates below these steep slopes. The simulated avalanches are based on a simple model parameterization in which the snow is transported via single stream paths, resulting in a few pixels with extreme accumulation of SWE. This is mainly visible in the northeastern part of the catchment. Modeling the divergence of transported snow might improve the extreme accumulation simulated for some pixels. 

For the climate sensitivity tests a delta-change method is used. This method has limitations as climate variability of future climate is not constant compared to the study period (Kobierska et al., 2013). In addition, Kobierska et al. (2013) showed that changes in runoff due to climate change are predicted differently by a physically based snow model and a parameterized snow model for a glacierized catchment. Parameterized snow models (such as the modified seNorge snow model that is used in this study) are calibrated to fit the current climate and not future climate and might therefore be incapable of predicting future states of the snowpack. However, the scope of this study is to show the sensitivity of the SWE and snowmelt runoff to changes in air temperature and precipitation and not that of a full-fledged climate impact study. Therefore, the use of a parameterized snow model and the delta-change method is suitable in this case. All climate sensitivity tests show a decrease in SWE, but the relative change is greatest at low elevations in the valley. We also observe a strong gradient of decreased relative change in SWE with increased elevation. An increase in temperature leads to an increase in melt and more precipitation in the form of rain instead of snow. Both processes result in decreased relative change of SWE with elevation. Near the catchment outlet there is an area with 100 % decrease in SWE, as precipitation will only fall as rain instead of snow. 

A slight deviation from the elevational trend in SWE change occurs between 3000 and 4000 m a. s. l., which is a zone that could be sensitive to changes in the elevation at which snowfall occurs. The combination of snowfall at higher elevations due to higher temperature and the monthly differing spatial patterns in precipitation are likely to explain the banded patterns. 

Changes in SWE and the spatial distribution of SWE will also be affected by changes in total precipitation. The influence of precipitation can be determined based on comparison of the two wet and dry climate sensitivity tests (Figs. 7 and 8 ). A decrease in precipitation results in decreased SWE as there is less snowfall. However, the increased precipitation for the wet/cold and wet/warm climate sensitivity tests (+12. 1 and +12. 4 % respectively) does not compensate for the temperature-related increase in melt and decrease in snowfall in the valley. 

Reduced warming under the wet/cold climate sensitivity test results in a smaller decrease of SWE compared to the wet/warm climate sensitivity test, even in the valley. At higher elevations, changes in SWE are weakly negative and in some areas positive. Snowpack sensitivity to temperature change decreases with elevation (Brown and Mote, 2009). The increased SWE under both wet climate sensitivity tests occurs in the southeastern part of the catchment where relatively large amounts of precipitation occur in winter (Collier and Immerzeel, 2015). Schmucki et al. (2015) showed similar results for the Alps. They showed that low-and midelevation stations are sensitive to temperature change but not to precipitation change. In contrast, at high-elevation stations an increase in precipitation partly compensates for an increase in temperature. The compensating effect of increased precipitation at high elevations is important for glacier systems and emphasizes the importance of accurate estimations of both change in precipitation and its spatial distribution. 

The modeled snowmelt and rain runoff at the catchment outlet is greatest during the monsoon and lowest during winter (Fig. 9 ). Peak snowmelt and rain runoff occur in June and July respectively. The snowmelt season starts in March when temperatures and insolation are rising and continues until October. Snowmelt runoff contributes most to total runoff during pre-monsoon and early-monsoon (March-June), which is in agreement with Bookhagen and Burbank (2010). Validation of the simulated runoff with observed runoff was impossible, because (i) there were no reliable runoff data available for the study period, as there was no reliable rating curve, and (ii) the model focusses on rain and snowmelt runoff; however, glacier runoff and delay of runoff due to groundwater and glacier storage is not incorporated in the model structure. 

The climate sensitivity of snowmelt and rain runoff is shown in Fig. 9. All climate sensitivity tests show an increase in snowmelt runoff from October to May. In contrast, snowmelt runoff decreases from June to September. Higher temperatures result in more snowmelt and less snowfall during winter and an early melt season which leads to a shift in the peak of snowmelt runoff. In other mountain regions similar changes in runoff patterns appear. Several studies in the Alps show that the peak in snowmelt runoff shifts from summer to late spring (Bavay et al., 2009 (Bavay et al.,, 2013;; Kobierska et al., 2013). Immerzeel et al. (2009) showed that in the upper Indus Basin, the peak in snowmelt runoff appears 1 month earlier by 2071-2100 as result of an increase in temperature and precipitation. However, Immerzeel et al. (2012) showed that total snowmelt runoff remains more or less constant under positive temperature and precipitation trends in the upper part of the Langtang catchment. In their study snowmelt on glaciers is not defined as snowmelt runoff and is therefore a minor component of total runoff, leading to different results. 

For the wet climate sensitivity tests, total runoff (i. e., the sum of snowmelt and rain runoff) increases throughout the year. The decrease in melt runoff during the late melt season is compensated by the increase in rain runoff as there is more precipitation. The future hydrology of the central Himalayas largely depends on precipitation changes, as it is dominated by rainfall runoff during the monsoon (Lutz et al., 2014). As we perturb the model with a percentage change in precipitation that is constant through the year, the absolute change in precipitation is greater in the monsoon than in winter. For climate sensitivity tests with decreased precipitation, total runoff from June to September decreases, but from October to May it increases as a result of increased snowmelt. Estimates of seasonal changes in precipitation are thus critical for determining whether rain and snowmelt runoff increases or decreases during monsoon.
section_name: Climate sensitivity of SWE and snowmelt runoff

section_content: Remotely sensed snow cover, in situ observations and a modified seNorge snow model were combined to estimate (climate sensitivity of) SWE and snowmelt runoff in the Langtang catchment. Validation of remotely sensed snow cover (Landsat 8 and MOD10A2 snow maps) shows high accuracies (85. 7 and 83. 1 % respectively) against in situ snow observations provided by surface temperature and snow depth measurements. Data assimilation of MOD10A2 snow cover and snow depth measurements using an EnKF proved to be successful for obtaining optimal model parameter values. Independent validations of simulated snow depth and snow cover against snow depth measurements and Landsat 8 snow cover show improvement after assimilation of snow depth and snow cover compared to results before data assimilation. The applied methodology of simultaneous assimilation of snow cover and snow depth allows for the calibration of important snow parameters and validation of the snow depth rather than snow cover alone. This opens up new possibilities for future snow assessments and sensitivity studies in the Himalayas. 

The spatial distribution of SWE averaged over the simulation period (January 2013-September 2014) shows a strong gradient of increasing SWE with increasing elevation. In addition, the SWE is considerably higher in the southeastern part of the catchment than the northeastern part of the catchment as a result of the spatial and temporal distribution of precipitation. 

Finally the climate sensitivity study revealed that snowmelt runoff increases in winter and the early melt season (December-May) and decreases during the late melt season (June-September) as a result of the earlier onset snowmelt due to increasing temperature. There is a strong relative decrease in SWE in the valley with increasing temperature due to more snowmelt and less precipitation as snow. At higher elevations an increase in precipitation partly compensates for increased melt due to higher temperatures. The compensating effect of precipitation emphasizes the importance and need for the accurate prediction of change in the spatial and temporal distribution of precipitation.
section_name: Conclusions

section_content: Acknowledgements. This project was supported by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 676819 ) and by the research programme VIDI with project number 016. 161. 308financed by the Netherlands Organisation for Scientific Research (NWO). The authors thank Hendrik Wulf and two anonymous reviewers for their constructive comments that helped in improving the manuscript. 

Edited by: Guillaume Chambon Reviewed by: Hendrik Wulf and two anonymous referees

section_content: Data availability. Data are available upon request to the corresponding author.

section_content: Competing interests. The authors declare that they have no conflict of interest.