Hydrological Modelling and data assimilation of Satellite Snow Cover Area using a Land Surface Model, VIC

The snow cover plays an important role in Himalayan region as it contributes a useful amount to the river discharge. So, besides estimating rainfall runoff, proper assessment of snowmelt runoff for efficient management and water resources planning is also required. A Land Surface Model, VIC (Variable Infiltration Capacity) is used at a high resolution grid size of 1 km. Beas river basin up to Thalot in North West Himalayas (NWH) have been selected as the study area. At first model setup is done and VIC has been run in its energy balance mode. The fluxes obtained from VIC has been routed to simulate the discharge for the time period of (2003-2006). Data Assimilation is done for the year 2006 and the techniques of Data Assimilation considered in this study are Direct Insertion (D.I) and Ensemble Kalman Filter (EnKF) that uses observations of snow covered area (SCA) to update hydrologic model states. The meteorological forcings were taken from 0.5 deg. resolution VIC global forcing data from 1979-2006 with daily maximum temperature, minimum temperature from Climate Research unit (CRU), rainfall from daily variability of NCEP and wind speed from NCEP-NCAR analysis as main inputs and Indian Meteorological Department (IMD) data of 0.25 °. NBSSLUP soil map and land use land cover map of ISRO-GBP project for year 2014 were used for generating the soil parameters and vegetation parameters respectively. The threshold temperature i.e. the minimum rain temperature is -0.5°C and maximum snow temperature is about +0.5°C at which VIC can generate snow fluxes. Hydrological simulations were done using both NCEP and IMD based meteorological Forcing datasets, but very few snow fluxes were obtained using IMD data met forcing, whereas NCEP based met forcing has given significantly better snow fluxes throughout the simulation years as the temperature resolution as given by IMD data is 0.5°C and rainfall resolution of 0.25°C. The simulated discharge has been validated using observed data from BBMB (Bhakra Beas Management Board) and coefficient of Correlation(R) measured for (2003-2006) was 0.67 and 0.61 for the year 2006.But as VIC does not consider snowmelt runoff as a part of the total discharge, snowmelt runoff has been estimated for the simulation both with and without D.A. The snow fluxes as generated from VIC gives basin average estimates of Snow Cover, SWE, Snow Depth and Snow melt. It has been observed to be overestimated when model predicted snow cover is compared with MODIS SCA of 500 m resolution from MOD10A2 for each year. So MODIS 8-day snow cover area has been assimilated directly into the model state as well as by using EnKF after every 8 days for the year 2006.D.I Technique performed well as compared to EnKF. R between Model SCA and MODIS SCA is estimated as 0.73 after D.I with Root Mean Square Error (RMSE) of +0.19. After direct Insertion of D.A, SCA has been reduced comparatively which resulted in 7% reduction of annual snowmelt contribution to total discharge.The assimilation of MODIS SCA data hence improved the snow cover area (SCA) fraction and finally updated other snow components.


INTRODUCTION
Snow is one of the key environmental parameter which not only maintains the earth's radiation balance but also plays an important role in river discharge.So proper snowmelt studies is required for efficient management of water resources.Snow Cover Area (SCA), portion of the watershed covered by the snow is one of the vital parameter in terms of snowmelt runoff estimation.It forms a firm basis for improving other snow components i.e.Snow water equivalent (SWE), Snow depth and Snow Melt.Over the past two decades, Satellite remote sensing came into notice as the conventional methods of acquiring snow properties becomes very much difficult in inaccessible mountainous regions.It has opened the possibility of data acquisition at regular intervals, monitoring the existing conditions and also has been considered useful for several runoff forecasting and water management systems (Jain and Rathore, 2001).So this study aims at assimilating the MODIS snow cover product into an Energy balanced Land Surface model known as Variable Infiltration Capacity (VIC).This is known as Data Assimilation.There are many physically based hydrological model as well as Land Surface Models like SWAT, WinSRM , VIC, Noah LSM etc. SWAT partitions the rain and snow on the basis of temperature.WinSRM works on the Degree Day Approach and it divides the whole basin area into zones based on elevation and gives the snowmelt runoff but these models does not allow for Data Assimilation.This is the key reason of choosing VIC for this study.(Huang, 2012) used MODIS Snow cover fraction (SCF) to improve the Snow Water Equivalent (SWE), Snow Depth (SD) and Snow Melt(SM) by EnKF Technique of Data Assimilation.Validation of MODIS snow products and its very reasonably good agreement with the ground observation data from various regions are well known from the literatures of (Klein and Bernett, 2003); (Parajka and Bloeschl 2006); (Wang et al, 2008).(Jain et al, 2008) concluded that MODIS data could be effectively used for SCA estimation under Himalayan rugged topographic and harsh climatic conditions.Review of many other works also reveals how the MODIS snow-cover products have been used and also includes details regarding snow cover mapping accuracy.Though hydrologic remote sensing help in providing valuable information about Land Surface Conditions like SWE, SD, SM and soil moisture etc. but sometimes it may not be enough for many other applications.(Reichle et al, 2008) stated these data may be considered useful when used in Data Assimilation System.Direct Insertion is one of the earliest and most simplistic approaches to data assimilation where the forecast model states are directly replaced with the observations.This approach has been generally applied in many snow data assimilation studies (Reichle et al., 2008).(Rodell et al., 2004) have used direct insertion to assimilate MODIS snow cover observations into the Global Land Data Assimilation System (GLDAS).(Andreadis and Lettenmaier, 2005) have carried out snow data assimilation via an Ensemble Kalman Filter but slight improvement was noticed over the unupdated simulation.(Jin and Miller, 2003) also used this technique to study the snowpack impact on climate variability and snowmelt processes in mountainous regions.(McGuire et al., 2006) used VIC, a macroscale hydrological model in to adjust model's initial snow data using MODIS imagery for winters (2000)(2001)(2002)(2003)(2004) and the insertion of the MODIS data resulted in forecast error reduction.The current study uses the MODIS snow cover product (MOD10A2) to update the snow components by using a land surface hydrological model, VIC by applying the D.I and EnKF technique.The study also aims at estimating snowmelt runoff contribution to the annual discharge at the Thalot site.MOD10A2 comes in a spatial resolution of 500*500 m 2 in every 8-day period that begins on the first day of each year and extends to the first few days of the next year, (Maskey et al., 2011).

STUDY AREA AND DATA USED
The study area chosen for this work is Beas river basin in North West Himalaya (NWH).Beas River in northern India is an important river of Indus River system.It rises in the Himalayas in central Himachal Pradesh, India, and flows for some 470 kilometers to the Sutlej River in the Indian state of Punjab.Its total length is 470 kilometers and its drainage basin is 20,303 square kilometers.The Beas river basin with its outlet at Thalot discharge station is chosen for this study.It extends from 76. 935-77.867°E and 31.507-32.414°N.The basin covers an area of 4883 km 2 out of which 70 km 2 is under snow.The climate here varies with the elevation and is affected by the monsoon.The elevation of the basin till Thalot extents from 938 m to 6500 m approximately where Rohtang covers the higher elevated area.The origin of higher streams and the Beas valley upstream to Kullu is the glaciated area of the basin.So, the discharge of the Beas River basin includes runoff from both melting of the snow and ice and the other resulting from the rainfall in the catchment.
The remote sensing data used in this study are 8-day 500 m resolution MODIS satellite Snow Cover Area data (MOD10A2 SCA product) for updating the model predicted fractional snow cover area.Other geo-spatial data such as Land Use Land Cover (LULC) map, DEM, soil map and the hydro-meteorological data has also been used.Cartosat Dem with resolution of 30 meters has been used for hydroprocessing and extraction of other topographic features like slope map, elevation zonal map, aspect map etc.The meteorological forcings were taken from 0.5 deg.resolution VIC global forcing data (Adam andLettenmaier, 2003 andAdam et al., 2006)   Observed discharge data from Bhakra Beas Management Board (BBMB) has been used.These are all essential required inputs for running the VIC model.

METHODOLOGY
This study uses the land surface model Variable Infiltration Capacity (VIC) which solves both water balance and energy balance, Liang et al., (1994).This model accounts for the heterogeneity of the surface such as subgrid variability of vegetation classes, soil moisture storage capacity as it divides the whole study area into number of grids.Each grid can also be divided into different elevation zones that results in more realistic hydrology in mountainous regions.Study area has been divided into high grid resolution of 1km (1km* 1 km).At first model setup is done and VIC has been run in its energy balance mode.The fluxes obtained from VIC has been routed to simulate the discharge for the time period of (2003)(2004)(2005)(2006).Direct Insertion of Data Assimilation with SCA has been implemented for the year 2006 in this study to improve the model forecasts of other snow components.Keeping in view of other advanced Data Assimilation Techniques, Ensemble Kalman Filter has also been implemented.

Methodology adopted for VIC model
The watershed boundary is delineated from Cartosat DEM of 30 meters resolution and grids were generated over the basin area of high resolution of 1km (0.01*0.01 degree approx.).The soil parameter file generated from NBSSLUP describes the characteristics of each of the considered soil layer for each grid.provides the VIC routing model (Lohmann et al, 1996) the fraction of each grid cell in the study area and the station file provides the exact position of the outlet.Flow direction file shows flow direction value for each grid cell in the format of VIC.Finally, the water and energy balance components which were obtained from the VIC fluxes has been routed to get the discharge.

Assimilation Method
Direct Insertion approach has been used in this study which is the earliest and simplest method for the Data Assimilation.(Evensen, 1994) who gives the idea to use kalman filter for the non-linear systems by using ensembles or statistical examples of the forecasted model states.In the forecast step, the ensemble members are each forecast individually by the full non linear model, and in the analysis step the update equations are based on the linear equations used in the kalman filter.
The main Equation of the Ensemble Kalman Filter is given below: In this study the observation operator was considered to be 1 as both the observation and model parameter were same i.e.Snow cover fractional area.

Model (VIC) Results
VIC has been run for consecutive four years (2003)(2004)(2005)(2006) and routing has been done for the same.Coefficient of Correlation (R 2 ) as compared with the observed discharge data from Thalot for these four years is 0.65.But VIC does not consider snowmelt runoff while routing discharge, rather only considers the runoff (surface water) that has entered the local stream channel.So, in order to differentiate the rainfall runoff and snowmelt runoff, Jain et al (2015) used the method as follows.
Runoff  VIC simulated snow cover fractions are overestimated when compared with 8-day 500 m resolution MODIS satellite Snow Cover Area data which gives average snow cover fraction as about 0.42 for the same four years.According to (Sheffield et al., 2003), the VIC model uses both the subgrid vegetation tiling and elevation banding and clearly assumes snow cover if Snow Water Equivalent is present in that grid otherwise snow cover is considered to be Zero.VIC has also been run using ERA-interim once for the year 2006 to compare SCA with the original simulated SCA as well as after D.A. R 2 between SCA from ERA-interim and original MODIS SCA is 0.71 and RMSE of +0.13 (13 % SCA).Though the trend is better than that simulated by D.I and EnKF technique, still there is a scope for Data Assimilation.

Data Assimilation Results
In order to improve the simulated snow water equivalent,snow depth and snow melt and river flow , MODIS snow cover fraction has been directly inserted into the model state file and also has been inserted using the EnKF technique as explained in Eq 1,2,3.8-day 500 m resolution MODIS SCA data has been resampled into 1 km Cloud cover has been removed and SCA has been modified using the snow line elevation.Although VIC is able to simulate all the water budget and energy balance components but inaccuracy and uncertainty of input meteorological forcing data directly affects its hydrological and snow outputs such as runoff, base flow, snow depth, SCA and SWE.Although SWE has been updated well for the snow free areas (Figure 4.14, 4.15) after D.I but inaccuracy still lies in snow fed areas where SWE should directly vary with the elevation.At some days the model detects more snow in the lower elevation than the higher elevated areas.These are mainly due to the meteorological datasets used.So, VIC global Meteorological Forcing data should be used with caution in head water river basins of Himalayas, where snow and glacier dominates the overall hydrology of basin.(Thakur et al.,2016).Simulation has been done with ERA-interim for 2006 to avoid the over estimation of SCA.But there is also a need for Data Assimilation as at some places SCA is overestimated and at some, underestimated (Figure 4.11).Another area which needs attention is the Data Assimilation Techniques.May be in case of proper meteorological datasets used, EnKF may become an effective feasible solution for improving snow properties prediction.There are also many advanced techniques available such as Extended Kalman Filter, 3D/4D Variance etc. which should be looked upon too.Apart from SCA validation, Snow Depth, SWE shall also be validated.

Figure
Figure 3.1.Overall Methodology Figure 2.2.Beas Aspect Map X a = updated estimate of the analyzed state.(N*1) X b = background Model forecast (N*1) I=1 to N, where N is the number of ensembles Y = Observations (Y*1) H = Observation operator that converts the states in the model into observation space K = Kalman gain that weighs the effect of the observations to the state update R= Observation error covariance with dimension (P*P) P b= Background error covariance (P*P)

Figure 4
Figure 4.2.Snow Water Equivalent (2003-2006) (VIC) It was needed to validate the VIC based snow cover with MODIS SCA.We can clearly see that VIC predicts full snow till the month of August whereas snow starts melting from the middle of May in case of MODIS.R 2 measured is 0.61 with Root Mean Squared Error (RMSE) of 0.25 for the year 2006.As the melting of snow is directly proportional to the temperature, VIC global meteorological forcing were validated with ERA Interim data (Balsamo et al., 2012) and

Figure
Figure 4.3.Snow Melt (2003-2006) (VIC) Figure 4.10 shows the MODIS SCA for the same period which has been assimilated in the model to mark the changes.The impact of assimilating MODIS data on the SCE of the basin for 2006 with the D.I technique has been examined first.The snow cover starts declining from the month of May which has updated the SWE in a consistent manner if compared with MODIS SCA.The Figures below shows comparison of the snow components before and after Direct Insertion.As can be seen from Figure 4.8 maximum snow occurs in the end of January to the middle of February hence two peaks of SWE can also be seen at that time.But in the original simulation a shift in the SWE peak is noted which is in the month of April.Enkf propagates through the errors of both model and Satellite observation so during the start of the simulation period, SWE values (300 mm) seems to lie in between original simulation (200 mm) and that of D.I (480 mm).This trend continues till the month of march and after that it becomes less than that of simulated from D.I.The snow melt when compared with the snow cover curve as shown in Figure 4.10 shows that the melting of snow well corresponds to the snow cover after both the assimilations (DI & Enkf) compared to the original simulation.It has been well observed that from the end of February to march, there is a sudden increase in Snow melt after assimilation which is quite satisfactory if the MODIS .SCA is observed at that time as well as the snow water equivalent in Figure 4.10.When the snow cover starts declining in the middle of March,Snow melt has also decreased a bit after D.A.As in Figure 4.11 Snow Cover reaches the minimum value in the month of August whereas

from 1979-2006 with daily maximum
temperature, minimum temperature from Climate Research unit, rainfall from daily variability of NCEP and wind speed from NCEP-NCAR analysis as main inputs.
MODIS Relative daily snow cover areas for the year 2006 were derived from MOD10A2 8 day's composite snow cover data.These snow products are geocoded gridded products covering large areas, therefore boundary of study area map was used to subset them to get basin wise SCA maps which contained areas of snow ,non-snow and cloud.Based on snow line elevation, cloud cover has been removed.Then the MODIS snow cover extent observations are directly inserted into the model inputs.This method clearly assumes that the observation data are more useful than the model predicted results.To, test the effect of assimilation of fractional snow cover area in hydrological model, this approach has been implemented on Beas river basin of NWH.At first VIC model has been run for 1st 8 days of 2006, fluxes are generated for 8 days, along with that, a model state file is generated by the VIC model for a particular day i.e. 8th day.Similarly MODIS SCA is available for that particular day.Then the fractional snow cover area given by MODIS is directly replaced in the model state file.Then again the model has been run for the next eight days.The process continues so on.