Parsimonious Modeling of Snow Accumulation and Snowmelt Processes in High Mountain Basins

Abstract

The success of hydrological modeling of a high mountain basin depends in most case on the accurate quantification of the snowmelt. However, mathematically modeling snowmelt is not a simple task due to, on one hand, the high number of variables that can be relevant and can change significantly in space and, in the other hand, the low availability of most of them in practical engineering. Therefore, this research proposes to modify the original equation of the classical degree-day model to introduce the spatial and temporal variability of the degree-day factor. To evaluate the effects of the variability in the hydrological modeling and the snowmelt modeling at the cell and hillslope scale. We propose to introduce the spatial and temporal variability of the degree-day factor using maps of radiation indices. These maps consider the position of the sun according to the time of year, solar radiation, insolation, topography and shaded-relief topography. Our priority has been to keep the parsimony of the snowmelt model that can be implemented in high mountain basins with limited observed input. The snowmelt model was included as a new module in the TETIS distributed hydrological model. The results show significant improvements in hydrological modeling in the spring period when the snowmelt is more important. At cell and hillslope scale errors are diminished in the snowpack, improving the representation of the flows and storages that intervene in high mountain basins.

1. Introduction

The success of hydrologic modeling of high mountain basins depends in most cases on the correct quantification of snow accumulation and melting processes. According to [1], snow accumulation in winter as well as spring snowmelt gives to mountain catchments a particular hydrological response that should be taken into account when modeling river runoff. Also, mountain catchments are very sensible to temperature changes; therefore climate change can drastically impact the hydrological cycle [2,3]. For example, snowmelt is one of the processes intervening in the hydrological cycle and interacting with many other processes [4]. According to [5], climate change is likely to impact the seasonality and generation processes of floods, which has direct implications for flood risk assessment, design flood estimation, and hydropower production management. Therefore, it is very evident the importance of modeling the accumulation and snowmelt in mountain basins where a very high percentage of water comes from snow [4]. However, mathematically modeling of these processes is not an easy task due to the high spatial and temporal variability of the snow accumulation by itself and strong relationship of snowmelting and with precipitation, temperature, and orographic effects [6,7]. The study their relationship with runoff coming from snowmelt is a recent and recurrent topic [8,9,10]. These processes and relationships are traditionally modeled with mathematical simplifications that are based on simple approaches such as degree-day models and more complex conceptualizations such as energy balance models [11]. These mathematical simplifications are generally implemented as modules executed in hydrological conceptual rainfall-runoff models. There are recent publications that use models developed to implement the snow-runoff process (i.e., [12,13,14,15]).

Unfortunately, the use of energy balance models is not the practical solution due to the high number of variables that can be relevant and can change significantly in space and, in the other hand, the low availability of most of them in a typical [16]. In addition, if used, most of the needed variables are indirectly estimated that increase the uncertainty in the results [16]. For the above mentioned reasons, we decided to the degree-day approach as an option for mountain catchments with limited observed input. According to [17] the classical degree-day model is the conceptualization most used for its parsimony and easy application in medium and large basins. Use only temperature and a spatially constant parameter denominated degree-day factor in their formulation, comprise one of these simplifications. However, according to Hock [18], these models incorrectly reproduce the spatial variability of the snowmelt, because the degree-day factor actually changes in time and space due to the influence of a series of variables such as: season, ground cover, topography, snow cover, snow contamination, atmospheric conditions and rainfall. Consequently, in recent decades there has been an attempt to introduce the variability of the degree-day factor using hybrid mathematical formulations that consider the distribution of the energy flow of longwave radiation (LWR) and shortwave radiation (SWR) as those proposed by [6,18,19,20,21]. Nevertheless, Hock [22] pointed out that there are two problems when using the distributed degree-day model: (a) despite working with long time periods, its precision diminishes when the temporal resolution increases and; (b) the spatial variability is not modeled with precision, because the snowmelt rates can vary substantially influenced by topographic effects such as hill shade, slope and orientation.

For these reasons, our research introduces the spatial and temporal variability of the degree-day factor, using a new modification of the equation of the classical degree-day method with different approaches that will be compared. Moreover, we evaluate the effects of the variability in the hydrological modeling and the snowmelt modeling at the cell and slope scale in two high mountain basins. In this research, we decided to implement the new snowmelt formulation in the distributed hydrological model TETIS, with physically-based parameters developed by the Research Group of Hydrological and Environmental Modeling of the Polytechnic University of Valencia, Spain. The TETIS model simulates the main processes of the hydrological cycle through a conceptualization of tanks [23]. The TETIS model is a free software available at http://lluvia.dihma.upv.es/ES/software/software.html. This model has been implemented in a high number of basins in Spain [24] and Latin America [25], and also there are applications in France [26], USA [27], UK [28], and China [29].

2. Case Studies

The basins selected for the case studies in this paper are the Carson River and the American River basins located in the Sierra Nevada, between the states of California and Nevada in the United States (longitude 118${}^{\circ }$–124${}^{\circ }$ W and latitude 38${}^{\circ }$–40${}^{\circ }$ N) (Figure 1). These two basins have also been used in the Distributed Hydrologic Model Intercomparison Project-Phase 2 (DMIP2) directed by the National Weather Service of the National Oceanic and Atmospheric Administration [30].

Despite being geographically quite close, the hydrologic regimes of these basins are very different due to the average elevation and the location of the dividing line of the Sierra Nevada [31]. The Carson River basin has an area of 922 km${}^2$ and is has an altitude of between 1548 m and 3408 m. This basin presents a hydrologic regime completely dominated by snow. The American River basin has an area of 886 km${}^2$ and has a lower altitude than that of the Carson, between 281 m and 2630 m. The hydrologic regime of this basin is mixed, influenced by rain and snow. According to NFDC1, GRDN2, and CMEC1 streamflow gauges (Figure 1), the most important runoffs in the basins are registered in spring and summer. The Carson River basin registers a mean annual accumulated precipitation of between 559 mm and 1244 mm, and the American River basin records a mean annual accumulated precipitation of between 813 mm and 1651 mm [32]. Furthermore, the mean annual temperature of the Carson River basin varies between 0 ${}^{\circ }$C to 14 ${}^{\circ }$C, while that of the American River basin is between 3 ${}^{\circ }$C and 18 ${}^{\circ }$C [32].

3. Methodology

To perform a better hydrological and snowmelt modeling of high mountain basins with limited observed input, this paper proposes a methodology that consists of several phases (Figure 2). In hydrological modeling we use the TETIS model, which is a conceptual model with physically-based parameters [23]. The TETIS model is described in detail later. On the other hand, in snowmelt modeling we use the degree-day method. This method has shown good results in a high number of basins which vary in size [33] and is widely used because it requires little information and it is easily adaptable to rainfall-runoff models. The choice of this method is basically due to the shortage of complete information of net radiation, sensible energy, latent energy, soil heat and energy advection. This information is necessary to estimate the snowmelt using the models which consider the energy balance applying the conservation laws [34].

For the foregoing, we propose a modification to the original equation of the degree-day model, to introduce the variability of the degree-day factor. The models are automatically calibrated and are compared with the models of the Distributed Hydrologic Model Intercomparison Project-Phase 2 (DMIP2), of the National Weather Service (NWS) of the National Oceanic and Atmospheric Administration (NOAA) [30].

The information and data used in the modeling have been obtained from the database of the DMIP2. This information is available to the public at the following link http://www.nws.noaa.gov/oh/hrl/dmip/2. The database of DMIP2 includes a Digital Elevation Model (DEM), soils and vegetation data, which are used to perform the basin characterization and calculation of TETIS model parameters. In addition, the hourly observed precipitation data, temperature, streamflow, daily snow water equivalent and snow cover data were used as inputs in the calibration and validation of the models. In the hydrologic modeling of these basins, hourly precipitation and temperatures and instantaneous hourly streamflow data have been used. For more information about this hydro-meteorological information you can consult [30]. The modeling has also used monthly potential evapotranspiration data, daily Snow Water Equivalent data (SWE) and a total of 254 satellite images of snowpack provided by NWS of NOAA, used to validate the spatial distribution of mathematically simulated snowpack.

3.1. Hydrological Model

The TETIS model represents the main processes of the hydrological cycle at different spatial and temporal scales, besides water storage and vertical movement on the slopes, in the riverbeds, and the aquifer (Figure 3). In its conceptual scheme, the vertical flow of water is distributed in six levels or tanks where runoff production and state variables are obtained by a simple hydrological balance in each cell (Figure 3a). The horizontal flow of the water is carried out laterally considering three outcomes: direct runoff, interflow and base flow (Figure 3b). These flows follow the flow directions until they reach the drainage network constituted by gullies and streambeds.

The propagation of the flows in the streambed networks is modeled using the so-called Geomorphological Kinematic Wave which is a simplification of the Saint Venant equations as it disregards inertial and pressure terms. To differentiate the aforementioned flow components, the TETIS model uses threshold areas. For example, when the interflow of a cell where the area drained is superior to the threshold area, the flow becomes superficial and flows through gullies (Figure 3c).

3.2. Snowmelt Model

In snowmelt modeling, the original TETIS model uses the classic degree-day method with homogeneous parameters. This model will be used to evaluate and differentiate the effects of the variability of the degree-day factor in the hydrological and snowmelt modeling. The degree-day snowmelt model with homogeneous parameters (MHO) uses a snow tank ($T_0$) in the conceptual scheme of the TETIS model (Figure 3). $T_0$ tank holds the solid precipitation (i.e., snow) differentiates from liquid precipitation using a base temperature. The snowmelt runoff stored in the $T_0$ tank is quantified using the following equation:

$$M_{\left(i\right)}=\left\{\begin{array}{cc}DD{F}_{1,2}\left[T_{a\left(i\right)}-T_b\right],\hfill & \mathrm{if}{T}_{a\left(i\right)}\ge T_b\mathrm{and}SWE>0\hfill \\ 0,\hfill & \mathrm{if}{T}_{a\left(i\right)}<T_b\hfill \end{array}\right.$$

(1)

where $M_{\left(i\right)}$ is the snowmelt runoff in each cell (mmd${}^{-1}$), $DD{F}_1$ is the homogeneous degree-day factor (mm${}^{\circ }$C${}^{-1}$d${}^{-1}$) which this model uses when there is no rain, $DD{F}_2$ is the degree-day factor, also homogeneous, that the model uses when there is energy input from rain (mm${}^{\circ }$C${}^{-1}$d${}^{-1}$), $T_{a\left(i\right)}$ is the mean air temperature in each cell (°C) and $T_b$ is the base temperature (°C). This model also uses a threshold temperature to differentiate solid from liquid precipitation, which usually is $T_b$. The inputs that the MHO model requires are precipitation and mean air temperature data. The snowmelt rates and SWE in each cell are the outputs given.

4. Spatio-Temporal Variability of the Degree-Day Factor

Our proposed mathematical formulation of snowmelt, identified as MDI, was designed taking into account the models developed by [6,18,19,21,35,36,37,38]. In particular, our proposal conserves the parsimony of the classic model described previously (i.e., it can be used in basins with limited observed input), but introduces the spatial and temporal variability of the degree-day factor using maps of relative radiation indices. We decided to use shortwave radiation (SWR) because according to Kustas et al. [19], this radiation better explains the variability of the snowmelt process. In addition, Ohmura et al. [39] identify shortwave radiation as the greatest source of energy for snowmelt. Therefore, our snowmelt model is different from previous ones, because it considers the spatial (scale of cell) and temporal (monthly) distribution of degree-day factor considering the position of the sun and the topography. The formulation that the MDI model uses is the following:

$$M_{\left(i\right)}=\left\{\begin{array}{cc}\left[M{f}_{1,2}·I_{R\left(i\right)}\right]\left[T_{a\left(i\right)}-T_b\right],\hfill & \mathrm{if}{T}_{a\left(i\right)}\ge T_b\mathrm{and}SWE>0\hfill \\ 0,\hfill & \mathrm{if}{T}_{a\left(i\right)}<T_b\hfill \end{array}\right.$$

(2)

where $M{f}_1$ is the degree-day factor (mm${}^{\circ }$C${}^{-1}$d${}^{-1}$) to be affected in space by $I_R$, $M{f}_2$ factor for rainy days (mm${}^{\circ }$C${}^{-1}$d${}^{-1})$, $T_{a\left(i\right)}$ is the mean air temperature in each cell (°C) at time i, $T_b$ is the base temperature (°C) and $I_{R\left(i\right)}$ is a dimensionless coefficient for time of the year i that introduces the energy values of shortwave radiation actually received in each cell compared with the values of shortwave radiation received on a flat surface:

$$I_{R\left(i\right)}=\frac{SW{R}_{v\left(i\right)}}{SW{R}_{h\left(i\right)}}$$

(3)

where $SW{R}_v$ is shortwave potential radiation (Whm${}^{-2}$) calculated considering the slope, orientation, relief shading and the zenith angle depending on the time of the year; $SW{R}_h$ is shortwave potential radiation received on a flat surface at the same altitude and time of the year (Whm${}^{-2}$). To calculate SWR, a DEM of 400 m resolution and the script Area Solar Radiation of ArcGIS, which uses the Viewshed algorithm developed by [40] and modified by [41], was applied.

Actually, the MDI model has been evaluated in two modalities: (i) MDI-1, using a single $I_R$ map computed as the mean annual value, and (ii) MDI-6, using six monthly $I_{R\left(i\right)}$ maps, calculated for the months from January to July, and applied by symmetry to each month of the year. In Figure 4, there are some examples of the $I_R$ maps obtained for both basins. The $I_R$ calculated for each cell from the SWR are introduced into the model by maps of the radiation indices in ascii format.

5. Calibration and Validation of Models

The hydrologic and snowmelt models have been calibrated automatically employing the optimization algorithm Shuffled Complex Evolution of the University of Arizona (SCE-UA) proposed by [42,43]. The decision to use the SCE-UA has to do with its popularity, being employed by various hydrologic models, such as the SAC-SMA [44], the NWSRFS-SMA [45], the MIKE11/NAM [46] and many others, but mainly because it works perfectly with the effective parameters of the TETIS [23]. Also, automatic calibration allows a fairer comparison between different models. The methodology followed in this research consists of a combined automatic calibration, that is, the parameters of the hydrologic model TETIS and the parameters of the MHO and MDI models have been calibrated together.

Table 1. Runoff production parameters of TETIS model calibrated automatically with the SCE-UA algorithm.

Parameter	Correction Factor	Effective Parameter
Static storage (mm)	$F{C}_1$	$H_{u\left(i\right)}^{*}=F{C}_1·H_u$
Vegetation cover index $(-)$	$F{C}_2$	${\mathsf{\Lambda }}^{*}=F{C}_2·{\lambda }_v$
Infiltration capacity (mmh${}^{-1}$)	$F{C}_3$	$k_{s\left(i\right)}^{*}=F{C}_3·k_s$
Overland flow velocity (ms${}^{-1}$)	$F{C}_4$	$u_{\left(i\right)}^{*}=F{C}_4·u$
Percolation capacity (mmh${}^{-1}$)	$F{C}_5$	$k_{p\left(i\right)}^{*}=F{C}_5·k_p$
Interflow velocity (mmh${}^{-1}$)	$F{C}_6$	$k_{ss\left(i\right)}^{*}=F{C}_6·k_{ss}$
Deep aquifer permeability (mmh${}^{-1}$)	$F{C}_7$	$k_{ps\left(i\right)}^{*}=F{C}_7·k_s$
Connected aquifer permeability (mmh${}^{-1}$)	$F{C}_8$	$k_{sa\left(i\right)}^{*}=F{C}_8·k_{ps}$
River channel velocity (mms${}^{-1}$)	$F{C}_9$	$v_{\left(t\right)}^{*}=F{C}_9·v_{\left(t\right)}$

shows the parameters of the TETIS model calibrated by correction factors (FC) which make a global correction of the initial parameters quantified with the information available [23].

Table 2. Parameters of the snowmelt models calibrated automatically with the SCE-UA algorithm.

Snowmelt Model	Parameter	Effective Parameter
	$DD{F}_1$ (mm${}^{\circ }$C${}^{-1}$d${}^{-1})$	$DD{F}_1^{*}=DD{F}_1$
MHO	$DD{F}_2$ (mm${}^{\circ }$C${}^{-1}$d${}^{-1})$	$DD{F}_2^{*}=DD{F}_2$
	$T_b$${(}^{\circ }$C)	$T_b^{*}=T_b$
	$M{f}_1$ (mm${}^{\circ }$C${}^{-1}$d${}^{-1})$	$DD{F}_{1\left(i\right)}^{*}=M{f}_1·I_{R\left(i\right)}$
MDI	$M{f}_2$ (mm${}^{\circ }$C${}^{-1}$d${}^{-1})$	$DD{F}_{2\left(i\right)}^{*}=M{f}_2·I_{R\left(i\right)}$
	$T_b$${(}^{\circ }$C)	$T_b^{*}=T_b$

shows the parameters of the snowmelt models calibrated automatically with the SCE-UA.

In the evaluation of the results obtained with the hydrological model and the different alternatives for snowmelt models, we use the Nash-Sutcliffe Efficiency index (NSE) as an objective function. This index is commonly used as a measure of efficiency in hydrologic models given that it involves the standardization of the residual variance and its expected value does not change with the length of the record [47,48,49,50]. The NSE contemplates efficiency between $-\infty$ and 1, one being the perfect adjustment between the observed flow and the simulated flow. In scientific literature, a calibration is acceptable when the NSE value is higher than 0.6 and excellent when the NSE value is higher than 0.8 [51]. The NSE calculation has been carried out using the following equation:

$$NSE=1-\frac{{\sum }_{i=1}^n{\left[q_{obs\left(i\right)}-q_{sim\left(i\right)}\right]}^2}{{\sum }_{i=1}^n{\left[q_{obs\left(i\right)}-{\overline{q}}_{obs\left(i\right)}\right]}^2}$$

(4)

where $q_{obs}$ is the observed flow, $q_{sim}$ the simulated flow and ${\overline{q}}_{obs}$ is the mean flow observed.

The calibration consisted of a warm up of a water year from which the initial moisture conditions and state variables used in basin modeling are obtained (

Table 3. Periods used in the warm up, calibration, and validation of the models in the Sierra Nevada basins.

Modeling Phase	GRDN2 (1)	CMEC1 (1)	NFDC1 (2)
Warm Up	Oct./1991–Sep./1992	Oct./1991–Sep./1992	Oct./1991–Sep./1992
Calibration	Oct./1992–Sep./1994		Oct./1992–Sep./1994
Temporal validation	Oct./1994–Sep./2000		Oct./1994–Sep./2000
Spatial validation		Oct./1992–Sep./1994
Spatio-temporal validation		Oct./1994–Sep./2000
Gauges Carson basin $\left(1\right)$ and American basin $\left(2\right)$

). The calibration period selected consisted of two water years, one wet and one dry, the former, with the objective of obtaining effective parameters that were robust and capable of representing these variations in annual events. We decided to use a short period of two years to reduce the computation times employed in the optimization, which when combined with a separated structure of effective parameters of the model can be sufficient.

The predictive capacity of the models in was evaluated through spatial, temporal and spatio-temporal validation (Table 3). In the validation of the snowmelt models implemented, six SNOTEL stations have been used (Figure 1): Blue Lakes, Ebbetts Pass, Poison Flats, Spratt Creek, Blue Canyon and Huysink.

In the evaluation and analysis of the results obtained with the different alternatives for snowmelt models, they have been: the Root of Mean Square Error (RMSE) and the SWE Centroid Date (SCD) has been used as in other studies such as by [52].

The RMSE is widely used and it is computed using the following equation:

$$RMSE={\left[\frac{{\sum }_{i=1}^n{\left[q_{obs\left(i\right)}-q_{sim\left(i\right)}\right]}^2}{n}\right]}^{1/2}$$

(5)

where n is the number of observations.

The SCD is computed using the following equation:

$$SCD=\frac{\sum t_{i}xSW{E}_i}{\sum SW{E}_i}$$

(6)

where $SWE$ is the daily observed or simulated $SWE$ in mm, t is the number of the day from the beginning of snow accumulation and i denotes an individual $SWE$ value.

Finally, simulated snowpacks were validated with a total of 202 satellite images at a resolution of 1 km. Obtaining the differences between the area of the simulated snowpacks and the area obtained from satellite images in percent. Moreover, the DDF variability effects in the hillside snow modeling use a set of cells located on the northern and southern hillsides was selected.

6. Results and Discussion

The first analysis carried out consisted of comparing the effective parameters of the hydrological models obtained with the SCE-UA algorithm. The effective parameters in both basins show significant differences, because the state variables and the flows in each cell have different results influenced by the hydrological regime of these basins (

Table 4. Average effective parameters of the TETIS model and effective parameters of the snowmelt models obtained in the Carson and American River basin.

Effective		Carson Basin			American Basin
Parameter	MHO	MDI-1	MDI-6	MHO	MDI-1	MDI-6
$H_u^{*}$	112.43	125.6	113.31	192.33	187.16	200.07
${\mathsf{\Lambda }}^{*}$	0.77	0.50	0.52	1.38	1.34	1.30
$k_s^{*}$	30.02	43.52	26.51	18.64	18.64	17.70
$u^{*}$	2.94	1.56	0.69	2.74	2.31	2.58
$k_p^{*}$	8.39	19.47	14.47	4.89	6.66	4.80
$k_{ss}^{*}$	16.1 × 10${}^3$	21.6 × 10${}^3$	12.3 × 10${}^3$	5.8 × 10${}^3$	6.3 × 10${}^3$	6.9 × 10${}^3$
$k_{ps}^{*}$	0.89	1.43	1.96	0.0	0.0	0.0
$k_{sa}^{*}$	41.70	301.88	58.01	33.65	85.30	32.30
$v^{*}$	1.36	0.87	1.29	0.99	1.32	1.04
$DD{F}_1^{*}$	3.39	1.2–3.7	0.34–3.90	3.51	0.93–4.8	0.23–3.21
$DD{F}_2^{*}$	3.17	1.3–4.0	0.7–7.3	7.96	1.19–6.2	0.24–3.40
$T_b^{*}$	2.54	1.80	2.91	2.05	1.64	1.17

).

Nonetheless, as can be observed in Table 4, the effective parameters of the snowmelt models do present significant differences directly related to the spatial and temporal variability of the $I_R$ used in each case (Figure 5). Table 4 also demonstrate that in both basins, the MDI-6 model has obtained minimum DDFs and is lower than the fixed range of 1 to 11.6 mm${}^{\circ }$C${}^{-1}$d${}^{-1}$, taking into consideration the results of [22,53,54]. Also, a greater variability is achieved for both basins when six $I_R$ maps are used and follows a pattern clearly influenced by the orientation of the hillside, the relief shading, and the zenith angle. The former coincides with that stated by [22] who considers it fundamental to introduce the effect of these variables in the DDF, which has been achieved with this new snowmelt model.

On carrying out the sensitivity analysis of the MDI-1 and MDI-6 models to the $DD{F}_1$ and $DD{F}_2$ parameters, it was observed that both snowmelt models are not very sensitive to the $DD{F}_2$, because the condition that combines rainfall and air temperatures is higher than the base temperature that activates the snowmelt mechanism and is not representative. Greater sensitivity was observed in the snowmelt models for the $DD{F}_1$ factor.

6.1. Rainfall-Runoff Modeling

The efficiencies obtained in the flow modeling observed in the calibration period has demonstrated that in both basins the DDF variability significantly influences the quantification of the flows at the streamflow gauges. On analyzing separately the efficiencies of a wet year and a dry year, it was observed that there are important differences in the snowmelt models. In

Table 5. Statistical of efficiencies obtained by the models in the simulation of the streamflow for the period of calibration in the Carson and American basin.

Gauges	Statistics		Wet Period			Dry Period
		MHO	MDI-1	MDI-6	MHO	MDI-1	MDI-6
GRDN2	NSE	0.926	0.877	0.900	0.351	0.425	0.737
	RMSE (m${}^3$s${}^{-1}$)	4.492	5.823	5.316	4.787	4.508	3.049
NFDC1	NSE	0.896	0.894	0.890	<0	<0	0.573
	RMSE (m${}^3$s${}^{-1}$)	11.999	12.117	12.314	9.942	8.104	4.994

, it can be observed that the MHO obtained high efficiencies for the wet period and worse results for the dry year in both basins. In the same tables, it can be seen that the MDI-1 also obtained worse results in the dry period, improving its efficiencies in the wet period. The MDI-6 model obtained very good efficiencies in the wet period and the best results in the dry period.

In Figure 6 we can better observe the effect of the day-grade factor on overestimations of the snowmelt, which produce large contributions to the melt flow in winter and decrease considerably in the spring thaw. Although the differences do not seem so significant, it can be seen that the MDI-6 model approaches better to the recession curve of the observed hydrograph (Figure 6). This indicates that introducing the variability of the degree-day factor allows a better approximation to the real snowmelt rates. The opposite case is observed in the models MDI-1 and MHO, where the rates are too high, causing a faster snowmelt even in winter months. This is more evident in the drought period, where the sensitivity of using homogeneous degree-day factor produces very large errors in hydrological modeling (Table 5).

The predictive capacity of the models evaluated by temporal, spatial and spatio-temporal validation is in the range of those achieved by the DMIP2 model [30]. The efficiencies in the temporal validation show that MDI-6 is the model that achieves the best prediction in both basins (

Table 6. Efficiencies obtained in the validations of the streamflow for the Carson and American basin.

Basin	Gauges	Model	NSE	RMSE (m${}^3$s${}^{-1})$
		MHO	0.756	29.975
American	NFDC1 (temporal)	MDI-1	0.743	34.814
		MDI-6	0.786	29.721
		MHO	0.757	8.967
	GRDN2 (temporal)	MDI-1	0.735	9.358
		MDI-6	0.810	7.930
		MHO	0.873	4.437
Carson	CMEC1 (spatial)	MDI-1	0.875	4.399
		MDI-6	0.881	4.289
		MHO	0.664	11.820
	CMEC1 (spatio-temporal)	MDI-1	0.682	11.495
		MDI-6	0.751	10.186

). In the case of the spatial validation, again the MDI-6 is the model that obtained the best efficiencies in the Carson River basin. In the spatio-temporal validation, the results show the best efficiencies distributed between the two models. In general, the MDI models show a prediction of between acceptable and very good where the effect of variability in the modeling of the flows observed is positive with respect to the MHO model.

In continuous simulations (Table 6), the efficiencies of the MDI-6 are also better and of the order of those obtained by models: HL-RDHM [55], NWSRFS [43], TOPKAPI [56] and GR4J [45] in the same basin (

Table 7. Efficiencies obtained by the DMIP2 models and TETIS-snowmelt models, in the period of temporary validation (Oct./1994–Sep./2000).

Models	Carson Basin (GRND2)		American Basin (NFDC1)
	NSE	RMSE (m${}^3$s${}^{-1}$)	NSE	RMSE (m${}^3$s${}^{-1}$)
HL-RDHM	0.91	4.85	0.89	17.92
NWSRFS	0.88	5.68	0.89	17.91
TOPKAPI	0.81	7.01	0.87	19.68
GR4J	0.80	6.57	0.73	28.89
TETIS-MHO	0.77	7.70	0.76	27.12
TETIS-MDI-1	0.75	8.04	0.77	26.87
TETIS-MDI-6	0.85	6.75	0.83	25.93

).

6.2. Modeling Snow Accumulation and Snowmelt

In snow modeling, the point accumulation of the snow is evaluated using daily SWE data from four SNOTEL stations located in the Carson River basin and two in the American River basin. The efficiencies achieved in the calibration period have shown that the models obtained the best results in the SNOTEL stations located at an altitude higher than 2000 m (

Table 8. Efficiencies obtained in the modeling of the accumulation of snow (SWE) in the SNOTEL stations in the calibration period for the Carson and American basin.

Basin	SNOTEL	Model		Calibration		Validation
			NSE	RMSE (mm)	NSE	RMSE (mm)
	Blue Canyon	MHO	0.659	80.899	0.284	99.950
		MDI-1	0.547	93.280	0.210	104.96
American	(elev. 1609 m)	MDI-6	0.273	118.103	0.079	113.313
	Huysink	MHO	0.789	174.596	<0	559.10
		MDI-1	0.847	148.763	0.427	538.85
	(elev. 2011 m)	MDI-6	0.900	126.255	<0	551.976
	Spratt Creek	MHO	<0	119.857	<0	90.600
		MDI-1	0.154	92.157	0.382	57.604
	(elev. 1863 m)	MDI-6	<0	106.279	<0	88.920
	Poison Flats	MHO	0.840	90.459	0.701	115.889
		MDI-1	0.821	95.948	0.691	117.939
Carson	(elev. 2357 m)	MDI-6	0.867	82.858	0.824	88.921
	Blue Lakes	MHO	0.778	175.00	0.885	151.293
		MDI-1	0.844	146.894	0.917	128.082
	(elev. 2455 m)	MDI-6	0.878	129.934	0.936	112.880
	Ebbetts Pass	MHO	0.871	137.388	0.902	165.09
		MDI-1	0.916	110.98	0.892	173.307
	(elev. 2671 m)	MDI-6	0.954	82.306	0.922	147.37

).

With regards to the data on the stations located at lower altitudes, the three snowmelt models overestimated the snow accumulation. The above results coincide with those registered by [57,58] who obtained better results in the SNOTEL stations located at a higher altitudes in the same basins with an increment of uncertainty for those stations at lower altitudes. The results of the MDI-1 and MDI-6 models are acceptable as those obtained in the DMIP2 Project by the SNOW-17 model (HL-RDHM and NWSRFS) and the TOPKAPI model that uses the energy balance [30]. In the calibration period, the best results in the modeling of snow point accumulation were obtained with the MDI-6 model. Therefore, the variability of the DDF positively influences and improves the approximation of the snow point accumulation. The results achieved in the validation period have again shown greater efficiencies in the Carson River basin in the SNOTEL stations located at a higher altitude. However, the biases show underestimations and overestimations of the three snowmelt models in both basins (Table 8). In the American River basin, the three snowmelt models do not adequately reproduce accumulation in the validation period. In general, an underestimation of snow accumulated at the monitor points is seen with respect to the observed data. On evaluating the models by SWE Centroid Data (SCD) it has been observed that the errors generalized in the two snowmelt models and for both basins is a product of the great heterogeneity of snow accumulation (Figure 7). Figure 7 shows that the MDI-6 model is the one that achieves a better approximation with lower SCD indices in the Carson River basin. In the American River basin, the three snowmelt models presented greater SCD indices with a greater uncertainty of the accumulation modeling in the two SNOTEL stations used (Figure 8).

In the case of the spatial distribution of the snow, it was seen that the DDF variability has a very strong effect improving the results positively. That is, the modeled snowpacks have an excellent approximation to the snowpacks obtained from the satellite images when the DDF variability is used (Figure 9). The results also showed that the positive effect is greater when the six $I_R$ maps are used with the MDI-6 model (Figure 9). The former has demonstrated the importance of the energy flow from the SWR and the variables considered in its quantification such as topography, relief shading, and the zenith angle. In the same figure, we show that the MHO that uses homogeneous DDF obtained the worse results.

On analyzing the rates of snowmelt in the cells on the north and south hillside, important differences were also found in the snowmelt models. The greatest differences were recorded in the American River basin in the cells on the north hillside and for the wet period (Figure 10). The effect of the hillside is better represented by the MDI-6 model. The former is influenced by the number of $I_R$ maps used by this model. In the same figures, the MHO is not capable of differentiating the hillside using only the mean air temperatures and the homogeneous DDF. In quantitative terms, it has been estimated that the MHO melts some 19% more than the MDI-6 in the first months of winter in the cells on the north hillside of the Carson River basin. In the American River basin, it has been estimated that the MHO melts some 15% more than the MDI-6 in the first months of winter in the cells on the north hillside.

In the case of mass balance carried out in the selected cells, the results have shown the same behavior as the snowmelt rates in both basins. The MDI-6 model accumulates a greater amount of snow on the south and north hillside than the MDI-1 and MHO models. In the cells on the north hillside of the Carson River basin, an average accumulation of 64 mm more than that modeled with MHO has been estimated. In the American River basin, the MDI-6 accumulated an average of 63 mm more than the MHO in the north hillside cells.

In doing this research, the authors intended mainly that the greater advantage of the model MDI with respect to MHO was the way to simulate the spatial distribution of the snow, according to the observed SWE data and the snowpacks obtained from satellite images. In the first case (Figure 7 and Figure 8), it was interesting to note that the MDI-6 in higher-altitude basins better reproduces the punctual accumulation of the snow observed in the SNOTEL stations. With lower elevations in the American River basin, its snowmelt modeling becomes more complex even for energy balance models (references: [30,57]). In the case of the snowpack modeling, the MDI-6 surpass clearly the MHO (Figure 9), which is why we believe that this is the main advantage of the MDI-6 model with respect to the MHO and the true positive effect to include the spatio-temporal variability of the degree-day factor through the use of the $I_R$ maps.

7. Conclusions

In this research, the spatio-temporal variability of the degree-day factor in snow modeling using the degree-day method has been introduced using radiation index maps. The evaluation carried out has shown that the introduction of the spatial and temporal variability of these factors have positive effects on the hydrological modeling at catchment scale (i.e., at the outlet) and the simulation of the spatial and temporal variability of snow accumulation and snowmelt.

In the first case, the simulation of the discharges at the outlet of the two case studies is significantly improved when maps of radiation indices are used. In the calibration period and for both basins, the best degree-day snowmelt model has been the MDI-6, which uses six monthly maps of radiation indices. The MDI-6 model has also demonstrated a predictive ability superior to that of the MHO model with constant snow parameters, with very reasonable approximations in the temporal, spatial, and spatio-temporal validation. In fact, the results with MDI-6 are in the range of the best reported by [30] for these same basins in the DMIP2 project.

In many hydrological and environmental applications, a good spatial description of snow accumulation and melting is very important. In this sense, the results have shown that the spatio-temporal variability of the degree-day snowmelt factors significantly improve the simulation of snow accumulation when compared to the observations in the SNOTEL stations. In particular, the three distributed snowmelt models achieved good results in the stations located at a higher altitude and a greater uncertainty in those located at lower altitudes. This spatial pattern of the snow model reliability coincides with the results reported by [30,57]. Another very significant positive effect of the spatio-temporal variability of the snowmelt factors was observed in the spatial distribution of the snowpack. In this case, the MDI-6 model reached the best approximations to the snow cover observed by satellite images. These positive spatial results are highly relevant, since previous published research about improvements of the traditional MHO model, such as previous published research similar improvements of the traditional MHO model, such as [6,18,19,20,21], do not report this type of spatial validations.

Finally, it should be noted that the proposed distributed hydrological modeling with variable in space and time degree-day factor, does not imply a considerable additional effort compared with the modeling benefits: the needed radiation index maps can be easily computed from a DEM without increasing the information demand for the model implementation in high mountain basins.