Incorporating snow model and snowmelt runoff model for streamflow simulation in a snow-dominated mountainous basin in the western Hindukush-Himalaya region

: A Snow Model (SM) using a temperature-index method was used to optimize the degree-day factor ( DDF ) and pre‐ cipitation gradient ( PG ) for the different elevation zones of the Panjshir sub-basin for snowmelt runoff modelling. The values derived for DDF and PG were calibrated and vali‐ dated by comparing observed snow cover area and snow cover area simulated by SM. The Snowmelt Runoff Model (SRM) was used to simulate daily runoff over the hydro‐ logical years 2009–2014 using the optimized values for SRM accuracy. The optimized DDF values were 0.3 to 0.9 (cm °C –1 d –1 ) for elevations from 1593 m to 5694 m. Mean‐ while the PG was +0.002 m –1 for elevations 1593–4000 m and 0 m –1 above 4000 m. The simulated runoff by SRM during the entire data period correlated very well with a Nash-Sutcliffe coefficient NS = 0.93 utilizing both observed and simulated snow cover area. This method not only evaluates the characteristics of snowfall and snowmelt in different elevation zones to obtain the DDF and PG , but can also estimate the snowmelt runoff.


INTRODUCTION
Snow and glacier melt are key elements of the water cycle and play a major role as the lifeblood for people in downstream areas of mountainous environments. Snow not only affects the radiation balance of the Earth, but also plays as a major role in stream flow (Jain et al., 2008). The Himalayan range and the Tibetan plateau are the origin of Asia's main river system (Immerzeel et al., 2009). Upstream snow and ice storage in these river catchments are important for maintaining the availability of seasonal water (Immerzeel et al., 2010). The economics and the livelihoods of most people living within the Hindukush-Himalayan (HKH) countries are extremely waterdependent, and agriculture accounts for nearly 90% of all water withdrawals (Shrestha et al., 2015a). Increasing populations, agricultural activity, and urban development and industrialization lead to higher water requirements (Liniger et al., 1998). Some studies suggest that the water resources in the Indus River Basin (IRB) derived from snow and ice melt contribute about 50-80% of the annual flow and come from the western Himalayan ranges (Jeelani et al., 2012;Shrestha et al., 2015b). In areas where the water supply is presently dominated by the melting of snow or ice, rising ground temperatures may have significant hydrological implications (Jeelani et al., 2012).
The Snowmelt Runoff Model (SRM) has been developed to simulate and forecast daily stream flows in mountainous basins (Martinec, 1975). SRM has been successfully applied in the Himalayan region and the Tibetan plateau (Immerzeel et al., 2009;Zhang et al., 2014). The spatial distribution of snowfall and snowmelt, including their altitudinal characteristics, is essential for snowmelt runoff modelling in high mountain catchments. Previous studies (Tahir et al., 2011;Zhang et al., 2014;Hayat et al., 2019) have assumed parameter values from the literature or during model calibration for runoff simulations. But the topography, geography and geology of regions differ, which affects climatic conditions and runoff generation. Therefore, alternative methods, such as the Snow Model (SM) are being explored to estimate the spatial snowfall and snowmelt characteristics using degree-day factor (DDF) and precipitation gradient (PG) by making use of remotely sensed data (Azizi and Asaoka, 2019). The SM calculates changes in snow water equivalent based on a degree-day approach (Asaoka and Kominami, 2013). SRM and SM can be used to estimate the snow water equivalent and to assess spatiotemporal changes in snow cover area (Kazama et al., 2008;Whitaker and Yoshimura, 2012). SM has the ability to assess the impact of reconstructed snowfall on the catchment snow hydrology (Asaoka and Kominami, 2012). In contrast to other HKH areas, the Hindukush mountain range has received little scientific focus in terms of spatial and temporal snow cover variation and streamflow simulation in Afghanistan river basins due to climatic datascarcity. On the other hand, no studies on the HKH region have been conducted to date to combine SM and SRM to evaluate snow characteristics and melt runoff. This study focused on: (1) optimization of the DDF and PG, and estimation of snow cover area, (2) simulation of the river flow using simulated snow cover area, and (3) evaluation of the capability of the SRM using optimized DDF and PG to carry out the simulations of stream flow for the study area. The novelty of this research includes the combination of two models to decrease the uncertainty of the parameters for runoff modelling in a data-deficit mountainous environment. Moreover, this approach can evaluate the altitudinal characteristics of snowfall and snowmelt over high elevation and data-scarce mountainous areas by obtaining DDF and PG. First this study integrates remotely sensed Snow Cover Area (SCA) and the SM in order to optimize the DDF and PG for different elevation zones to simulate runoff with good accuracy. Second it estimates daily melt runoff utilizing the integration of SRM with both observed and simulated SCA.

Study area
The Kabul river basin, with a semi-arid and strongly continental climate, has several tributaries including small and large rivers originating from high altitude areas of the basin ( Figure 1a). Its main tributaries include the Logar, Panjsher (with its own major tributary the Ghorband), Laghman and Kunar rivers. The basin also includes the Kabul urban area, which hosts around 20% of the country's population solely in the capital city (CSO-IRoA et al., 2019) and consumes a large fraction of the installed energy generation capacity as well as withdrawing considerably more water than the rest of the country. The Panjshir sub-basin located within the Kabul basin was selected for this study (Figure 1a). The study area has a total area of 3540 km 2 with an elevation range from 1593 m to 5694 m (a.s.l.) and a mean catchment elevation of 3498 m (a.s.l.) estimated based on a hypsometric curve (Figure 1b). The Panjshir basin physiography includes high Hindukush mountain ranges dominated by seasonal snowpack and glaciers.
The data set used consists of remotely sensed satellite data (ASTER GDEM and MODIS cryosphere product) and field observations (stream flow and climate data), with a summary description included in text S1.

Snow Model
The Snow Model has been used to simulate the distribution of snowfall, SCA, snowmelt, and snow water equivalent spatially and temporally. The spatial and temporal resolution of the SM output is based on the input DEM and climate variables. This model simulates snowmelt based on a degree-day approach. We can estimate Snow Water Equivalent (SWE) by combining snowfall and snowmelt, and it is possible to calculate the temporal changes in SWE utilizing the following equation: where SWE is the snow water equivalent (mm), SF is the daily snowfall (solid precipitation) at each pixel (mm d -1 ), and SM is the daily snowmelt rate (mm d -1 ) and t represent the date of snow water equivalent, and dt is the time interval of one day. Daily snowmelt rate was calculated based on the degree-day method (Martinec, 1975) as: where SM is the daily snowmelt rate (mm d -1 ), DDF is the degree-day factor (mm °C -1 d -1 ), T a is the average daily temperature (°C). A base temperature is defined as 0°C, below which snowmelt does not occur. Daily precipitation at the measuring points were interpolated over each grid by: where P is daily precipitation (snowfall/rainfall) (mm d -1 ) at the elevation h (m), PG is the coefficient of the vertical gradient (m -1 ) of precipitation with respect to elevation h, and elv obs is observation elevation. The "obs" subscript indicates the ground observation station. To separate snowfall and rainfall events we used the critical temperature. The algorithm and parameters of SM have been discussed in detail by Azizi and Asaoka (2019), and a summary can be found in text S2.

Snowmelt Runoff Model
Snowmelt Runoff Model (SRM) is a conceptual, semidistributed and degree-day-based hydrologic model developed by Martinec (1975). SRM was developed to simulate and forecast daily flows resulting from snowmelt, glacier melt, and rainfall in mountainous environments where snowmelt is the principal contributor to runoff (Martinec, 1975;Martinec et al., 2008). The SRM has proven to be efficient for use in Himalayan river basins where hydrometeorological networks are limited (Immerzeel et al., 2009). The following equation is used in SRM to calculate the daily discharge (Zhang et al., 2014).
where Q is the mean daily discharge (m 3 s -1 ), Qs n and Qr n are the proportions of snowmelt and rainfall runoff in the nth day respectively; C s and C r are the runoff coefficients for snow and rain respectively; a is the degree-day factor (cm °C -1 d -1 ); T is the degree-days (°C d); ΔT is the adjustment by temperature lapse rate when extrapolating the temperature from observation points to the average hypsometric elevation of the zone (°C d); S is the ratio of the snow covered area to the total area; P is the measured precipitation on that day (cm); A is the area of basin or zone (km 2 ); k is the recession coefficient (derived from parameters X c and Y c ); and n is the sequence of days during the computation period. The factor of 0.1157 converts data from cm km 2 d -1 to m 3 s -1 . Martinec et al. (2008) provided a full description of the factors and parameters used in the formula (4). Mean daily temperature, daily precipitation, and daily observed and simulated SCA percentage were used as input variables for model operation. The daily observed SCA was obtained from a linear interpolation based on adjacent days. However, daily simulated SCA was found through SM as discussed in detail in text S2. ASTER GDEM was used to delineate the basin characteristics such as zonal areas and the hypsometric/area elevation curve of the study area. To obtain the spatial climate heterogeneity common for mountainous basins, the Panjshir sub-basin was divided into eight different altitudinal zones based on the hypsometric curve ( Figure 1b). Daily mean temperature data was extrapolated for each elevation zone using a calculated lapse rate value of -0.007°C as discussed in text S2. On the other hand, daily precipitation in each zone was estimated using the precipitation gradient values presented in Table SI. SCA derived from SM (simulated) and MODIS (observed) are then extracted for the eight zones, respectively. The observed and simulated SCA details are given in text S1 and text S2, respectively. The model is calibrated (October 2009-September 2011) and validated (October 2011-September 2014) using daily discharge data measured in the outlet of the Panjshir sub-basin. Calibration parameters required for SRM operation are given in Table SII. SRM parameter values for the degree-day factor were estimated through SM for each elevation zone and rest were either defined from the available observed data (Text S2) or at the time of model calibration ( Table I). The model efficiency was evaluated using the Nash-Sutcliffe determination coefficient (NS) and the volume difference (DV). Further evaluation of the SRM was carried out with the Pearson correlation coefficient and the Root Mean Square Error (RMSE) to determine the correlation between the measured and simulated daily discharges (Table II).

A.H. AZIZI AND Y. ASAOKA
Landsat and MODIS satellite remote sensing products are widely used as input for snowmelt runoff modelling for current and future runoff simulations (Martinec and Rango, 1986;Dey et al., 1989;Tahir et al., 2011;Adnan et al., 2017;Azmat et al., 2017;Hayat et al., 2019), although these products have limitation in snow cover estimation due to their spatio-temporal resolution and cloud cover. There are also limitations with assuming the future SCA for runoff simulation. This will cause errors for future runoff projections. To avoid ambiguities, our approach can reduce uncertainties in runoff modelling with the simulation of daily SCA and future SCA using climate model projections through SM.   (Tahir et al., 2011;Hayat et al., 2019). In addition, we found that the PG is positive at altitudes 1593-4000 m (+ 0.002 m -1 ) while becoming zero in high elevation zones above 4000 m. A previous study (Immerzeel et al., 2012) also found a strong positive vertical precipitation gradient value of 0.0021 m -1 in the upper Indus basin. Since the climatic stations in the Panjshir subbasin are valley-based, the precipitation in high mountain regions can be underestimated (Buda et al., 2016). Therefore, many studies have applied equation (3) under the assumption that snowfall increases with elevation. However, we found that the snowfall is constant with elevation in higher altitudinal zones. We suspect this might be due to the topography and steep slopes in high elevations (above 4000 m), which may cause the redistribution of snowfall due to the wind velocity and snowdrift.

Snow cover simulation using SM
Daily snow cover distribution and SCA estimated by SM were evaluated through the comparison between the simulated SCA and observed snow cover values over the calibration and validation periods 2009-2010 and 2010-2014 respectively, as presented in Figure 2. Commonly, simu- lated SCA is in good agreement with the SCA obtained from the MODIS 8-day data set with RMSE = 162 km 2 . In September 2010 and 2012, SM simulated some SCA with step shape changed over time due to some precipitation events assumed as snowfall. In order to validate the calibrated DDF and PG values for each elevation zone, we further examined the simulated and observed SCA to evaluate the correlation between these two snow cover values in the eight different altitudinal zones. The SM overestimated snow coverage in high elevation zones above 4000 m, during the end of the melt season (July-August) in the years 2010, 2012, 2013 and 2014. However, in the year 2011, the snowpack was underestimated at both high and low elevations by the SM. Some overestimation may be due to the simulations of rainfall as snowfall in high elevation zones in the spring and summer seasons ( Figure S2). On the other hand, these over-and under-estimations may be attributed to the valley-based climate stations and observed SCA limited by an 8-day interval or influenced by cloud cover. Snow covers nearly the entire basin in late winter and early spring. In summer, at the end of the melt season when temperature rises, the SCA depletes to a minimum ( Figure S3). Overall, we could successfully optimize zone-wise DDF and PG, which leads to better simulation of daily SCA.

Runoff simulation using remotely sensed snowcovered area
To simulate the daily discharge, the SRM was calibrated and validated over the hydrological years October 2009 to September 2011 and October 2011 to September 2014, respectively. The zone-wise estimated values and parameters used for runoff simulation are given in Table I. Figure  3 shows the observed and simulated flow hydrographs in the calibration and validation period. The model performance was examined using statistical analysis (Difference of volume (D v ), Nash-Sutcliffe coefficient (NS), Root mean square error (RMSE) and Pearson Correlation coefficient) derived from the comparison of the observed daily discharge values with simulated values over both the calibration and validation periods using observed SCA presented in Table II. Overall, the assessment criteria showed that the SRM is able to successfully simulate the daily discharges in the Panjshir sub-basin. The observed and simulated flow hydrographs showed that the model simulated high and low flows very well during the calibration and validation period (Figure 3). Figure S4a shows a scatter plot of the daily observed and simulated discharge, which resulted in better agreement between the low and mid-level flows, whereas

Runoff simulation using simulated snow-covered area
To simulate the daily discharge, the SRM was calibrated and validated over the same hydrological periods as in the observed SCA. The common values and parameters for each elevation zones provided in Table I were used for discharge simulation using the simulated SCA. Figure 4 shows the observations and simulations of the daily discharge for all years from October 2009 to September 2014. Table II shows the relationship between the measured and simulated daily flow using different efficiency criteria. The model performance showed good correlation between measured and simulated discharge for all years using simulated SCA. The mean NS coefficient and DV were found to be 0.91 and -3.34%, respectively. The Pearson Correlation coefficient was found to be ≥0.95, and the lowest and highest values for RMSE were 10.7 and 28.5, respectively, over the annual simulation periods. The NS and the Pearson Correlation coefficient values were ≥0.64 and ≥0.88 over the snowmelt season, respectively, but the lowest NS and Pearson Correlation coefficient values were found to be -0.09 and 0.77 during the extreme events period. Figure  S4b shows daily measured and simulated flow plotted against each other to evaluate their correlation, which indicated a very good relationship over the entire streamflow duration. However, in contrast with the use of observed SCA, high flow periods were overestimated in 2010 and  These results indicate that SRM simulated both the base flow and discharge peaks well with observed SCA, but some peak flows are underestimated. Conversely, with simulated SCA, the SRM simulated the high flows very well compared to the base flows.
In the case with simulated SCA, the discharge overestimation is attributed to the overestimation of SCA in high elevations zones, whereas the underestimation of discharge during the winter season may be due to the low ability of SRM to simulate slow melting as a consequence of ground heat exchange (Whitaker and Yoshimura, 2012).

CONCLUSIONS
In this study a combination of SM and SRM was used to reduce uncertainty in the estimation of degree-day factor and spatial precipitation distribution, including SCA, in the simulation of daily streamflow in a data-scarce region of Panjshir sub-basin in Kabul basin in the western Hindukush-Himalaya region. This combination approach can evaluate the spatial characteristics of snowfall and snowmelt by obtaining the DDF and PG. Some significant conclusions can be found from the results obtained in this study.
Daily temperature lapse rate was estimated using daily temperature data at five stations located at different elevations. Since it showed high daily variation, we made use of the mean lapse rate for this study. The critical temperature was calibrated and a value of 0°C was used. The DDF and PG are adjusted for the distribution of snowmelt and snowfall predictions in the eight altitude zones using the assimilation technique between the simulated SCA and the 8-day SCA obtained from the observed data.
We obtained DDF and PG values that varied with changes in elevation. The DDF values ranged from 0.3 to 0.9 (cm °C -1 d -1 ) from the low to high elevation zones.
The simulated SCA by SM was in good agreement with the SCA measured from the 8-day interval observed data. The SRM can efficiently simulate daily discharge in datascarce mountainous environments such as Panjshir subbasin. With consideration of the calculated lapse rate and calibrated critical temperature including optimized DDF and PG daily stream flow simulated, the SRM performed well with the integration of both observed snow cover product and the simulated SCA by SM. The efficiency of the model was never less than 0.89 and 0.87 for the NS coefficient for both the observed and simulated SCA, respectively. The high flow periods were more efficiently modeled using simulated SCA than observed SCA during both the calibration and validation periods. To reduce the uncertainties associated with the assumption and calibration of the parameters and increase the accuracy of the simulations, it is essential to accurately estimate the temperature and precipitation parameters and to reduce observation errors.
Studies have shown that global warming can strongly influence snowfall, snowmelt, and SWE distribution, which can adversely affect the snow and glacier storage. To precisely evaluate the impact of this phenomenon, future work is needed to estimate snowfall distribution and runoff simulation using climate model projections for better water resources management. The method presented in this research may effectively estimate future water resources in high elevation data-scarce mountain regions.

SUPPLEMENTS
Text S1. Data sets Text S2. Snow Model (SM) Text S3. Uncertainties and Limitations Figure S1. Daily lapse rate variation calculated using daily observed temperature Figure S2. Evaluation of the observed and simulated SCA in different elevation zones during the calibration period, the x-and y-axis represents the SCA (km 2 ) and elevation (m), respectively Figure S3. Evaluation of simulated spatial distribution of SCA with observed snow cover data during calibration period Figure S4. Scatter plot of daily observed and simulated discharge values using; (a) Observed SCA, and (b) Simulated SCA Table SI. Optimized values of DDF and PG using SM Table SII.SRM and SM parameters used for model operation