Impact of Catchment Discretization and Imputed Radiation on Model Response: A Case Study from Central Himalayan Catchment

: Distributed and semi-distributed hydrological modeling approaches commonly involve the discretization of a catchment into several modeling elements. Although some modeling studies were conducted using triangulated irregular networks (TINs) previously, little attention has been given to assess the impact of TINs as compared to the standard catchment discretization techniques. Here, we examine how different catchment discretization approaches and radiation forcings inﬂuence hydrological simulation results. Three catchment discretization methods, i.e., elevation zones (Hypsograph) (HYP), regular square grid (SqGrid), and TIN, were evaluated in a highly steep and glacierized Marsyangdi-2 river catchment, central Himalaya, Nepal. To evaluate the impact of radiation on model response, shortwave radiation was converted using two approaches: one with the measured solar radiation assuming a horizontal surface and another with a translation to slopes. The results indicate that the catchment discretization has a great impact on simulation results. Evaluation of the simulated streamﬂow value using Nash–Sutcliffe efﬁciency (NSE) and log-transformed Nash–Sutcliffe efﬁciency (LnNSE) shows that highest model performance was obtained when using TIN followed by HYP (during the high ﬂow condition) and SqGrid (during the low ﬂow condition). Similar order of precedence in relative model performance was obtained both during the calibration and validation periods. Snow simulated from the TIN-based discretized models was validated with Moderate Resolution Imaging Spectroradiometer (MODIS) snow products. Critical Success Indexes (CSI) between TIN-based discretized model snow simulation and MODIS snow were found satisfactory. Bias in catchment average snow cover area from the models with and without using imputed radiation is less than two percent, but implementation of imputed radiation into the Statkraft Hydrological Forecasting Toolbox (Shyft) gives better CSI with MODIS snow.


Introduction
Accurate runoff prediction is one of the fundamental challenges for effective and sustainable water resource management in the Himalaya region. Reliable meteorological forcing data and accurate land cover information at a range of spatial and temporal scales are critical to effective runoff prediction, and for water resource management [1]. The main challenge to estimate (or measure) meteorological forcing variables comes from the fact that they vary both in space and time as a function of meteorological inputs. Spatial variability of topography and land use cover types also adds to the complexity of the problem [2].
In the past, several hydrological models have been developed for use in various applications. Some of these models include SRM [3], HEC-HMS [4], J2000 [5], Statkraft's Hydrological Forecasting Toolbox (Shyft) [6], etc. Unlike other hydrological models, Shyft provides an optimized platform for the implementation of many well-known hydrological models from conceptual to physically based distributed hydrological models. Most of the distributed hydrological models discretize the catchment and derive terrain attributes from digital elevation models (DEMs). Most commonly used catchment discretization methods in hydrological applications are grid-based [7][8][9][10], and elevation band (hypsography) based [11][12][13]. Recently, triangulated irregular network (TIN)-based discretized hydrological models are also used for discharge simulations [14,15], but to date TIN-based distributed hydrological models are not tested in the Himalaya catchments. Although the use of grid and elevation based discretized models are more prevalent due to their simplicity and ease of implementation [16,17], TIN-based discretized models have a better potential of accurately representing topographic features than grids and hypsography [14]. In this context, Shyft can be a valuable tool to evaluate different models as it have the ability to take the terrain heterogeneity into account and provide distributed output variables. Shyft can be used to test different functioning hypotheses from which simplified and/or predictive models can be derived for more operational purposes. If reliable output variables are expected, Shyft parameters can be estimated a priori from available information. Verification of these models behavior can be performed not only on the discharge at the outlet but also at intermediate variables (such as snow cover area, snow water equivalent, etc.) leading to a multiobjective verification. The main drawback of the current version of the Shyft is that it uses conceptual models.
The topographical information such as snow cover area is one of the most important variables for the application of hydrological model [18][19][20], and can dominate local and regional climate, and hydrology. Snow and glacier strongly influence the regional radiation budget with its high albedo and acts as an insulation, controlling snow accumulation and melt [21]. In this regard, accurate representation of snow cover at micro-scale of spatial variation is required to provide accurate and timely information on water availability and its management. Many distributed hydrological models are commonly employed to quantify snow cover area and snow water equivalent (SWE). All models have strengths and limitations [22,23]. For example, current model-based snow simulation approaches in Shyft are limited to the radiation received on a horizontal plane surface, which introduce a certain degree of uncertainty in the radiation received by the surface. In view of this, we have implemented the imputed radiation algorithm following a similar approach as in Allen et al. [24] in order to account for the effect of slope and aspect on incident radiation.
Remote sensing is an alternative powerful tool that offers the ability to examine quantitatively the physical properties of snow [25]. A well-known example is a MODerate Resolution Imaging Spectroradiometer (MODIS) instrument which can provide daily snow cover fraction with nearly global coverage at the resolution of 500 km [26,27]. For regional snow cover mapping, the MODIS satellite sensors are particularly appealing due to their high spatiotemporal resolution [28]. In this regard, validation of snow cover area simulated from the Shyft with MODIS will also add the most important aspect to hydrological simulation.
Despite a number of previous modeling studies, the conversion of a highly variable meteorological field into a spatially-distributed snow response remains an unanswered, open question. This is primarily due to the fact that meteorological processes operating in a watershed are highly heterogeneous and interconnected in the system [29]. Although some modeling studies were conducted using TIN previously, little attention has been given to assess the impact of TIN as compared to the standard catchment discretization technique over the Himalayan catchment. In this regard, we have implemented imputed radiation following a similar approach to Allen et al. [24] on different catchment discretized models (i.e., elevation band, regular spaced grid, and TIN) for calculating snow cover area using PT_GS_K stack of Shyft (https://shyft.readthedocs.io). The primary method for the simulation and distribution of snow cover area in PT_GS_K is in the matrix format of equally spaced elevation points. A data point in a matrix does not allow meaningful analysis of the terrain. Hence we hypothesize that the Shyft using the TIN for discretizing a watershed should enable a better prediction of the change in a runoff pattern caused by a change in a watershed characteristics. Finally, we tried to answer the question, what value does the TIN format offer over a matrix format for this particular application. This study also presents and evaluates the different catchment discretization models and the newly added translated radiation algorithm on the snow cover area and discharge simulations. As of our best knowledge, TIN-based discretized model with imputed radiation has not been applied before to simulate discharge in the Himalayan river catchments.
This manuscript is novel providing a capability to simulate discharge over the Himalayan catchment using community modeling framework. This paper will help in an understanding of downscale and translated radiation measurements onto the inclined surface and its impacts on hydrological response. An understanding of the relationships between catchment discretization and hydrological response and a cross-validation of hydrological model snow simulation using MODIS derived snow cover data will help to fill existing knowledge gap over the region. Hydrological model calibration and validation was carried out between 2000 to 2009. Selection of calibration and validation years was based on the observations data quality (less missing data), particularly observed discharge. In the paper, the materials and method, study area, hydrological models, model forcing datasets, and validation datasets are presented in Section 2. Results are presented in Section 3, while the discussion is presented in Section 4. Finally, conclusions are presented in Section 5.

Hydrologic Model Framework
The Statkraft Hydrological Forecasting Toolbox (Shyft, https://gitlab.com/shyft-os) was used in this study. Shyft provides an optimized platform for the implementation of well known hydrological models [6]. The high-performance generic time series framework in Shyft allows for rapid calculations of hydrologic response at the regional to a cell scale. The standard models in Shyft, i.e., PT_GS_K, PT_SS_K, and PT_HS_K, mainly differ in their snow mass balance methods calculations. The main model forcing variables are atmospheric temperature, precipitation, wind speed, relative humidity, and shortwave radiation, given in Table 1. Shyft works down from the regional level to the cell level by distributing the input variables and parameters to the georeferenced grid cell (described in catchment discretization technique). In the model, Bayesian kriging approach is used to distribute point temperature data over the catchment, while the inverse distance weighting approach is used for other forcing variables. The PT_GS_K model is the widely used as compared to the remaining models of Shyft (Teweldebrhan, 2019). Potential evapotranspiration in this model is calculated by Priestley-Taylor method [30]. Actual evapotranspiration is assumed to take place from snow-free area and is estimated as a function of potential evapotranspiration, and a scaling factor. Gamma snow routine (Section 2.1.1) is used for snowmelt, sub-grid snow distribution, and mass balance calculation. A simple storage-discharge function is used for catchment response calculation [31,32]. The catchment response function is based on the storage-discharge relationship concept described in Kirchner [32] and represents the sensitivity of discharge to storage changes given by Equation (1). As discharge is generally nonlinear and varies by many orders of magnitude, the recommended approach is to use log-transformed discharge values to avoid numerical instability.
where P, Q, and E are the rate of precipitation (from snow-melt and rain), discharge, and actual evapotranspiration (from snow-free area) in units of depth per time. The concept behind this method is that the catchment sensitivity of discharge to changes in storage i.e., g(Q), can be estimated from the time series of the discharge alone through fitting empirical functions to the data such as the quadratic equation, and is given by Equation (2) [32].
g(Q) = e c1+c2(ln(Q))+c3(ln(Q)) 2 (2) where c1, c2, and c3 are the catchment specific outlet parameters (hereafter called Kirchner parameters) obtained during the model calibration. There is no routing into the model, but Kirchner parameters response delays outflow from precipitation and snow/glacier melt. The Gamma snow routine uses an energy balance approach (Equation (3)) for snow ablation, and the snow depletion curve represented by gamma distribution for snow distribution. The energy balance approach in the Gamma snow routine follows similar approach as used in DeWalle and Rango [33]. The net energy flux (∆E) at the surface available for snow ablation is expressed as where S is the net shortwave radiation; L in and L out are the incoming and outgoing longwave radiation; respectively, H SE and H L represent sensible and latent heat fluxes, respectively; and E G is the net ground heat energy flux and is calculated using bulk-transfer approach. Two parameters defining wind profile-intercept (wind constant) and slope (wind scale)-are determined through either calibration or prescribed (Table 2). For a given time step (t), the snow albedo (α) at each grid cell depends on the minimum (α min ) and maximum albedo (α max ) values as well as the albedo decay rate, temperature, and snowmelt e.g., Hegdahl et al. [9]: In Equation (4), FDR and SDR denote fast and slow snow cover decay rates, respectively. In this study, α min and α max are estimated during the model calibration.
Within the gamma snow routine, precipitation falling in each grid cell is classified as solid or liquid depending on a threshold temperature (tx), and the actual grid cell temperature. Snow distribution to each cell is estimated by using a three-parameter gamma probability distribution [7]. Snowmelt depth is calculated by multiplying ∆E (available energy) with the latent heat of fusion for water.
The temperature index model [34] is used in the model for glacier melt calculation. The glacier reservoir is assumed to be infinite and glacier area is assumed to be constant throughout the study periods. In a grid cell, snow is assumed to melt first from glacier cover/free areas, and glacier melt only happens from glacial areas that are snow-free. Obtained runoff depth per grid cell is summed with snowmelt depth from the gamma snow routine and converted to volume by multiplying the area of each grid cell before Kirchner routine. Table 2. Model calibration parameters with upper and lower limits. Prescribed parameters are marked with *. Parameters only used in radiation model are marked with **.

Parameter
Description and Unit Lower Limit Upper Limit kirchner parameter 3 (see Equation (2)

Study Catchment and Data
The study has been conducted in the Marsyangdi-2 River catchment ( Figure 1), Nepal. It is located between 28 • 10 21 N to 28 • 54 11 N and 83 • 47 24 E to 84 • 48 04 E. The catchment has a total area of 3001 sq. km, of which 80% lies above the elevation of 3000 m asl. The catchment elevation ranges from 641 m asl to 7924 m asl with a mean elevation of 4377 m asl, and mean slope of 29.34 • . Marsyangdi-2 river is a tributary of the Narayani river catchment ( Figure 1). It originates from the southern edge of the Tibetan Plateau, flows through Nepal to India, and drains into the Ganges River. Approximately 24% of area is covered by glaciers [35]. The hydrological regime is heavily influenced by seasonal weather variations. The seasons are defined as monsoon (June to September), post-monsoon (October to November), winter (December to February), and pre-monsoon (March to May) [36]. Approximately 80% of the total annual precipitation occurs during the monsoon period [37,38]. The Narayani River catchment lies also in the primary hydropower development region of Nepal, providing 44% of the total national electric supply [39]. Some of the operating hydropower development projects in this catchment include Marsyangdi (69 MW) and middle marsyangdi (70 MW). The upper marsyangdi (600 MW), lower Manang marsyangdi (100 MW), and Nyadi (30 MW) are some of the hydropower project currently under different stages of development. Therefore, the study area is an important catchment in Nepal from hydropower perspective.

Meteorological Forcing and Observed River Gauge Data
Standard model forcing data are temperature, precipitation, relative humidity, wind speed, and incoming short wave radiation. All these model forcing were obtained from Water and Global Change (WATCH) Forcing Data ERA-Interim (WFDEI). Selection of the forcing data sets for the hydrological model was based on the prior study by Bhattarai et al. [10]. WFDEI is a global meteorological forcing data set at 0.5 • × 0.5 • horizontal resolution obtained by bias-correcting ERA-I data [40,41]. It covers the period 1979-2016 and contains eight meteorological variables at a 3 h time step [42]. WFDEI then corrects ERA-Interim for precipitation biases using data from the Climatic Research Unit (CRU) or the Global Precipitation Climatology Center (GPCC). In this study, GPCC products were preferred to CRU because of their higher resolution and data quality [40,41]. The WFDEI data set is freely available online from ftp.iiasa.ac.at.
Daily observed river discharge data (observation periods: 2000-2010) at Marsyangdi_2 River 28.2036 • N, 84.403 • E) are obtained from the Department of Hydrology and Meteorology (DHM), Government of Nepal (GoN). The catchment is highly influenced by seasons as more than 70% of total annual discharge flows from June to September. Highest daily mean discharge is observed during August (386.7 m 3 /s) and low mean discharge is observed in February (30.5 m 3 /s), and daily mean observed discharge is 133.85 m 3 /s. Thus, the obtained discharge data from the DHM GoN was used for the model validation purpose.

MODIS Snow Product, Topographical, and Land Cover Data Sets
In this study, the standard eight-day MODIS Terra and Aqua snow cover products, MOD10A1 and MYD10A1, are chosen for the analysis. They contain snow cover, snow albedo, fractional snow cover, and quality assessment (QA) data in a compressed HDF-EOS format along with corresponding metadata. Only snow cover information is used in the study. Both Aqua and Terra products are with 463.3 m spatial resolution in sinusoidal projection. Frequent cloud cover is one of the major challenges when using MODIS and other optical remote sensing data in the Himalayan region [43]. MOD10A1 was successfully used in earlier studies [44][45][46]. However, Rittger et al. [47] and Toure et al. [48] showed that the MODIS data can be systematically biased compared to Landsat ETM+ snow cover data in the mountainous regions. The average root mean square error (RMSE) between Landsat ETM+ and MOD10A1 snow cover fraction over the Nepal Himalaya was 0.23 [48]. Thus, a composite data set was formed, following Muhammad et al.'s [49] approach, using data retrieved from the Aqua and Terra products in order to minimize the effect of clouds and other sources.
A land cover map (0.3 • × 0.3 • , resolution) providing forest cover, lake cover, glacier cover, and reservoir cover of the study area is extracted from the GlobCover Project [50]. Digital Elevation Model (DEM) of 90 m spatial resolution from the Shuttle Radar Topography Mission (SRTM) [51] is used in the study.

Catchment Discretization Technique
The study catchment was delineated using an automatic catchment delineation method described by Akram et al. [52]. The delineated catchment was further discretized into modeling elements using three methods.

Triangulated Irregular Network (TIN)
TIN is a vector-based data structure comprising a triangular network of vertices with associated coordinates (x, y, z) connected by edges. TIN utilizes the original sample points to constitute many nonoverlapping triangles that cover the entire region. Spacing and shape of the triangles are determined by the terrain and by the desired degree of fit. The mathematical basis for a TIN construction is the Delaunay triangulation methods described in Field [53] and Hjelle and Daehlen [54], and is based on a principle of maximizing the minimum angle of all triangles produced by connector lines to nearest neighbor points. The TIN interpolation method is an altitude interpolation algorithm, this makes it possible to distribute hydrometeorological forcing variables more correctly for hydrological simulations [55].
Rasputin (https://github.com/expertanalytics/rasputin) is used for generating TIN terrain models of the study area. It can convert a point set of coordinates to a TIN using Delaunay triangulation methods. Specifically, it has been developed to convert Digital Elevation Models (DEMs) into simplified triangulated surface meshes. Rasputin generates catchment TIN using the provided DEM and also calculates physical parameters of each TIN facet, such as slope, area, and aspect, from the TIN geometry. The Rasputin software requires wkt-file describing the catchment polygon, which is easily obtained from shapefiles generated by Qgis. Rasputin is freely available under GNU GPLv3 license.
In the study, TINs with different numbers of cells were generated with the help of Rasputin and DEM [56]. The generated TINs with different number of cells are used for the model calibration and validation. The Nash-Sutcliffe Efficiency (NSE; see Section 2.4.1) value was used for evaluating the performance of the models with particular numbers of TIN cells. Implementing the spatially distributed parameters for the TIN with different number of cells yielded different simulation performance when compared to the observations over the calibration (2000)(2001)(2002)(2003)(2004) and validation (2005-2009) periods. As expected, model efficiency in terms of NSE increased with increasing number of TINs ( Figure 2). However, after reaching a certain number of cells no further improvement in NSE was found, indicating the limit for Shyft to utilize terrain data. The number of cells corresponding to this limit was considered as an optimum i.e., 671 cells (with median area of 4 km 2 ) (see Figure 2).

Regular Space Grid
Vector creation tools under Qgis geoalgorithms were used for creating regular space grid cells. Due to the ease of comparison with TIN, a grid mesh with an area of 4 km 2 was selected for this study. Unlike TIN, each grid cell is connected to either four or eight adjacent neighbors, the single resolution inherent in raster grids. The physical parameters (slope, aspect, elevation, latitude, and longitude) of each cell at the centroid (center point) are calculated by using vector geometry tools under Qgis geoalgorithms and further used for distributing the hydrometeorological forcing.

Hypsography (Elevation Zone)
In this study, the catchment area was divided into ten elevation zones. Mean elevation, aspect, and slope of each zone were calculated and further used for the distributing the hydrometeorological forcing. Figure 3 shows discretization of the study catchment using TIN, regular space grid, and elevation bands.

Implementation of New Radiation Algorithm in the Shyft
Daily incoming shortwave radiation (from forcing data) falling over the inclined surface with specified slope and aspect was calculated using the methods described by Allen et al. [24]. The procedure assumes each grid/TIN surface has uniform slope. Thus, effects of protruding surrounding terrain were not considered. This simplification facilitates implementation of the procedure within Shyft in energy balance calculation for snow mass balance and in Priestley-Taylor-based potential evapotranspiration calculations.

Model Construction, Calibration, and Validation
Combination of different catchment discretization methods and radiation algorithms resulted in different models (Table 3). Model forcing datasets was interpolated to each geometric center of cells/TIN using Kriging (for temperature) and IDW methods (for precipitation, wind speed, relative humidity, and incoming shortwave radiation). Details about the interpolation methods are explained in detail in Bhattarai et al. [10]. Each model was calibrated and validated separately. Calibration and prescribed parameters are given in Table 2. Manual calibration can be time-consuming and subjective; therefore, an automatic calibration was carried out using CREST v.2.1, the Shuffle Complex Evolution University of Arizona (SCE-UA) [57]. As described in Bhattarai et al. [10], the calibration procedure involves the selection of samples in the parameter space through the use of competitive evolution schemes, such as the simplex scheme, to reproduce better the observations. After several iterations, either due to convergence or when the maximum number of iteration is reached, the best set of parameter values based on the Nash-Sutcliffe coefficient of efficiency (NSE) is determined. The model was validated with the observed discharge. A cross-validation technique was used, where the optimal model identified in the first periods, i.e., years 2000-2004 was validated during the second period, i.e., 2005-2009, and vice-versa. For convention, different names were given to the model-discretization combinations (Table 3).

Nash-Sutcliffe Efficiency (NSE)
The NSE [58] was used (Equation (5)) to quantitatively evaluate the degree of agreement between observed and simulated discharge. NSE is commonly used to assess the predictive power of hydrological models and its value can range from −∞ to 1. An efficiency of 1 corresponds to a perfect match of simulated to observed discharge values. However, the NSE calculated using raw values tends to overestimate model performance during peak streamflow and to underestimate during low-streamflow conditions [23,59]. To partly overcome this problem, NSE is often calculated with log-transformed observed and modeled values. In this study, both NSE and NSE with log-transformed streamflow values (LnNSE) were thus employed as streamflow efficiency metrics.
where Q sim , Q obs , andQ obs are simulated discharge, observed discharges, and arithmetic mean for observed discharges, respectively.

Critical Success Index (CSI)
We determine the agreement between observed (MODIS) and simulated snow via CSI [60]. It is also known as the Threat Score (TS), combines different aspects of the POD and FAR, describing the overall skill of the simulation relative to reference observation [61]. The value of the CSI score ranges between 0 and 1, with 1 being the ideal score. The formula for CSI is defined as Equation (6), and two-dimensional contingency table fit for evaluating observed and simulated binary events is given in Table 4 [23]. Table 4. Set-up of the 2 × 2 contingency table for binary snow cover data comparison. O and S, respectively, represent observed and simulated binary snow cover and the subscripts refer to a snow-free (0) and snow-covered (1) grid cell.

S.N.
S 1 S 0 Sum O 1 n 11 n x1 n x1 O 0 n 10 n 00 n x0 Sum n 1x n 0x n xx Shyft generates fractional snow cover area (fSCA) as one of its main output variables. Thus, in order to convert the fractional snow cover area into binary values, a threshold fSCA value of 0.4 was used to distinguish between snow and snow free cells.

Intercomparison of the Catchment Discretization Methods
The different catchment discretization methods resulted in different number and distribution of computational cells with in the study catchment. The fixed median area of TIN cells (4 km 2 ) has yielded to lower number of computational cells (671) as compared to the regular square grid (839). The distribution of cells across the elevation range of the study catchment shows that most of the cells fall within the elevation range of 4000 to 6000 m asl. Most of the reductions in the number of computational cells when using TIN were also obtained in this elevation range (Figure 4). TIN cell areas vary from 0.16 to 15.3 km 2 , while in the hypsographic discretization method, zonal areas range from 6.9 to 726.8 km 2 (Table 5).  Remarkable differences are observed in grid/TIN area distributions with respect to aspect. In regular space square grid, a higher percentage of area is facing towards southwest direction, while in TIN a higher percentage of area is facing towards northeast direction ( Figure 5). Differences in area distribution with respect to aspect in different discretization methods ultimately result in differences in snow/glacier cover area in the catchment.   Low flow simulations are relatively well captured by PTGSK-TIN and RPTGSK-TIN as compared to the other models. Modeled discharge during pre-monsoon season from all models was higher than observation, which is caused by the positive biases in the WFDEI precipitation during this season [10]. From the daily simulated hydrograph presented in the Figure 6, it can be seen that there is no notable difference in the simulated discharge from the model with and without using the imputed radiation.

Evaluation of Catchment Discretization and Imputed Radiation on Streamflow Efficiency
Effect of the catchment discretization methods and radiation algorithms on model performance was further evaluated using NSE and LnNSE (Figure 7). The NSE and the LnNSE functions are a pair of conflicting objective functions, where the NSE function evaluates the ability to reproduce all stream flows, but it is known to be biased to predict peak flows [23], while the LnNSE function emphasizes low flows [62]. The model efficiency metrics were computed for each hydrologic year constituting the validation period. As such, each factor had 10 samples with a total size of 60 samples (two radiation algorithms and three discretization types). For model efficiency under high-flow condition as estimated using the raw-data (NSE) (Figure 7a), a more significant difference in model efficiency was observed between the catchment discretization types as compared to the radiation algorithms. The intercomparison between discretized models shows that TIN was superior to both the square grid cells (SqGrid) and the elevation bands (i.e., HYP), while HYP has yielded better result as compared to SqGrid. A similar order of performance between the three catchment discretization types was observed under the two radiation algorithms. Although a slight improvement was observed in SqGrid as a result of using RPTGSK against PTGSK, a slight deterioration in model performance was observed when using HYP and TIN with RPTGSK. A superior result was also obtained in model performance during low-flow conditions (LnNSE) when using TIN as compared to other catchment discretization types (Figure 7b). Unlike to the result obtained using NSE, a relatively better simulation result was obtained under SqGrid than HYP when assessed using LnNSE. No significant difference in LnNSE was also observed between the results obtained from PTGSK and RPTGSK.

Effects of Imputed Radiation on Snow Simulations
Results presented in the previous sections show that RPTGSK-TIN and PTGSK-TIN performed better than the remaining models. Therefore, these models were selected for snow cover analysis and validation of their snow cover simulations through comparison against the MODIS snow product. Figure 8a demonstrates monthly average percentage of snow cover area from MODIS, PTGSK-TIN, and RPTGSK-TIN. MODIS observations and model simulations show higher percentage of snow cover area during winter periods (January-April). The observed time series with strong seasonality in snow cover are adequately represented with both models. However, simulated values tend to underestimate the observed values during the low snow cover period i.e., June-October. To find the agreement between the model simulation and MODIS, monthly correlation coefficients were calculated. Highest monthly correlation coefficients based on five years (2005-2009) simulation were observed between MODIS and RPTGSK-TIN snow during February (0.85), April (0.76), and May (0.63). Similar to RPTGS-TIN, highest correlation coefficients between PTGSK-TIN and MODIS were observed during February (0.75), April (0.64), and May (0.60). Negative correlation coefficients between MODIS and both models were observed during December and January. It is also important to note that the snow simulations from both models at the beginning of 2006 was near to zero. However, we can see such abrupt decrease in snow cover at the beginning of every years. It is due to the fact that in the gamma snow model, snow cover area from the model start accumulating from the beginning of the calendar year. Snow cover areas from the MODIS during July, August, and September are higher than from the models. Higher percentage of snow cover from MODIS might be connected to the fact that MODIS does not distinguish between snow and glacier, and considers both as a snow cover. Unlike to MODIS observations, glacier area is not simulated in Gamma snow, hence lower percentage of snow cover area during the snow ablation periods was found. The other possible source of uncertainty in the snow observations from MODIS during this period is related to cloud cover, which is obviously very high during monsoon period as mention by Riggs et al. [63]. Figure 8b shows five year daily average snow cover area percentage differences between the two models simulations and with the MODIS. With the implementation of imputed radiation in to the models, significant differences were observed in snow simulation. With in the two model snow simulation, highest differences (up to 12%) were observed during the winter season. However, when compared with the MODIS the difference was up to 35%. This highest difference between MODIS and model simulation is due to the fact that MODIS is not able to differentiate between snow and glaciers but from the model it only simulate snow. Observed highest peak snow from models are related with the model forcing datasets particularly precipitation. Because we have used relatively coarse forcing datasets. Therefore, when it interpolate it might give precipitation to the grid/TINs with no precipitation. The highest CSI from both models was obtained during pre-monsoon season (March-June) and the CSI score from both models are nearly the same. Median CSI score for both models during the pre-monsoon was approximately 0.75, indicating that 75% of the simulated snow from PTGSK-TIN and RPTGSK-TIN matches the remote observations. Similar to the correlation coefficient values, the lowest CSI scores were observed during the winter and monsoon periods. There is a slight increase in CSI score for RPTGSK-TIN simulation relative to PTGSK-TIN during June, July, and August. The scatterplots between the 8-day snow cover area fraction from models and MODIS in different elevations are shown in Figure 10. As can be seen from this figure, the topography seems to have a dominating effect on snow cover distribution. In the elevation ranges up to 3000 m asl, very few snow events were detected from both models and almost no snow events were detected from the MODIS. The ranges from 3000 to 6000 m asl experience greater snow coverage, whose spatial distribution is relatively better captured by both data sets. A slightly higher correlation coefficient between MODIS and RPTGSK-TIN simulation was observed as compared to PTGSK-TIN ( Figure 10). During the snow ablation period, the high mountain ranges (>6000 m asl) have permanent glaciers which are not simulated by Gamma snow but detected in MODIS giving lower correlation coefficients for both models.

Discussion
The basic assumption in this study was that different model discretization methods were the reason for the differences in model performance. The results obtained through a variety of statistical models analysis undeniably follow the catchment discretization methods. We found all six models performed well (NSE > 0.75) for both calibration and validation. In general, all discretization methods provide an excellent representation of the general flow pattern and the overall water balance, while maintaining the significant interannual variability satisfying. As higher efficiency discharge simulation from the TIN-based discretized models also means higher efficiency snow cover simulation, only RPTGSK-TIN and PTGSK-TIN are used for snow cover simulation.
The predicted hydrographs for the calibration/validation, shown in Figure 6, confirm that the overall hydrograph pattern is predicted quite well by all models. Similar to the previous studies in [14,64,65], better performance was observed for the models with TINs. The poor performance of the grid-based models are partly attributed to complex topography of the region and the associated difference in distribution of land area ( Figure 11). According to the NSE and the LnNSE of the outlet discharge, the rank of discretized model efficiency, in descending order, is TIN, HYP, and SqGrid models both in the calibration and validation periods. SqGrid are better than HYP when assessed using LnNSE. Theoretically, the ability to account for spatial variability of meteorological forcing and physical features within a catchment should lead to better simulations [66]. However, similar to the prior studies by Reed et al. [67] and Smith et al. [68], we found that grid based distributed modeling approaches do not always provide improved discharge simulation compared to hypsography discretized conceptual models. As presented by Birhanu et al. [69] and Xu et al. [70], land cover types are highly dependent on elevation change, thus the catchment discretized by the elevation bands represent catchment better than the grid. As a result, in our study, better land representation in hypsography-based models results in higher model efficiency compared to grid based model. The Delaunay triangulation is a widely appreciated and investigated mathematical model used to represent the topography and is highly efficient for building TINs. Vivoni [65] mention that in the regions of high terrain variability, the catchment discretized with TINs can be modeled more precisely due to higher flexibility of the mesh. The areas with little variation in topography are represented by larger triangles and the areas with more variation are represented by smaller ones leading to better representation of the ground shape [15]. Figure 11 shows the distribution of land cover area in grid and TIN. Remarkable differences are seen for the glacier and forest cover area. Another reason for finding higher amount of forest area in the TIN-based discretization was because of assuming more vegetation types as a forest. Therefore higher efficiencies from the models with TIN are also attributed to the better representation of land cover area. One of the most useful characteristics of a TIN for hydrologic system is the ability to define stream in terms of triangle boundary segments. This allows a more continuous description of stream paths and networks in conjunction with the topography. By comparison, grid data tend to produce zig-zag meandering paths for stream on upslope portions of a watershed [64]. As indicated by Jenness et al. [71], Raster datasets, such as DEMs and surface-area grids, are inherently less accurate and precise than vector data sets such as TINs and polygons. Grid values do not reflect the actual measurements. The grid model is not adaptive, and although TINs will naturally represent an area with detailed relief information with a denser triangle pattern than area with a smoother relief, grids will be far less flexible to cope with variable levels of details [15].
The overall performance of the CSI estimated by the PTGSK-TIN is somewhat reduced compared to RPTGSK-TIN 5 years (2005-2009) simulation runs for the Marsyangdi_2 catchment as expected (Figures 8 and 9). However, biases in catchment average snow cover area from PTGSK-TIN and RPTGSK-TIN are less then 1%. The central estimate (median) from RPTGSK-TIN is slightly better (Figure 9; July-August-September) than PTGSK-TIN for studied catchment, giving better representation of observed simulated snow cover area.
Significant differences between snow simulation were observed by the implementation of imputed radiation algorithm in to the Shyft. Differences were higher during the winter season (i.e., up to 12%). When we compared with the MODIS snow, the difference between the simulation and the model as compared with MODIS was found to be higher (up to 4%) during summer. This is because MODIS does not distinguish between snow and ice, but models only simulate snow. With the imputed radiation, notable differences were also observed in evapotranspiration (see: Table 6). As seen in Figure 12A, unlike to Figure 12B, radiation is not distributed uniformly, and more radiation was received by south-east and south-west aspect. Lower percentage of the radiation is distributed to the Northern aspect. It should be noted that small percentage of higher radiations are also seen on northern aspect, but similar to the study by Burtscher [72], higher radiation in the northern aspect is because of the high elevation area there. Although the overall efficiencies in terms of CSI and NSE are not significantly different (Figures 6 and 8), effect of implemented imputed radiation can be clearly seen in water balance components ( Table 6). As seen in Figure 12A, a lower percentage of radiation is distributed towards the East-South aspect resulting in lowering the total glacier melt. Furthermore, the lower percentage of evaporation is also observed for the model with imputed radiation and could be caused by reduced radiation.

Conclusions
This study was focused on evaluation of the effect of three catchment discretization approaches: hypsography (HYP), square grid (SqGrid), and TIN on hydrological model simulation results, i.e., streamflow and snow cover. Different discretization methods are not compared for snow cover in this paper. Use of two different algorithms for imputed radiation and its impact on model simulation results was also assessed. The different configurations were successfully applied to discharge simulation in Marsyangdi_2 river catchment. Generally, a good agreement between the observed and modeled discharge was obtained from the different configurations. The model simulation results based on the different discretization approaches, i.e., TIN, HYP, and SqGrid, were evaluated using graphical and certain efficiency metrics and they have shown a decreasing order of performance.
Selection of the catchment discretization methods depends upon the availability of computational resources and acceptable model efficiency. Regular SqGrid models are less flexible, but offer higher speed, lower memory requirements, and easier implementation algorithms as most important assets, making them to be preferred when the studied area is more flat. For area with high relief, typical for hydrological modeling, TINs are a priori the preferred option.
Direct use of radiation from WFDEI is a good option for hydrological modeling in a catchment. Snow simulation results did not show any changes when using imputed radiation; the snow cover determined by the direct use of radiation approach in the high mountain region does not represent the real snow cover area. Additional work is still needed to test and validate the suitability of snow simulation determination at different spatial and temporal scales using imputed radiation. Snow simulation with the models with TIN and imputed radiation should also be validated by other models with a much wider range of terrain types.
Finally, our study area represents a very complex regime in the central Himalayas, Nepal, and the results are encouraging. The findings presented in the paper provide valuable new regional knowledge and contribution in the catchment hydrology of the Himalayan region. The TIN-based models are particularly suited to the Himalayan catchments with steep descent/ascent because of the uniform slope along each triangle facet. With including values of TIN and the bias corrected WFDEI dataset is the best performing methods for the discharge simulation and water resource-related model application in the Himalayan region. This methods can be used in the other Himalayan catchments. If the study is constrained by time and computational memory, hypsography-based discretization methods represent a better option than the grid-based method.