Hydrological Modeling of Highly Glacierized Basins ( Andes , Alps , and Central Asia )

The Soil and Water Assessment Tool (SWAT) was used to simulate five glacierized river basins that are global in coverage and vary in climate. The river basins included the Narayani (Nepal), Vakhsh (Central Asia), Rhone (Switzerland), Mendoza (Central Andes, Argentina), and Central Dry Andes (Chile), with a total area of 85,000 km2. A modified SWAT snow algorithm was applied in order to consider spatial variation of associated snowmelt/accumulation by elevation band across each subbasin. In previous studies, melt rates varied as a function of elevation because of an air temperature gradient while the snow parameters were constant throughout the entire basin. A major improvement of the new snow algorithm is the separation of the glaciers from seasonal snow based on their characteristics. Two SWAT snow algorithms were evaluated in simulation of monthly runoff from the glaciered watersheds: 1) the snow parameters are lumped (constant throughout the entire basin) and 2) the snow parameters are spatially variable based on elevation bands of a subbasin (modified snow algorithm). Applying the distributed SWAT snow algorithm improved the model performance in simulation of monthly runoff with snow-glacial regime, so that mean RSR decreased to 0.49 from 0.55 and NSE increased to 0.75 from 0.69. Improvement of model performance was negligible in simulation of monthly runoff from the basins with a monsoon runoff regime.


Introduction
The two basic snowmelt approaches generally used in hydrologic modeling are energy balance models and temperature-index models [1].Temperature-index models are widely used in hydrologic studies because of the models' performance and simplicity as well as the wide availability of temperature data [2][3][4].
In temperature-index based runoff models such as SWAT, melt rates only vary as a function of elevation, because of the air temperature gradient [5].To overcome this drawback, a modified snow process was incorporated into SWAT in order to consider spatial variation of snowmelt/accumulation parameters by elevation band across subbasins.Previous studies using SWAT kept these parameters constant for an entire basin [6][7][8][9][10].While this method was successful in simulation of snowmelt flow, simulation of runoff from glaciered watersheds demands a distributed model for distinguishing seasonal snow from glacier.The new approach allows for the separation of seasonal snowmelt from glaciers' melt based on the spatial variability of associated melt parameters.The SWAT model has been applied worldwide, and its hydrologic components have been successfully tested in cases where streamflow was predominantly generated from rainfall events [11][12][13][14][15].The model less frequently has been applied in mountainous watersheds, and a few recent studies have been conducted to test and improve SWAT's snow hydrology component.Fontaine et al. [16] incorporated the elevation bands method with SWAT's original snowmelt algorithm (temperature-index model), which improved the Nash-Sutcliffe coefficient of monthly runoff simulation from −0.70 to 0.86 in the 4999-km 2 Rocky Mountain Basin in Wyoming.Debele and Srinivasan [4] incorporated a modified version of SNOW17 [17] into SWAT and compared its performance with the temperature-index model in three watersheds (ranging from 22.28 to 7106.82 km 2 ), the results of which showed that the temperature-index model performed better than the SNOW17 model.Debele et al. [18] incorporated the distributed process-based energy budget SNOWEB in the pixel and elevation band scales into SWAT and compared its performance with the temperature-index model.In this method, it was assumed that solar radiation varies not only with latitude and altitude of subbasins, which is applied in the current version of SWAT, but also with land surface inclinations (aspect and slope).The temperature-index based snowmelt computation method had overall Nash-Sutcliffe model efficiency coefficient ranging from 0.49 to 0.73 in simulation of monthly streamflow while the energy budget based approach had Nash-Sutcliffe efficiency coefficient ranging from 0.33 to 0.59.Zhang et al. [8] applied SNOW17 in SWAT at the pixel scale.The SWAT model with temperature-index plus elevation bands performed as well as the SWAT model with SNOW17.Rahman et al. [19] applied the modified temperature-index method [5] to simulate the runoff from fully and partially glaciated subbasins of the Rhone River Basin.The applied method showed improvement in simulated streamflow from the main outlet of the basin even though the snowmelt parameters were lumped across the basin.
Distributed, process-based energy budget models have been tested in SWAT, but we focused on another component of frequently used glacier/snowmelt models.A major difference between SWAT melt processes and the melt routine of hydroglacial models arises from the associated melt parameter distribution (melt factor).In previous studies using SWAT, associated snowmelt parameters were considered to be uniform within each subbasin, but in the hydroglacial models, melt factors are spatially variable in pixel or band scales.A common approach is to assign two different melt factors to ice and snow.This enables the user to treat seasonal snow and glaciers differently; melt factors are generally reported to be higher for glaciers and lower for snow [20].The applied approach in this study separates seasonal snow from glaciers by setting the accumulation/melt parameters, including maximum and minimum melt factors, melt lag factor, melt temperature and snow fall temperature for the elevation bands of each subbasin.The details of the proposed approach are available at Omani et al. [21] for three benchmark glaciers in Asia and Europe.The objective of this study was to extend the applied method to macro-scale river basins that are global in coverage and vary in climatic condition.Two snow process algorithms (lumped and distributed) were examined for their ability to simulate flow in five river basins that provide global assessment and feature contrasts in climatic conditions.

Study Watersheds
This study focuses on areas where climate variability has had a strong impact on highly glaciated areas.Five river basins with a total area of 85,000 km 2 , in the northern and southern hemispheres, for which sufficient hydrologic information is available were selected for this study.These river basins have different spatial scales and climatic situations, from extreme maritime to extreme continental climates.They include: Vakhsh in Tajikistan, Narayani in Nepal, Mendoza in Argentina, five individual glaciated watersheds in the central dry Andes of Chile, and Upper Rhone in Switzerland (Figure 1).

Hydrologic Data
For this study, a 90-m resolution Global SRTM Digital Elevation Model contributed by USGS/NASA was used to describe topographic conditions such as slope and slope length, to create flow direction and accumulated flow, and to delineate watershed boundaries and channel networks [23].Land use information was adopted from USGS Global Land-Cover Characteristics (GLCC).GLCC was developed using satellite data collected from 1992 to 1993 with 1 km spatial resolution for the entire globe [24].Land use information can be found in Omani et al. (2016b).The soil properties dataset with 5 arc-minutes resolution (Ver.1.2) was taken from ISRIC-WRS [25].The major soil type in the upper areas of the river basins is Lithosols, which are very shallow and occur mainly on steep slopes, often with exposed rock debris.Weather data were obtained from the Climate Forecast System Reanalysis (CFSR) global meteorological dataset [26].Global CFSR data are available going back to 1979 at a 38-km resolution.Daily discharge records for model calibration were collected from local hydrologic administrators and online databases (Table 1).

Model Simulation
This section focuses on glacier modelling in SWAT.More information about SWAT snow processes and theoretical details of the applied approach in this study can be found in Omani et al. [21].

Subbasin Delineation and HRU Definition
The watershed boundary of the river basins was determined using the SWAT interface automatic delineation procedure (Figure 2).Based on the thresholds selected, there were a total of 675 subbasins in the five river basins.SWAT further divides each subbasin into smaller areas called Hydrological Response Units (HRUs).There were a total of 9878 HRUs in the five river basins.

Hydrologic Data
For this study, a 90-m resolution Global SRTM Digital Elevation Model contributed by USGS/ NASA was used to describe topographic conditions such as slope and slope length, to create flow direction and accumulated flow, and to delineate watershed boundaries and channel networks [23].Land use information was adopted from USGS Global Land-Cover Characteristics (GLCC).GLCC was developed using satellite data collected from 1992 to 1993 with 1 km spatial resolution for the entire globe [24].Land use information can be found in Omani et al. (2016b).The soil properties dataset with 5 arc-minutes resolution (Ver.1.2) was taken from ISRIC-WRS [25].The major soil type in the upper areas of the river basins is Lithosols, which are very shallow and occur mainly on steep slopes, often with exposed rock debris.Weather data were obtained from the Climate Forecast System Reanalysis (CFSR) global meteorological dataset [26].Global CFSR data are available going back to 1979 at a 38-km resolution.Daily discharge records for model calibration were collected from local hydrologic administrators and online databases (Table 1).

Model Simulation
This section focuses on glacier modelling in SWAT.More information about SWAT snow processes and theoretical details of the applied approach in this study can be found in Omani et al. [21].

Subbasin Delineation and HRU Definition
The watershed boundary of the river basins was determined using the SWAT interface automatic delineation procedure (Figure 2).Based on the thresholds selected, there were a total of 675 subbasins in the five river basins.SWAT further divides each subbasin into smaller areas called Hydrological Response Units (HRUs).There were a total of 9878 HRUs in the five river basins.

Initial Glacier Storage
The equilibrium line altitude (ELA) is a theoretical line where climatic mass balance is zero on average over a number of years and determines the boundary of the accumulation zone and ablation zone on a glacier.The ELA of a glacier is sometimes denoted ELA0 (the balanced-budget ELA).The

Initial Glacier Storage
The equilibrium line altitude (ELA) is a theoretical line where climatic mass balance is zero on average over a number of years and determines the boundary of the accumulation zone and ablation zone on a glacier.The ELA of a glacier is sometimes denoted ELA 0 (the balanced-budget ELA).The ELA is determined by climate and the aspect of the glacier.It is not influenced by glacier dynamics, extent and hypsometry, and thus reveals a largely unfiltered climatic signal [27].Therefore, it can be a representation of the lowest boundary of climatic glacierization [28].
In this study, it was assumed that the ELA 0 represented the glacier boundary.The physical boundaries of glaciers were then corrected based on glacier inventory data, GLIMS glacier outlines and MODIS products for the glacier-free areas at elevations higher than the ELA 0 , and the areas climatically suited for glacier formation at altitudes lower than ELA 0 .Elevation bands below the ELA 0 were considered to be seasonal snow cover regardless of glacial tongues extending into lower elevations.Debris-covered areas of glaciers, accumulated wind-blown snow, and small isolated glaciers with an area of 0.1 km 2 or less were ignored when estimating the total glaciated area of the river basins.ELA 0 values were obtained from the literature and are available in Table 2. Glacier covered areas were extracted from the Global Land Ice Measurements from Space (GLIMS) [37,38] dataset and the World Glacier Inventory (WGI).If the GLIMS dataset did not completely cover the entire study area it was complimented with glaciated area by the Moderate Resolution Imaging Spectroradiometer (MODIS) maximum snow extent product (MOD10A2) with eight-day composites at 500-m resolution.The minimum snow cover area at the end of the melting season, from February 2002 to 2010, was considered to be the glacier/permanent snow cover area (Figure 3).The percentage of the glacier covered area of each river basin by each model, MODIS and GLIMS, is presented in Table 2.It is assumed that GLIMS glacier outlines represent the glacial extent at the end of the reference period as a starting point for the future simulations of glacial extent.The mean ice thickness in each elevation band was calculated by averaging the measured ice thickness values from the WGI.When no data were available for high elevations, ice thickness was assumed to be 1000 m water equivalent (m.w.e) (SNOEB = 999,999 mm H 2 O).Figures A1 and A2 in Appendix A show the measurement locations of glacier depth and area from the WGI data inventory for the river basins.
ELA is determined by climate and the aspect of the glacier.It is not influenced by glacier dynamics, extent and hypsometry, and thus reveals a largely unfiltered climatic signal [27].Therefore, it can be a representation of the lowest boundary of climatic glacierization [28].
In this study, it was assumed that the ELA0 represented the glacier boundary.The physical boundaries of glaciers were then corrected based on glacier inventory data, GLIMS glacier outlines and MODIS products for the glacier-free areas at elevations higher than the ELA0, and the areas climatically suited for glacier formation at altitudes lower than ELA0.Elevation bands below the ELA0 were considered to be seasonal snow cover regardless of glacial tongues extending into lower elevations.Debris-covered areas of glaciers, accumulated wind-blown snow, and small isolated glaciers with an area of 0.1 km 2 or less were ignored when estimating the total glaciated area of the river basins.ELA0 values were obtained from the literature and are available in Table 2. Glacier covered areas were extracted from the Global Land Ice Measurements from Space (GLIMS) [37,38] dataset and the World Glacier Inventory (WGI).If the GLIMS dataset did not completely cover the entire study area it was complimented with glaciated area by the Moderate Resolution Imaging Spectroradiometer (MODIS) maximum snow extent product (MOD10A2) with eight-day composites at 500-m resolution.The minimum snow cover area at the end of the melting season, from February 2002 to 2010, was considered to be the glacier/permanent snow cover area (Figure 3).The percentage of the glacier covered area of each river basin by each model, MODIS and GLIMS, is presented in Table 2.It is assumed that GLIMS glacier outlines represent the glacial extent at the end of the reference period as a starting point for the future simulations of glacial extent.The mean ice thickness in each elevation band was calculated by averaging the measured ice thickness values from the WGI.When no data were available for high elevations, ice thickness was assumed to be 1000 m water equivalent (m.w.e) (SNOEB = 999,999 mm H2O).Figures A1 and A2    Subbasins over 2000 m in altitude were divided into 10 elevation bands with 100-m to 200-m intervals, depending on the elevation range of the subbasin.Smaller elevation band intervals enable SWAT to model the glacier boundaries more accurately.It was assumed that the glacier boundary in a subbasin matched the lowest altitude of the elevation band if more than 50% of the elevation band area was covered by a glacier.

Model Calibration and Validation
The model was calibrated using monthly streamflow data at the gauge stations listed in Table 1, focusing on glaciated subbasins.The non-glaciated and snow-free gauged watersheds located at lowlands and glaciated subbasins with insufficient flow data, as well as dammed subbasins, were ignored in the model calibration.The model was validated using the monthly streamflow from the Rhone, Mendoza, and Chile for two periods: recent years from 2008 to 2010, and early years from 1982 to 1992.This was to confirm that the model works well for different periods.For the Vakhsh and Nepal River Basins, model validation periods were variable based on data availability.
Monthly calibration was performed using the SWAT Calibration and Uncertainty Program (SWAT-CUP) [39].SWAT-CUP is a generic interface to provide a link between any automatic calibration/uncertainty or sensitivity program and SWAT.Sequential Uncertainty Fitting Version 2 (SUFI2) incorporated in SWAT-CUP allows for calibration and validation of the model at multiple gauging stations and multiple observed datasets simultaneously.SUFI2 calculates a weighted multi-component objective function (weighted summation of the square errors) based on simulated variables and observed time series of all gauged watersheds and then minimizes it by searching the Latin Hyperbolic [40] generated parameters for the best solution.
In order to perform automatic calibration by SUFI2, the initial parameter values and ranges were determined.A list of all snowfall/melt parameters and associated groundwater parameters and   Subbasins over 2000 m in altitude were divided into 10 elevation bands with 100-m to 200-m intervals, depending on the elevation range of the subbasin.Smaller elevation band intervals enable SWAT to model the glacier boundaries more accurately.It was assumed that the glacier boundary in a subbasin matched the lowest altitude of the elevation band if more than 50% of the elevation band area was covered by a glacier.

Model Calibration and Validation
The model was calibrated using monthly streamflow data at the gauge stations listed in Table 1, focusing on glaciated subbasins.The non-glaciated and snow-free gauged watersheds located at lowlands and glaciated subbasins with insufficient flow data, as well as dammed subbasins, were ignored in the model calibration.The model was validated using the monthly streamflow from the Rhone, Mendoza, and Chile for two periods: recent years from 2008 to 2010, and early years from 1982 to 1992.This was to confirm that the model works well for different periods.For the Vakhsh and Nepal River Basins, model validation periods were variable based on data availability.
Monthly calibration was performed using the SWAT Calibration and Uncertainty Program (SWAT-CUP) [39].SWAT-CUP is a generic interface to provide a link between any automatic calibration/uncertainty or sensitivity program and SWAT.Sequential Uncertainty Fitting Version 2 (SUFI2) incorporated in SWAT-CUP allows for calibration and validation of the model at multiple gauging stations and multiple observed datasets simultaneously.SUFI2 calculates a weighted multi-component objective function (weighted summation of the square errors) based on simulated variables and observed time series of all gauged watersheds and then minimizes it by searching the Latin Hyperbolic [40] generated parameters for the best solution.
In order to perform automatic calibration by SUFI2, the initial parameter values and ranges were determined.A list of all snowfall/melt parameters and associated groundwater parameters and their ranges [23] are presented in Omani et al. [22].Since most of the watersheds studied are glacier dominated and above the tree lines, the land use and soils in the upper catchment had little or no significant impact on hydrology such as evapotranspiration (ET).So, none of the land use or soil parameters were found to be sensitive in the model calibration process, which confirmed that the roles of these parameters are not important to the water budget.The groundwater delay (GW-DELAY) and baseflow recession constant (ALPHA-BF) were initially set for simulation of low flow during the winter months.The parameter values then were optimized during the automatic model calibration.The baseflow recession constant is directly proportional to groundwater flow response to changes in recharge.
The parameters were adjusted for the elevation bands above the ELA 0 as an accumulation zone or glacier and below the ELA 0 as an ablation zone with seasonal snow.The melt parameters for the elevation bands were adjusted in order to match the observed and simulated average monthly flow curves and then the parameters were optimized by automatic model calibration using SUFI2.Briefly, model calibration consisted of two main steps.First, parameters were automatically calibrated for an entire basin so that the snow parameters were uniform throughout the entire basin.In the next step, the results were improved by calibrating the model for snow parameters at the elevation band-subbasin scale.Model performance was tested in some of the smaller subbasins first and then extended into the larger subbasins to test the hypothesis that the hydrologic significance of meltwater may be negligible at the macro scale despite the presence of large glaciers in the headwaters area [41,42].
Calibration and validation were evaluated using the root mean square error-observations standard deviation ratio (RSR), the Nash-Sutcliffe Coefficient of Efficiency (NSE) [43] and the Percent Bias (PBIAS).RSR is a ratio of the RMSE and standard deviation of measured data.It ranges from the optimal value of 0 to a large positive value [44].NSE indicates how well the plot of observed versus simulated data fits the 1:1 line.PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts [45].The optimal value of PBIAS is 0.0, with low values indicating accurate model simulation in term of magnitude.

Results and Discussion
Tables 3-7 show an average of adjusted values and range of the calibration parameters throughout the river basins for lumped and modified snow algorithms.The symbol "s" in the tables stands for glacier-free elevation bands or snowy elevation bands and "g" stands for glaciered elevation bands over the ELA 0 altitudes.The calibration parameter values are averaged for the elevation bands and subbasins and the parameter range shows minimum and maximum parameter values based on elevation band-subbasin scale throughout the entire river basin.Mean SMFMX (Maximum melt factor) throughout all river basins is between 2.80 and 5.38, and 3.30 to 7.71, for snow and glacier, respectively.TIMP (Snow temperature lag factor) values are lower in the high elevation bands where the glaciers exist and range between 0.010 to 1.000 with average values between 0.025 for high elevation bands and 0.740 for low elevation bands.It can be seen from Table 3 that the average values for SMFMX and SMFMN (Minimum melt factor) for the Narayani River Basin are larger than the obtained values for Rhone River Basin (Table 4) as also was indicated by Kayastha et al. [46].Notes: a l: elevation bands below 3300 m; b s: elevation bands between 3300 and 4300 m; c g: glaciered elevation bands over 4300 m.
Figure 4 shows how the simulated monthly flow and mean monthly flow fit the observed mean monthly flow cycle using the modified algorithm for some of the selected reaches.This degree of accuracy in simulation of the seasonal pattern of monthly flow is only achievable by adjusting melt parameters based on elevation bands.As demonstrated in Figure 5, SMFMX changes in lower elevation bands had more influence on the rising limb of the flow curve, whereas the descending limb was more sensitive to SMFMX changes in the higher elevation bands (glaciered areas).This indicates that the late spring flow is under the influence of snowmelt at lower elevations and that late summer flow is controlled by glacier melt at higher elevations.This can also be investigated by analyzing the melting lag time at the elevation bands.Seasonal melting from October 1998 to October 1999 from each of the elevation bands of a typical subbasin with a wide range of elevations is presented in Figure 6a.Melting for glacier-free elevation bands with seasonal snow cover are presented in Figure 6b for high elevation bands with permanent snow cover or glacier.In Figure 6a, the snowpack in elevation bands 1, 2, 3, 4, and 5 has completely vanished by mid-August, whereas in Figure 6b the glacier ablation reaches its peak in August and continues through October.The melt was decreased at lower altitudes by decreasing the SMTMP (snow melt base temperature) and SFTMP (snowfall temperature) at lower elevation bands.Calibration and validation results for all river basins are given in Tables 8 and 9.In the watersheds with highly glaciered drainage areas (Reaches 2 and 23) in the Rhone River Basin, PBIAS was improved when the modified snow algorithm was used instead of lumping all of the parameters for entire watershed area Both RSR and NSE also were improved using the modified model, although the improvement is more significant in simulation of magnitude of flow.There is negligible improvement in the simulated flow from watersheds with small areas of glacier (Reach 4) when using the modified snow algorithm in comparison to the lumped snow algorithm.For Narayani River Basin, PBIAS and NSE values in Table 8 reveal similar model performance by both lumped and modified snow algorithm methods, possibly as a result of the diminishing influence of glacier/snowmelt on runoff resulting from coincident summer monsoon precipitation in the Narayani River Basin.This monsoonal influence leads to accurate prediction of seasonal flow by rainfall-runoff models without explicitly considering the snow hydrology process.In Table 8, positive PBIAS in simulation of flow magnitude from all gauge stations (except Reach 96) means that the model has underpredicted the volume of flow.This systematic bias might be due to an underestimation of the NCEP reanalysis spring/summer precipitation (May-September) in the Narayani River basin of 20%.
The modified snow algorithm performed very well in simulation of monthly flow from the glaciered areas of the Vakhsh River Basin with RSR smaller than 0.5, PBIAS smaller than 10 percent and NSE greater than 0.70.Major difficulties in simulation of monthly flow from the Vakhsh River Basin were a lack of data for calibration and the low quality of the climate data, especially in the northern half of the river basin.The short calibration period of 2.5 years for Reach 72 was not long enough to capture the long-term variability of the monthly flow.The simulated mean monthly flow volume and cycle from the glaciated areas (drainage area to Reaches 72 and 133) by the modified snow algorithm during the summer months is in better agreement with the observed data in comparison to the mean monthly flow simulated using the lumped model for Reach 133 (Figure 4).At the main outlet (Reach 109), only the monthly flow variation shows improvement.Consequently, for an entire river basin it cannot be concluded that the modified model has better performance than the lumped model.Comparison between observed and simulated monthly flow from the main outlet of the basin during the validation period indicated that there was good agreement, verified by RSR, PBIAS, and NSE of 0.52, −8.15, and 0.73, respectively.
Among the gauged watersheds of the Mendoza River Basin, the drainage area of Reach 84 contains the largest area covered by glaciers.Subbasin 90 is a single gauged small subbasin, so it is a good example of the streamflow response to melt parameter distribution.There was a considerable improvement in the model performance when using the modified model instead of the lumped model in Reaches 79, 84, and 90. Figure 4 shows that the peaks simulated by the modified model match the observed peak flows.This level of accuracy is only achievable by adjusting the melt parameters for seasonal snow and permanent snow for each elevation band since the lumped model was not able to capture the peaks by calibrating the model for many different sets of melt parameters.Although the results from the modified model show significant improvement in highly glaciated areas and negligible improvement in less glaciated areas, the simulated monthly flow at the main outlet (Reach 86) does not show any improvement with the modified model.It can be seen in Table 8 that the model performance in simulation of downstream flow (e.g., Outlet 86 in the Mendoza River basin) declines when applying the modified snow algorithm.The new approach enabled calibration of the model versus the flow from upland catchments and optimization of the simulation results by adjusting the parameters separately for subbasin/elevation bands, while in the lump method the parameters were adjusted in order to get the reasonable result and were not necessarily the best at the gauged streams.It seems that the parameter value optimization for upland basins negatively affects the downstream flow in comparison to the lump parameter adjustment method.In other words, factors other than snowmelt could affect the downstream flow, and further model setup and parameter adjustment (groundwater, plant growth and soil moisture parameters, etc.) are necessary for lowlands and flat areas, which were out of the scope of this study.Another reason for declining the accuracy of downstream flow simulation could be the use of inaccurate reanalysis precipitation data.While snowfall/melt is the dominant hydrologic process in dry highlands (highland areas of Mendoza are very dry), flow is less affected by precipitation events, and snow parameter adjustment improves the flow simulation considerably.In lowlands, where the rainfall-runoff model is dominant, the flow is directly under the influence of precipitation events and any problem in input precipitation data could directly decline the accuracy of downstream flow simulation.So, when the results for snowy basins show improvement, flow in lowland areas may not necessarily improve.Comparisons between observed and simulated monthly streamflow during the validation period indicate range of model performances from good to unsatisfactory.
Comparisons between observed and simulated monthly flow from Chilean river basins during the calibration period indicate poor to satisfactory simulation.The RSR and NSE values are generally lower than those values obtained in simulation of monthly flows from the other river basins in this study.Winter precipitation (October to March) inter-annual variability in the Andes is linked to ENSO (El Niño-Southern Oscillation) events and consequently reveals a more complex response of streamflow [47], which may explain poor model performance in simulation of seasonal flow variability in the Mendoza River Basin and the Chilean river basins in the dry central Andes.

Conclusions
Treating the glaciered and un-glaciered areas in the watersheds separately improved SWAT model performance significantly in simulation of volume and variation of runoff in glaciered areas.Spatial and temporal variations of melt rates mainly depend on the spatial and temporal variations of melt factors in hydro-glacial models.While temporal variations of melt factors have been considered in the SWAT model in the past, there has been no consideration of spatial variations in melt factors and lag time factor, which are directly influenced by surface type (snow and ice).In this study, these spatial variations were specifically taken into account.
Model performance using the snow algorithms also depends on the climate of a river basin.Significance of melt water may be negligible in watersheds where the melt season coincides with monsoon precipitation, resulting in no significant difference in the simulation results from the two melt algorithms when monsoons were a factor.For the river basins in the central Andes, applying the modified distributed snow algorithm considerably improved the model performance in simulation of runoff volume in comparison with the lumped snow algorithm, while variation of runoff did not show any improvement.This may be due to high inter-annual and annual variability of flow in these regions, which are more dominated by rainfall-runoff relationships than snowmelt.We can conclude that considering the spatial variations of associated melt parameters significantly improves the SWAT performance in simulation of runoff volume and its seasonal variation in highly glaciated river basins.
While distributed process-based energy budget models have been tested in the SWAT model, no studies have been done to incorporate enhanced temperature-index models into SWAT.Incorporation of an enhanced temperature-index model into the modified snow algorithm of SWAT should be an objective of future studies on enhancing the SWAT snow hydrologic process.The modified snow algorithm did not enhance SWAT model performance in simulation of the streamflow at the main outlet of the river basins.Based on this, one might be tempted to draw the conclusion that at the macro scale the lumped model is preferable to the distributed model.This, however, requires further investigation into modeling such hydrologic processes as glacier surface mass balance and transient snow line altitude.Therefore, it is suggested that the model not only be calibrated using runoff but also be calibrated against the individual snow hydrology components.

Figure 2 .
Figure 2. Subbasin delineation, glacier outlines, and locations of dams and flow gauge stations in (a) Rhone River Basin; (b) Narayani River Basin; (c) Chilean River Basins; (d) Mendoza River Basin; and (e) Vakhsh River Basin.

Figure 2 .
Figure 2. Subbasin delineation, glacier outlines, and locations of dams and flow gauge stations in (a) Rhone River Basin; (b) Narayani River Basin; (c) Chilean River Basins; (d) Mendoza River Basin; and (e) Vakhsh River Basin.

Figure 3 .
Figure 3. Eight-day snow cover area variations extracted from MODIS products (MOD10A2) from 2000 to 2010 are represented as a percentage of the total area of river basins: (a) Upper Rhone; (b) Vakhsh; (c) Narayani; (d) Central Chile; and (e) Mendoza.

Figure 3 .
Figure 3. Eight-day snow cover area variations extracted from MODIS products (MOD10A2) from 2000 to 2010 are represented as a percentage of the total area of river basins: (a) Upper Rhone; (b) Vakhsh; (c) Narayani; (d) Central Chile; and (e) Mendoza.

Figure 4 .
Figure 4. Observed and simulated monthly runoff and mean monthly runoff using SWAT with different melt algorithms for the calibration period at selected reached.

Figure 4 .
Figure 4. Observed and simulated monthly runoff and mean monthly runoff using SWAT with different melt algorithms for the calibration period at selected reached.

Figure 4 .
Figure 4. Observed and simulated monthly runoff and mean monthly runoff using SWAT with different melt algorithms for the calibration period at selected reached.

Figure 4 .
Figure 4. Observed and simulated monthly runoff and mean monthly runoff using SWAT with different melt algorithms for the calibration period at selected reached.

Figure 6 .
Figure 6.Seasonal melting from the elevation bands from October 1998 to October 1999: (a) glacier free elevation bands with seasonal snow cover; (b) high elevation bands with permanent snow cover or glacier.

Figure A2 .
Figure A2.Glacier snow depth and area across the (a) Mendoza, (b) central Chile, and (c) Vakhsh river basins.Minimum snow cover area from MOD10A2 at the end of the melting season from February 2002 to 2010 was considered as glacier or permanent snow cover for Mendoza, central Chile, and part of the Vakhsh basin.Green dots show locations of glaciers with no depth information.

Figure A2 .
Figure A2.Glacier snow depth and area across the (a) Mendoza, (b) central Chile, and (c) Vakhsh river basins.Minimum snow cover area from MOD10A2 at the end of the melting season from February 2002 to 2010 was considered glacier or permanent snow cover for Mendoza, central Chile, and part of the Vakhsh basin.Green dots show locations of glaciers with no depth information.

Table 1 .
All available flow gauge stations with drainage areas for each river basin.

Table 2 .
The calculated area of glaciers using GLIMS, MODIS, and modeled glacier area in this study.

Table 2 .
The calculated area of glaciers using GLIMS, MODIS, and modeled glacier area in this study.
in Appendix A show the measurement locations of glacier depth and area from the WGI data inventory for the river basins.

Table 3 .
Range and mean values of calibration parameters in subbasin-elevation band scale, Narayani River Basin.

Table 4 .
Range and mean values of calibration parameters in subbasin-elevation band scale, Rhone River Basin.
Notes: a s: elevation bands below 2900 m; b g: glaciered elevation bands over 2900 m.

Table 5 .
Range and mean values of calibration parameters in subbasin-elevation band scale, Vakhsh River Basin Notes: a s: elevation bands below 4300 m; b g: glaciered elevation bands over 4300 m.

Table 6 .
Range and mean values of calibration parameters in subbasin-elevation band scale, Mendoza River Basin.
Notes: a s: elevation bands below 4700 m; b g: glaciered elevation bands over 4700 m.

Table 7 .
Range and mean values of calibration parameters in subbasin-elevation band scale, Chile River Basin.

Table 8 .
Model performance by both lumped and modified snow algorithms for calibration period.