Climate sensitivity of snow water equivalent and snowmelt runoff in a Himalayan catchment

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 this provides no 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 meteorological 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. Landsat 8 and MOD10A2 snow cover maps 15 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. The approach of modelling snow depth in a Kalman filter framework allows for data-constrained estimation of SWE rather than snow cover alone and this has great potential for future studies in the Himalayas. Climate sensitivity tests with the optimized 20 snow model show a strong decrease in SWE in the valley with increasing temperature. However, at high elevation a decrease in SWE 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. Finally the climate sensitivity study revealed that snowmelt runoff increases in winter and 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. 25


Introduction
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., 2009Immerzeel 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 (Immerzeel et al., 2010).So far, the main focus has been on effect of climate change on The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -216, 2016 Manuscript under review for journal The Cryosphere Published: 26 October 2016 c Author(s) 2016.CC-BY 3.0 License.
the glaciers and the resulting runoff.However, snow is an important short-term water reservoir in 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 lower 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).Studies about snowmelt in the Himalayas have predominantly relied on remotely sensed snow cover and a modelled melt flux estimating melt runoff resulting from this snow cover (e.g.Bookhagen and Burbank, 2010;Immerzeel et al., 2009;Tahir et al., 2011).However, this provides no information on snow water equivalent (SWE), which is an important hydrologic measure as it indicates the actual amount of water stored in a snowpack.Currently there is no reliable information on SWE for the Himalayas (Lutz et al., 2015).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.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 modelling snow depth allows to validate the quantity of simulated snow rather than snow cover alone and is a new approach in Himalayan snow research.

Study area
The study area is the Langtang catchment, which is located in the central Himalaya approximately 100 km north of Kathmandu (Figure 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 (asl) at the catchment outlet to 7234 m asl for Langtang Lirung, which is the highest peak in the catchment.The climate is monsoon-dominated and 68% to 89% of the annual precipitation falls during the monsoon (Immerzeel et al., 2014).The spatial patterns in precipitation are seasonally contrasting.During the monsoon precipitation mainly accumulates at the southern slopes and near the catchment outlet at low elevation.However, during the winter precipitation mainly accumulates along high-elevation southern slopes (Collier and Immerzeel, 2015).
Winter westerly events can 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 and less snow cover.
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).

Calibration and validation strategy
Remotely sensed snow cover, in situ meteorological 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.

MOD10A2
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 short-wave infrared (SWIR) and visible wavelengths and takes advantage of the properties of snow i.e. snow strongly reflects visible light and strongly absorbs SWIR.The NDSI is calculated with spectral bands 4 (0.545-0.565 μm) and 6 (1.628-1.652μm) following Eq.( 1): (1) 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 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).

Landsat 8
Landsat 8 imagery from 15 April 2013 to 5 November 2014 was downloaded from http://earthexplorer.usgs.gov/.Cloud free scenes, 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%.

In situ observations
Different types of snow and meteorological observations were available for the study period (January 2013 -September 2014; Table 1, Figure 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 temperature with 10 minute sampling intervals.Snow depths were measured with sonic ranging sensors at 4 locations along the transects at 15 minute intervals.Hourly measurements of snow depth were also made at the Kyangjin automatic weather station (AWS; 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 was also acquired at several locations with 10 and 15 minute recording intervals.

Model forcing
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 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.The daily observed precipitation and temperature lapse rates were corrected in the modified seNorge snow model with the correction factors  and   respectively to account for measurement uncertainty (Table 2).Collier and Immerzeel (2015) modelled 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 modelled precipitation fields from this study were therefore normalized and used to distribute the observed precipitation at AWS Kyangjin.Similarly, a radiation model (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.

Modified seNorge model
The seNorge snow model (Saloranta, 2012(Saloranta, , 2014(Saloranta, , 2016) ) was rewritten from its original code into the environmental modelling software PCRaster-Python (Karssenberg et al., 2010) to allow spatio-temporal modelling of the SWE and runoff within the catchment.The snow is modelled 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.

Water balance and snowmelt
Precipitation in the model is partitioned as rain or snow based on an air temperature threshold ℎ  (°C).The snowpack 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 (; °C) exceeds the temperature threshold for melt onset (; °C) the potential melt (  ; mm d -1 ) is calculated by Eq. (2): where  (m 2 mm W -2 d -1 ) is a radiative melt factor,  (mm °C-1 d -1 ) is a temperature melt factor,  (-) is the albedo of the snow cover and   (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.

Albedo decay
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  max (°C).When maximum air temperature is above 0 °C the air temperature is summed as long snow is present and no new snow has fallen.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).

Avalanching
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 using the SnowSlide algorithm (Bernhardt and Schulz, 2010).For each cell a maximum snow holding depth   (m), depending on slope  (°), is calculated using an exponential regression function following Eq.( 3): where  1 and  2 are empirical coefficients.If SWE exceeds   and the slope exceeds the minimum slope   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 component should be transported downwards.Separate transport and redistribution of both the liquid and ice component is beyond the scope of this study.Avalanches in 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.

Compaction and density
The compaction module is described in detail in Saloranta (2012Saloranta ( , 2014Saloranta ( , 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).

Sensitivity analysis
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 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 sensitivity of the modelled mean snow extent and mean snow depth were compared to the changes in parameter values.All the parameters were varied independently per run, except for the melt factors ( and ) as these are known to be dependent on each other (Ragettli et al., 2015).Therefore,  and  were varied simultaneously in the sensitivity analysis using a linear relation between these melt factors.

Parameter calibration
Using the Ensemble Kalman filter (EnKF; 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).An advantage of the EnKF calibration framework is that it allows obtaining an uncertainty estimate for the calibrated parameters.The EnKF obtains the simulation uncertainty by using a MC framework, where the spread in the ensemble members represent the combined uncertainty of parameters and input data.Unfortunately, the EnKF does not allow 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 to accurately simulate the processes, when provided with the correct parameters.Besides the parameter and model uncertainty there is uncertainty in the observations which are assimilated.The EnKF finds the optimal solution for the model states and parameters, based on the observations and modelled 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 6 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 asl.Snow maps with more than 30% cloud cover and with obvious miss-classification of snow were exempted from assimilation (3 snow maps out of 88).Only for cloud free pixels comparisons were made between modelled and observed snow extent.Two snow depth observation locations (Pluvio Langshisha and AWS Kyangjin; Figure 1) were also assimilated.
The EnKF framework allows for the inclusion of an uncertainty in the assimilated observations.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 determined with the in situ snow observations (Sect.3.1.2Remotely sensed snow cover).The uncertainty is dependent on the snow extent  (m 2 ), i.e. an increase in uncertainty for an increase in snow extent.To prevent the uncertainty to become zero when there is no snow cover, the minimum variance for each zone was restricted to the average snow extent  ̅̅̅̅  (m 2 ) times the  (-).
Therefore the variance  2 per elevation zone is defined following Eq.( 4): The four most sensitive parameters (,   ,  and  6 ) resulting from the sensitivity analysis were optimized based on the assimilation of snow depth and MOD10A2 snow extent.The first three parameters (,   and ) influence both snow depth and snow extent, and were optimized by assimilating MOD10A2 snow extent.The fourth parameter ( 6 ) only influences snow depth, and 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.

Climate sensitivity
Climate sensitivity tests were performed to investigate changes in SWE and snowmelt runoff as a result of temperature and precipitation changes.Climate sensitivity was tested by perturbing daily average air temperature and daily cumulative precipitation.The changes in temperature and precipitation were based on an envelope of projected changes in temperature and precipitation for RCP4.5 for the Langtang catchment as in Immerzeel et al. (2013).Four climate sensitivity tests were performed ranging from dry to wet and cold to warm (Table 3).

In situ snow observations
Surface temperature is an indirect measure of presence of snow.due to the isolating effect of snow (Lundquist and Lott, 2008).An optimal threshold for distinguishing between snow/no snow was determined to be 2°C difference between daily minimum temperature and maximum temperature.The use of a larger temperature interval as threshold value was 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.

Remotely sensed snow cover
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).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 more extreme relief and consequently larger 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.

Model calibration
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 (), precipitation bias (), temperature lapse rate bias (  ) and the coefficient for conversion for viscosity ( 6 ) are the most sensitive parameters.For the snow compaction parameters, snow depth is most sensitive for changes in  6 which is in agreement with Saloranta (2014).The melt parameters  and  influence melt directly but show small sensitivity as these parameters are dependent on each other.A higher value for  coincides with a lower value for  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   ,  , and  6 show a narrow posterior distribution (i.e.small standard deviation) indicating that parameter uncertainty is small.  and  account for measurement uncertainty of the model input.After calibration the modelled precipitation is increased and the temperature lapse rate is slightly steeper (more negative) than derived.The calibrated value of  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  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). is reported to be as negative as -6 °C for Pyramid, Nepalese Himalayas (Pellicciotti et al., 2012).Here  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  is negative the melt algorithm (Eq.( 1)) can give negative values.This is prevented by restricting the melt algorithm to a minimum value of zero.This can lead to no melt or refreezing at negative temperatures higher than .The restriction makes the algorithm therefore insensitive for very low temperatures and results in absence of convergence in the parameter solutions.The EnKF however does not restrict the parameter values, which allows  to become too negative.

Snow cover
Both the modelled and MOD10A2 snow extent show strong seasonality of snow cover in the catchment (Figure 3).After calibration modelled snow extent shows notable improvement in elevation zone 3500-4000 m asl 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 modelled snow extent in elevation zone 3500 -4000 m asl are indeed observed (Figure 3).The root mean square error (RMSE) decreased from 14.2 to 11.2 km 2 after calibration.
The model classification accuracy of snow cover after calibration is 85.9% based on pixel comparison between modelled 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 slopes in the northern part of the catchment (Figure 4).This area contains 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 modelling may be required.In addition the spatial resolution of the remotely sensed snow cover is likely to be insufficient to detect the avalanche dynamics.
After calibration the accuracy increased in the lower area of the northern slopes which corresponds to the improvement of snow extent in elevation zone 3500-4000 m asl (Figure 4).
Landsat 8 derived snow extent is lower in winter than the modelled snow extent and the MOD10A2 snow extent (Figure 3).Distinct differences between the Landsat 8 instantaneous snow cover observations and the MOD10A2 8-day maximum snow cover extents (Figure 3) can be attributed to (i) the sensitivity of the Landsat 8 snow cover maps to misclassified snow pixels in 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.3).The overestimation of snow cover by the model on 20-12-2013 is particularly large and it can be explained by a small snow event a few days before the acquisition day of 2.3 mm measured at Kyangjin.Given the below zero temperatures the model reveals a large snow cover extent, but based on a very small amount.Snow redistribution by wind, a patchy snow cover and/or sublimation may therefore explain the mismatch with the Landsat 8 snow cover in this particular case.

Snow depth
The observed and modelled snow depths at three locations are shown in Figure 5.At all locations the modelled snow depth decreased after calibration due to decrease in the threshold temperature for melt onset after calibration.In addition, the parameter relating snow density to snow depth is lowered after calibration, leading to reduced snow depth.Comparison of modelled and observed snow depth at Langshisha shows good agreement.Especially after calibration the timing of the melt onset during spring is improved.For Yala the agreement between modelled and observed snow depth is also good, though improvement of the timing of melt onset is limited.For Kyangjin the modelled snow depth agrees less well with observed snow depth in spring 2013, but it improves in 2014.
Yala is the only location which serves as an independent validation of snow depth as this station is 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 asl (Shea et al., 2015;Ragettli et al., 2015).This explains 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 is offset from the observations, the modelled quantity of snow is in the same order of magnitude for both modelled and observed time series.Hence there is no substantial overestimation or underestimation of SWE.While this study shows an approach in using snow depth observations for assimilation and validation, only three locations with snow depth observations were available.This is 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.

Snow processes
The snow model contains modules with snow processes such as avalanching, albedo decay and compaction.These processes are shown and discussed in this section.The cumulative basin-wide mean snowfall is 1222 mm for the simulation period.
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.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 simulated compaction and albedo decay at Yala (location of snow depth observation) for an accumulation and ablation period (left and right panel respectively) are shown in Figure 6.During the ablation period, snow depth declines through time while snow density shows an inversed similar trend.Reduction of snow depth is a result of both melt and compaction.A decrease in snow depth on days without melt is a result of viscous compaction i.e. compaction due to weight of overlying snow, which increases the snow density.Increases in snow density are greater on days with snowmelt, as snowmelt influences the viscosity of snow due to presence of liquid water in the snow.When liquid water is present within a snowpack the viscosity of the snowpack decreases and enhances compaction (Vionnet et al., 2012).Therefore the decrease in snow depth on days with melt is a result of both compaction and melt.Snow depth increases during the accumulation period due to snowfall.New snow has a lower density than the snowpack and therefore lowers the density of the single layer snowpack.The snow density is for example low in January due to series of snow events, whereas in November a prolonged melt period resulted in a higher snow density.
At low latitudes and high elevation incoming shortwave radiation is an important cause of snowmelt (Bookhagen and Burbank, 2010).The albedo of snow there has a strong impact on snow melt.For the ablation period the albedo initially decays rapidly (Figure 6).Next the albedo remains constant before the albedo declines again.The albedo decay of snow is modelled as a function of cumulative maximum temperature above 0 °C (Sect.2.5.2Albedo decay).The albedo remains constant on days with a maximum air temperature below 0 °C.The rapid decay after snowfall is caused by maximum air temperature above 0 °C and a logarithmic decay function.The logarithmic function also explains the more gradual decay in albedo from 23 November onwards.An initial rapid decrease in albedo after a snow event as well as a more gradual decay over time is characteristic for the decay of albedo (Brock et al., 2000).For the accumulation period the albedo increases occasionally as result of snowfall.

Climate sensitivity of SWE and snowmelt runoff
The simulation of the SWE for the reference scenario shows a pattern of increasing SWE with increasing elevation (Figure 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 Lansgshisha are positioned at approximately the same elevation (Ragettli et al., 2015), though the SWE is considerable higher at Langshisha glacier (Figure 7).Total precipitation is highest along the southern ridge of the catchment and at high elevation in winter (Collier and Immerzeel, 2015) and therefore explains the difference between the SWE at Langtang and Langshisha glacier.Also, some areas at higher elevation 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 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.Modelling the divergence of transported snow might improve the extreme accumulation simulated for some pixels.
Figure 8 shows the results of the change in SWE for different climate sensitivity tests.All climate sensitivity tests show a decrease in SWE, but the change is greatest at low elevations in the valley.We also observe a strong gradient of decreased SWE change with increased elevation.An increase in temperature leads to an increase in melt and more precipitation in form of rain instead of snow.Both processes result in decreased 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 in elevation zone (3000-4000 m asl) that could be sensitive to the transition of 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.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 snow falling 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).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 modelled snowmelt and rain runoff at the catchment outlet is greatest during the monsoon and low during winter (Figure 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 was 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 Figure 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.Immerzeel et al. (2009) showed that in the upper Indus Basin the peak in snowmelt runoff appears one 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 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 to determine whether rain and snowmelt runoff increases or decreases during monsoon.

Conclusions
Remotely sensed snow cover, in situ meteorological 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) show 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 proves to be successful for obtaining optimal model parameter values.The applied methodology of simultaneous assimilation of snow cover and snow depth allows for the calibration of important snow parameters and validation of the SWE 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
Figure 2 shows observed surface temperature for two locations.Snow cover is distinguishable based on the low diurnal variability in surface temperature when snow is present The Cryosphere Discuss., doi:10.5194/tc-2016-216,2016 Manuscript under review for journal The Cryosphere Published: 26 October 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 2 :
Figure 2: Observed surface temperature with 10 minute interval at two locations.The blue vertical lines indicate the start and end of the snow cover.

Figure 3 :
Figure 3: Modelled 8-day maximum snow extent before and after calibration (ensemble mean), Landsat 8 snow extent and MOD10A2 snow extent per elevation zone.The RMSE (km 2 ) is given per zone for the fit between modelled (before and after calibration) and MOD10A2 snow extent.5

Figure 4 :
Figure 4: Accuracy per pixel based on comparison between MOD10A2 snow extent and modelled 8-day maximum snow extent a) before calibration and b) after calibration.

Figure 5 :
Figure 5: Observed snow depth and modelled snow depth before calibration and after calibration (ensemble mean) at three 5

Figure 6 :
Figure 6: Snow processes at Yala.SD is the snow depth, ρ is the density, α is the albedo of the snow cover, M act is the actual snowmelt.The plotted variables represent the ensemble mean.

Figure 7 :
Figure 7: Spatial distribution of ensemble mean annual average snow water equivalent (SWE).5

Figure 8 :
Figure 8: Change in SWE averaged over the simulation period and all members for each climate sensitivity test: a) dry and cold b) dry and warm c) wet and cold d) wet and warm.

Figure 9 :
Figure 9: Modelled runoff at catchment outlet for the reference scenario (January 2013 -September 2014) and change in runoff compared to the reference scenario for the climate sensitivity tests.
The model classification accuracy, based on pixel comparison with Landsat 8 snow maps, increased from 74.7% to 78.2% after calibration.In Table7 individualmodel classification accuracy is given based on comparison with each Landsat 8 snow map.Relative low accuracies occur in winter (especially at 20-12-2013 and 05-01-2014) and the model overestimates snow cover compared to the Landsat 8 snow maps (Figure The CryosphereDiscuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc--216,2016Manuscript under review for journal The Cryosphere Published: 26 October 2016 c Author(s) 2016.CC-BY 3.0 License.

Table 2 : Parameters in the snow model. Value indicates the uncalibrated parameter value and the value range indicates the range which is used for the sensitivity analysis. Sensitivity of snow depth (SD) and snow extent (SE) represent the difference between the 90 th and 10 th percentile of mean snow depth and snow extent resulting from the sensitivity analysis.
The CryosphereDiscuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc--216,2016Manuscript under review for journal The Cryosphere Published: 26 October 2016 c Author(s) 2016.CC-BY 3.0 License.