Uncertainty evaluation of SWAT model for snowmelt runoff in a Himalayan watershed

Uncertainty is inherent in modeling studies. However, the quantification of uncertainties associated with a model is a challenging task. Furthermore, snowmelt estimation is a crucial part of the Soil and Water Assessment Tool (SWAT) model in watersheds where spring runoff is strongly affected by melting snow. The SWAT model for the snow dependent Kunhar basin in Himalayan watershed was calibrated (2001 - 2005) and validated (2006 - 2009) using Sequential Uncertainty Fitting Algorithm (SUFI-2). For the model uncertainty, two indices P-factor and R-factor along with frequently used objective functions R 2 , NSE, PBIAS, were taken into consideration. For calibration, multisite daily and monthly simulation results of SUFI-2 revealed that percentage of data enveloped by 95% confidence interval was 85% (monthly) to 87% (daily) at upstream calibration point and 63% (monthly) to 73% (daily) data at the downstream calibration point. Model validation by the usage of elevation bands indicated better model performance, enveloping 15 - 20% more observed data than the validation without elevation bands together with the other statistical standards. Equifinality in the model parameters was observed, and it was discovered that the model uncertainty lie inside the model parameters. It is recommended that critical model parameters correspondence with the watershed characteristics should be checked. The calibrated version of the model could be further used for the analysis and impacts of climate and land use changes on stream flows, water quality and sediment yield.


INTRODUCTION
The hydrological modeling has emerged as an important development of watershed hydrological simulation; with improvement and integration of new tools such as remote sensing/ geographic information system (RS/GIS). In current years, many such programs have been developed such as, System Hydrologic European (SHE) (Beven et al. 1980) and Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998). Traditionally, the hydrological models were calibrated and optimized with effective objective functions (e.g., Nash efficiency coefficient) and then used to assess water resources practices reasonably. But, uncertainty within model output is a major issue, predominantly when modeling results are used to set policy. The model predictions are not reached at a certain values, and be interpreted with a confidence range due to uncertainties associated with model input, structure, parameter, and output (Beven 2001;Van Griensven et al. 2008). The uncertainties in hydrological models could arise due to (1) inherent randomness of natural processes (e.g., weather, flood, precipitation variations in space and time), called natural uncertainty; (2) model structural error that reflects the inability of a model or design technique to represent precisely the system's true physical behavior, called model uncertainty; (3) model parameter value error that result from the inability to accurately quantify model parameters and inputs and also caused by the inherent inconsistency of model inputs and parameters in space and time is called parameter uncertainty; (4) data error due to measurement errors, nonhomogeneity and inconsistency of data, transcription and data handling errors, inadequate representation of data sample due to space and time limitations, is called data uncertainty. Operational uncertainties are associated with manufacture, construction deterioration, maintenance and human errors (Karamouz et al. 2012;Uniyal et al. 2015). According to Yang et al. (2008) model uncertainty is the over-simplification of natural processes and due to unknown activity in the watershed which is not considered by the model.
Snowmelt hydrology is an important component in mountainous watersheds where the stream flows are predominantly generated from melting of snow. Snowmeltrunoff modeling in a mountainous basin is perceived as difficult due to the complexity of simulation. Additionally, mountainous areas typically lack sufficient data, which leads to computational demand during simulation (Hartman et al. 1999;Fontaine et al. 2002). There are two basic approaches mostly used in snowmelt runoff modeling, i.e., temperature index approach and energy balance approach. In temperature index approach, snowmelt processes are based on temperature taken as major driving force (Tanasienko and Chumbaev 2008;Neitsch et al. 2011) while, the energy balance approach depends on the amount of energy added to the system (Valeo and Ho 2004). In energy balance approach temperature is not considered as a major driving force to explain the snowmelt processes (Marks et al. 1998(Marks et al. , 1999Zhang et al. 2007). Temperature index approach is common, simple and easy to use (Hock 2003) while, the energy balance approach is data intensive and sometimes cannot be done because of inadequate data or unwarranted detail for the work at hand (Debele et al. 2010).
There have been numerous attempts in modelling the snowmelt-runoff mechanism using SWAT model but none of them studied associated uncertainty (Grusson et al. 2015). Meng et al. (2015) developed an energy balancebased distributed snowmelt runoff model and coupled with the SWAT model. They also compared the performances of temperature index method an energy balanced methods in SWAT model. The performances of energy balance method and temperature index method was also compared in SWAT model (Debele et al. 2010). Kang and Lee (2014) studied the effects of snowmelt and temperature in the North Fork American River basin using the snowfall-snowmelt routine in SWAT. They concluded that the snowmelt -runoff model gives better results in terms of ENS, R 2 , and RMSE. Wang and Melesse (2005) studied the performance of a SWAT model on snowmelt-runoff simulation, and concluded that SWAT model shows acceptable accuracy in simulating mean discharges. Tahir et al. (2011) found that, in the high altitude basin, temperature change is one of the most sensitive input parameters in estimating snowmelt runoff. Fontaine et al. (2002) modified the original snowfall-snowmelt routine in SWAT by incorporating snow accumulation, snowmelt, areal snow coverage, and an option to input precipitation and temperature as a function of elevation bands.
In Pakistan the northern part of the country, with heterogeneity of elevation and diverse forests and glacier dy-namics present challenges as well as potential research opportunities to understand mountain hydrological processes. Mangla Dam is the first controlled structure on the River Jhelum and most of the annual Jhelum River influx stored at Mangla reservoir is derived from the mixed runoff generated from rainfall and snowmelt in the Himalayan and Pir Panjal ranges. The stored water of Mangla Reservoir is then delivered to the downstream irrigated agricultural lands through a network of barrages, canals and small watercourses. Nearly 60% of the whole catchment is covered by seasonal snow cover during the winter months (October-March). Keeping in view the relative importance of snow cover for river flow within the study area, the snowmelt runoff modeling at daily and monthly time scale was done using physicallybased semi-distributed hydrological SWAT model. Also the quantification of SWAT model uncertainty for Kunhar watershed was performed. To the best of our knowledge, this study is performed first time on SWAT model uncertainty associated in snowmelt-runoff modeling.

Study Area
SWAT model was tested on Kunhar basin with outlet points Naran and Gari-Habibullah covering the catchment area of 2442 km 2 at Gari-Habibullah. The contribution of this basin to the total annual flow of Jhelum River is 11% on an average (De Scally 1994). The major portion of this catchment is covered by snow and glaciers during the winter season. The river Jhelum originates from Verinag Spring situated in between Himalaya mountain ranges and Pir Panjal ranges in Jammu and Kashmir. Large tributaries of Jhelum River are Neelum River and Kunhar River (Fig. 1). The flow of Jhelum River enters Mangla reservoir in the Mirpur district. Kunhar watershed has two peak flow events; one occurs in June and the other in July-September. The higher June inflow is attributed to the increased quantity of snowmelt (due to temperature rise) while the July-September peak is a combination of rainfall and snowmelt. The high elevation Naran and Astore stations receive more snow than rainfall. The runoff from these sub-basins is controlled by temperature variations and is increased during summer due increase in temperature (Archer and Fowler 2008). The information in Table 1 shows the annual temperature variations along with location of weather stations and flow gauges.

SWAT Model and Snow Package
The SWAT is a comprehensive, process based semi distributed hydrological model (Arnold et al. 1998(Arnold et al. , 2012. It subdivides the whole watershed into sub-watersheds based on topography and then identifies Hydrological Response Units (HRUs) within each sub-watershed, based on land No reference! E a r l y R e l e a s e  E a r l y R e l e a s e use, soil and slope. The HRUs are non-spatially distributed assuming there is no interaction and dependency (Neitsch et al. 2011). The HRUs are used to compute a water balance based on snow, soil, shallow aquifer, and deep aquifer. Water balance computation is performed at the HRU level of spatial discretization, and contributions of each HRU are then aggregated at subbasins level. The water yield is then routed toward the reach and outlet of watershed. The snow module of SWAT consists of snowfall, snowpack, and snowmelt processes and it provides snow output in term of snowmelt. The snowpack increases with additional snowfall, but decreases with snowmelt or sublimation. The mass balance of the snowpack is computed as: The snowfall temperature (SFTMP) is the threshold parameter which defined the precipitation as snowfall or rainfall. If mean daily air temperature is lower than SFTMP, the precipitation is considered as snowfall and is added to snowpack. The snowpack will not melt until the snowpack temperature exceeds the snowmelt base temperature (SMT-MP) threshold value. The snowpack temperature of the current day is calculated as: where T snowi is the snowpack temperature on the current day i (°C), T snowi 1 -is the snowpack temperature on the previous day i -1 (°C), TIMP snow is the snow temperature lag factor, Tavi is the mean air temperature on the current day i (°C).
When the snowpack temperature reaches a higher than SMTMP, the snowpack start melting. The amount of snowmelt is estimated on the basis of following relationship. The snowmelt rate is based on two boundaries, i.e., the rate is maximum on 21 st June and minimum on 21 st December.
Where SNO mlti is the amount of snowmelt on day i (mm H 2 O), b mlti is the melt factor for day i (mm H 2 O °C -1 d -1 ) as defined in Eq. (4), SNO covi is the fraction of HRU area covered by snow as defined in Eq. (5), T maxi is the maximum air temperature (°C) on current day, b mlt, max is the melt factor for 21 st June (mm H 2 O °C -1 d -1 ), b mlt, min is the melt factor for 21 st December (mm H 2 O °C -1 d -1 ), SNOCOVMX is the threshold depth of snow at 100% coverage (mm H 2 O), and cov 1 and cov 2 are coefficients that define the shape of the curve. The more description of equations can be found in Neitsch et al. (2011) and Arnold et al. (2012).

Model Data Input
The input data files necessary for SWAT modeling and analysis in this study consist of geospatial data and temporal data. The geospatial data includes topography (DEM), landuse and soil maps while the temporal data consist of climate and streamflow data. The data were collected from different sources/agencies (Table 2). A 30 m × 30 m digital elevation model (DEM) with vertical accuracy of ±20 m at 95% confidence was obtained from USGS National Elevation Dataset (Gesch et al. 2002;Gesch 2007) as shown in Fig. 2a. The DEM was used to delineate the sub-basins and drainage network. The slope classification map was also prepared using the DEM data. The landuse data at a spatial resolution of 300 m was obtained from European Space Agency (ESA) database ( Fig. 2b). Landuse classification was done into eight major classes: irrigated croplands (14.06%), urban areas (10.80%), forests deciduous (0.58%), forests evergreen (6.38%), forest mixed (33.14%), rangeland (2.88%), grasslands (19.4%), and water bodies (12.75%) ( Table 3). These landuse classes were used to define the corresponding landuse type in the SWAT model. The soil map was obtained from the global IPCC soil classes (5 km resolution) of United Nation Food and Agriculture Organization (FAO) regional scale soil database. In the study area 4 soil classes were extracted (Fig. 2c) having the physical attributes of available water capacity, texture, saturated conductivity, bulk density, organic carbon and soil albedo.
The daily weather data including precipitation, maximum and minimum air temperature, solar radiation, and relative humidity of four climate stations (Balakot, Muzaffarabad, Naran, and Astore) as shown in Fig. 1

SWAT Model Setup
ArcSWAT 2012 with an interface with ArcGIS was used to delineate subwatersheds, HRUs and to setup the model in this study. Based on topography, the watershed is subdivided into 20 subwatersheds ( Fig. 3) with drainage areas extending from 9.41 -326.8 km 2 . Upon the overlapping of land use, soil and slope maps these 20 subwatersheds were subdivided in 421 HRUs. The HRUs are defined as homogeneous spatial units, characterized by similar geomorphologic and hydrological properties. Number of HRUs must be adequate to reason for the different hydrologic conditions in the watershed, and limited enough to reduce the data input requirements and improve modeling efficiency (Qiu and Wang 2014). For defining HRUs, areal thresholds for landuse, soil and slope maps were set to 4, 6, and 8%, respectively. Any class value less than the threshold was reassigned to predominant land use, soil or slope class in that subwatersheds. Soil Conservation Service (SCS) curve number method (Arnold et al. 1998) was utilized to calculate surface runoff. Muskingum method was used for the estimation of channel routing. Penman-Monteith method coupled with Simplified Pant Growth Model was used for the calculation of Potential evapotranspiration (PET). The snow package along with elevation band method was applied in SWAT. The elevation band method is helpful in consideration of orographic effects on precipitation and temperature in mountainous areas. SWAT allow maximum 10 elevation bands for each subbasin. In this study, five elevation bands were applied within 20 subbasins (Fig. 3). Snow cover and snowmelt simulations were done separately for each elevation band. The SWAT band selection was done based upon Fontaine et al. (2002) and Pradhanang et al. (2011) studies, in which they concluded that using three or five elevation bands improved simulation. The bands were setup in all the snow dominated subbasins keeping equal vertical distance from the mean elevation of the centroid of the subbasins. The missing values on daily precipitation, temperatures, along with solar radiation, wind speed, and relative humidity has been simulated by the weather generator, incorporated in the SWAT model package (Neitsch et al. 2011).

Model Sensitivity Analysis, Calibration, and Validation
The SWAT model was simulated for time period 2001 to 2009 in which first four years were set as a warm-up period in order to establish appropriate initial parameter values for the 20 subbasins. The calibration of hydrological parameters along with snow parameters was performed using multisite daily, and monthly stream flow data during 2001 to 2005. The model validation was setup at a daily and monthly time scale with and without elevation bands for 2006 to 2009. First, the sensitivity analysis of the snow and hydrological parameters was performed using Latin-hypercubeone-factor-at-a-time (LH-OAT) method (Van Griensven et al. 2006) in SWAT-CUP (Abbaspour et al. 2015). The aim of this process was to identify the parameters that considerably affect the stream flow (Winchell et al. 2009). Among the top ranked sensitive parameters, the parameters effecting the stream flow were selected for the automatic calibration. The auto-calibration technique was undertaken using SWAT-CUP and its Sequential Uncertainty Fitting SUFI-2 algorithm (Abbaspour 2011) for calibration of parameters. The SWAT-CUP is an external software tool which provides the advantage to SWAT users to do multiple iterations by using various algorithms until the best fit results obtained (Arnold et al. 2012). The SUFI-2 algorithm was selected because it is a semi-automatic parameter optimization algorithm allowing users to perform good calibration iteratively in a limited number of iterations. Moreover, this method recently has been increasingly used due to having a high computational efficiency (Rostamian et al. 2008;Yang et al. 2008;Abbaspour et al. 2015;Narsimlu et al. 2015).

Model Statistical Evaluation Criteria
Several researchers have proposed standard statistical model evaluation criteria. Santhi et al. (2001)   E a r l y R e l e a s e R 2 and NSE have been used by many researchers as a measures for the performance evaluation of SWAT (Santhi et al. 2001;Rahman et al. 2013;Magali and Daniel 2014;Rahayuningtyas et al. 2014). Another factor followed in this study was Percent Bias (PBIAS) (Moriasi et al. 2007 Where O i = observed value; O = mean observed value; S i = simulated value; and S = mean simulated value (Santhi et al. 2001).

Uncertainty Prediction Criteria
The performance of the calibrated models should be evaluated before use. The model predicted results doesn't reached to a certain value due to uncertainties associated with model input, model structure, parameters, and output. The results are represented with a confidence range (Gupta et al. 1999;Beven 2001). So in this study, the strength of a calibration and uncertainty analysis is quantified by using two additional statistical uncertainty prediction criteria referred as the P-factor and R-factor. The P-factor is defined as the percentage of measured data enveloped by the 95PPU (95% prediction uncertainity), while the R-factor is the average thickness of the 95PPU envelop divided by the standard deviation of the measured data. Theoretically, the values for the P-factor and R-factor range between 0 and 100%, and between 0 and infinity, respectively. The value of P-factor equal to 1 (100%) and that of R-factor is close No reference! E a r l y R e l e a s e to zero indicate that the simulated results exactly matching with the observed values (Abbaspour 2011). The P-factor and R-factor are computed as follows.

P factor N ny ti
Where ny ti is the number of measured values enveloped by the 95PPU, dx is the mean of the x variables, x v is the standard deviation of the measured variable x, N the total number of measured values.

Model Performance and Calibration
The sensitivity analysis identified fifteen most influential parameters in which ten were hydrological parameters while the five were snow parameters. The selected sensitive parameters with description are given as follow: Soil Conservation Service (SCS) curve number II (CN2), temperature lapse rate (TLAPS), base-flow alpha factor (ALPHA_ BF), available water capacity of the soil layer (SOL_AWC), Snowpack temperature lag factor (TIMP), snowmelt base temperature (SMTMP), maximum melt factor (SMFMX), threshold depth of water in the shallow aquifer (GWQMN), minimum melt factor (SMFMN), manning "n" for overland flow (CH_N2), average slope length (SLSUBBSN), average slope steepness (HRU_SLP), groundwater delay (GW_DE-LAY), plant uptake compensation factor (EPCO), surface lag coefficient (SURLAG). The sensitivity ranking of these parameters from the order of high to low along with corresponding range is given in Table 4.
It was observed that among the top ten sensitive parameters, there are five snow parameters. Temperature lapse rate (TLAPS) ranks second most sensitive parameter, for the reason that it is related directly with the snow or glacier melt. TIMP is also among the sensitive parameters because it considers the snowmelt for the preceding day and it indicates the impact of daily temperature variations on the daily snowpack temperature. SMTMP is the threshold temperature above which the snowmelt process starts and ranked as third most sensitive parameter among snow parameters. It controls the peaks and shape of the simulated hydrograph. SMFMX and TIMP are two interrelated parameters and control the snowmelt estimation. During the modifications of parameter values, it was observed that TIMP parameter has great impact on the model efficiency as compared to SMFMX parameter. So the TIMP is the second most sensitive parameter while SMFMX is on 4 th position among the snow parameters.
The other most dominant hydrological parameters are CN2 and ALPHA_BF. CN2 is rainfall runoff conversion coefficient and controls the runoff generation. It depends on soil type, land use type, and slope. It varies from HRU to HRU, reach to reach and sub-basin to sub-basin. ALPHA_BF and E a r l y R e l e a s e GWQMN which are ground water parameters, control the interchange between the stream flow and the ground water evaporation in the unsaturated zone (Smedema and Rycroft 1983). Moreover, by modifying the SLSUBBSN, HRU_SLP, and OV_N parameters, the lag effect on simulated flows were reduced and reshaped the simulated hydrographs. During the model sensitivity analysis and calibration, the baseline values allocated to each spatially varying parameter were changed by multiplying the baseline value by adding the sampled values to the baseline value shown in (Table 4).
The calibrated and model predictive performances for Kunhar River on daily and monthly flows are summarized in Table 5 for all calibration and uncertainty analysis method. The model parameter optimization was performed on the bases of five objective functions which are: (1) NSE (2) R 2 , (3) PBIAS, (4) P-factor, and (5) R-factor. The comparison between the observed and simulated stream flow indicated that there is a good agreement between the observed and simulated discharge which was verified by higher values of coefficient of determination R 2 and NSE. Overall, the performance of the model for both calibration points was good for daily (NSE = 0.68 and NSE = 0.74) and monthly (NSE = 0.73 and NSE = 0.78) time scale (Table 5). The model performance of the daily time scale was weak as compared to the monthly. Many researchers have reported SWAT lower performance for the simulation of daily stream flows (Yuan et al. 2001;Van Liew and Garbrecht 2003;Saleh and Du 2004;Fernandez et al. 2005;Van Liew et al. 2007). PBIAS value varied from 3.2 -4.1% for the calibration at Naran and 9.1 -12.7% at monthly time scale. The possible reason may be the weakness of SWAT model in modeling watershed flows which are predominately generated from snowmelt (Fontaine et al. 2002). A close analysis of five years calibration period on a daily scale indicated that model underestimated runoff for two years at Naran and four years at Gari-Habibullah while overestimated runoff for one year at Naran and three years at Gari-Habibullah as shown in Fig. 4.

Model Validation and Influence of Elevations Bands
The results over the validation period showed similar performance as the calibration period. Model simulation for the validation period captured the overall dynamics of the flow with greater accuracy at both calibration points. The model validation was performance by two techniques, with and without the use of elevation bands. The validation which was done without inclusion of elevation band is named as validation technique-1 (VT-1) while with band inclusion is VT-2. Results of SUFI-2 during the validation periods at two flow stations are shown in Fig. 5, the shaded area contains the uncertainty from different sources.

Validation Technique-1 (VT-1)
The validation technique-1 (VT-1) was performed without inclusion of elevation bands. The results of Naran gauging station for daily time step (Fig. 5a) showed an acceptable value of R 2 (0.71) and NSE (0.63). However, the flows were underestimated (PBIAS = 12.3%) higher than the calibration period. The uncertainty at this station increased significantly compared to the calibration period. Only 64% of the observations were enveloped by corresponding 95% confidence interval. The R-factor at this station was 0.70 (Table 6). For monthly time step, the model performance improved in terms of selected objective functions. The values of R 2 and NSE were greater than 0.80, which are higher than the daily flows calibration. Moreover, the P-factor enveloped 71% observed stream flows with R-factor of 0.90. The average tendency of the data improved from 12.3 to -2.6% as compared to daily flows calibration.
Gari-Habibullah flow station (Fig. 5b) for daily flows validation demonstrated R 2 and NSE > 0.65, which are acceptable as per Santhi et al. (2001). 95% prediction uncertainty (95PPU) band captured 72% observed data and R-factor equaled 0.89. From monthly time step, the model indicated large uncertainty capturing only 63% of observed data with R-factor of 0.74. However, PBIAS improved from 11.03 -8.6% (Table 6). Overall results pointed out that model improved performance for monthly time step in terms of R 2 , NSE, and PBIAS.

Validation Technique-2 (VT-2)
In validation technique-2 (VT-2), the model was validated by the inclusion of elevation bands. Model performance

Naran
Gari-Habibullah No reference! E a r l y R e l e a s e E a r l y R e l e a s e measures were higher in the validation (VT-2) as compared to VT-1. At Naran gauging station for daily time scale 74% of observed data have been captured by corresponding 95PPU band while R-factor equaled 0.78. Compared with the VT-1, P-factor improved from 64 -74%, with values of R 2 and NSE > 0.70 (Table 6). This improvement in model performance indicates the significance of elevation bands for the snowmelt modeling. Monthly flow simulation at this station indicated similar model performance as of daily. At Gari-Habibullah flow station for daily time scale, 89% of flow data was falling under the 95% confidence interval, which was higher than the P-factor in VT-1 and signifying less uncertainty. The values of R 2 and NSE were also higher, as compared to the validation VT-1 (Table 6). For monthly flow validation, the P-factor was improved, which indicated less degree model output uncertainty. The PBIAS also decreased from 8.6 -6.3%. The simulated flows matched the observation flows quite well with the optimum objective function values. On an average 10 -12% more data for the daily flow stimulation was captured by 95PPU band and 15 -25% for the monthly flow validation by using elevation bands. Moreover, some of the extreme flow events were also simulated in a better way. For example, at Gari-Habibullah flow station (Fig. 5d), the high flows were simulated better as compared to the simulation without elevation bands.

Uncertainty Analysis
Uncertainty in SWAT model was estimated through the SUFI-2 method. SUFI-2 was started with wide but meaningful ranges of sensitive parameters so that the measured data initially fall within the 95PPU; this uncertainty is then reduced by several iterations. The 95PPU was used to quantify all types of uncertainties. It combined outcome of the uncertainties in the input data, hydrological model, and the parameters. These uncertainty sources were not estimated independently but ascribed as total model uncertainty in model parameters. Normally, P-factor and R-factor were used to assess the uncertainty. P-factor is the percentage of observed data enveloped by the 95PPU modeling result. R-factor is the thickness of the 95PPU envelop. The reasonable values for P-factor and R-factor are > 70% and around 1, respectively for better result of discharge simulation (Abbaspour 2011). It was assessed that at Naran flow station the P-factor and R-factor were 0.87 and 1.08 at daily time steps. The P-factor indicated that 87% of the observed data were captured by 95PPU. At Gari-Habibullah flow station the P-factor and R-factor were observed to be 0.73 and 0.97. This reduction indicated the higher uncertainties as the stream flows moved towards downstream. The 95PPU captured 14% less data than the upstream calibration point which indicated large uncertainty prediction in driving input variables such as rainfall (Table 5).
For monthly time scale simulation at Naran station the P-factor was 0.85 and R-factor = 0.91. The 95PPU captured 2% less data than at the daily time scale. Also the R-factor decreased from 1.08 to 0.91. While at Gari-Habibullah the P-factor and R-factor were 0.63 and 0.78, respectively. This large reduction in P-factor and R-factor at the monthly time simulation is indicative of large uncertainties (Table 5). It was predicated that the model simulation on monthly time scale was good at both stations based upon R 2 and NSE, while the P-factor and R-factor values indicated that SWAT model exhibited a certain degree of uncertainty on monthly time step simulation. Close analysis of 95PPU band concealed that some of the low flows were not covered properly, this may be due to the limitation of SWAT model for simulating groundwater flow (Rostamian et al. 2008). Many researchers also reported difficulty in achieving good SWAT simulations for low flow conditions (Van Liew and Garbrecht 2003;Sudheer et al. 2007). This could be attributed to the inability of the runoff curve number (CN) to adequately account for hydrologic abstractions for various antecedent soil moisture conditions. Some of the high flow events at both gauging stations (Fig. 4b) were also not attributed by 95PPU which indicated the inability of SWAT model for extreme flow simulations.
It was observed that the model uncertainties varied at both flow stations. This may be due to inadequate distribution of weather stations inside the catchment. Since it is a mountainous watershed with an elevation difference of 4403 m, and covered with only 4 weather stations of which 2 station were inside the watershed. This spatial distribution variability of climate data is likely to be captured insufficiently. As only one weather station is located in Naran at E a r l y R e l e a s e upstream part of the catchment and the other one station at Astore which is 15 km outside the catchment boundary. At downstream, only one weather station is also in the catchment while the other one at Muzaffarabad, 8 km outside the catchment. This spatial distribution could be the cause of error in input weather data within this catchment because according to SWAT module, the weather condition in a subbasins are determined by the nearest weather station and rain gauge station. Thus, the input data may not represent the real weather conditions in this study catchment.

Parameters Equifinality and Identification
The dotty plots (Fig. 6) mapped the model parameter (xaxis) and NSE (y-axis) to illustrate the relative sensitivity of each parameter associted with NSE and also demonstrate the distribution of sampling points (Wagener and Kollat 2007). By observing the dotty plot Fig. 6, it was evident that the main sources of streamflow uncertanity for snowmelt modling were hydrological parameters (CN2.mgt, ALPHA_BF.gw, SOL_AWC.sol, GWQMN.bsn, and CH_N2.rte) and snow parameters (TLAPS.sub, TIMP.bsn, SMTMP.bsn, SMFMX. bsn, and SMFMN.bsn). All the parameters are concentrated within highest values of NSE (0.6 -0.8). Among the all the parameters, Ground water governing parameter ALPHA_BF and snowmelt parameter SMFMX were the most identifiable parameters for the study watershed. The scattered and/ or haphazard points on dotty plot indicated low sensitivity of parameters while if the points do follow a trend, the sensivity is higher. Most of the parameters were showing similar pattern, generating almost similar model performance. This is called parameter equifinality. However, it should be noted Fig. 6. Dotty plots depicting the variability of identifiability of SWAT parameters. E a r l y R e l e a s e that non-identifiability of a parameter does not mean that the model is not sensitive to these parameters (Shen et al. 2012). The results of dotty plots indicated that the sensitivity of SWAT parameters did not vary between the different parameter values. Even though many of the sensitive parameters were affecting the stream flow, but only a few sensitive parameters were identifiable. This equifinality is an indication of no unique parameter were estimated. This may be the indication that uncertainty in model output lies in the parameters, similar results were reported by Demaria et al. (2007) and Narsimlu et al. (2015). This may occured due to following reasons: (1) correlations between parameters, (2) spatial and temporal scales of model residual, (3) sensitivity or insensitivity of parameters (Wagener and Kollat 2007). Consequently, to deal with non-identifiable parameters in SWAT calibration, one must be careful because it may lead to equifinality of the parameter values. Under these circumstances, the user must check, if the final calibrated model values correspond to watershed characteristics and its underlying processes. To obtain more credible results of watershed management, a monitoring task must be done for the important parameters.

CONCLUSION
This study assessed the performance of SWAT model in the high altitude watershed where runoff is a subtle mixture of both rainfall and snowmelt processes. The approach starts with the initial range of parameters based on the understanding of the physics of the watershed. Then after the iterations, the value of parameters altered with respect to best fit objective function. The model was successfully calibrated (2001 -2005) and validated (2006 -2009) for the two gauging stations at daily and monthly time scale. The statistical evaluation of model results at daily time scale indicated good results verified by the NSE > 0.72 and R 2 > 0.68, whereas at monthly time scale both R 2 and NSE were high significantly. The model validation was performed with and without the use of elevation bands. Model results in terms of selected objective functions confirmed that model results improved significantly when elevation bands were used. The SUFI-2 algorithm was used to quantify the associated uncertainty. For calibration, P-factor and R-factor computed by SUFI-2 indicated good results by enveloping more than 73% of the observed data at daily time scale. While on a monthly time scale more than 63% observed data was enveloped by 95% confidence interval, indicating higher uncertainties. For validation, overall the P-factor enveloped 15 -25% more observed data when elevation bands were used. During the model parameter identification, the equifinality in the parameters was observed, generating almost similar model performance for different set of parameters. Greater parameter sensitivity does not mean that the parameter is also identifiable. The study suggested that the final cali-brated model should be in correspondence to the watershed characteristics and its underlying processes, and implementation of monitoring task for the important non-identifiable model parameters. This study showed that SWAT model can produce reliable estimation of the different components of the hydrological cycle in high mountain watersheds. The calibrated model can be used for further analysis of the effect of climate and land use changes, as well as to investigate the effect of different management scenarios on stream flows, water quality and sediment yield.