Multi-station calibration strategy for evaluation and sensitivity analysis of the snowmelt runoff model using MODIS satellite images

In this study, the snowmelt runoff model (SRM) was employed to estimate the effect of snow on the surface flow of Aji-Chay basin, northwest Iran. Two calibration techniques were adopted to enhance the calibration. The multi-station calibration (MSC) and single-station calibration (SSC) strategies applied to investigate their effects on the modeling accuracy. The runoff coefficients (cs and cr) were selected as calibration parameters because of their uncertainty in such an extended basin. To determine the most substantial input of the model which is the snowcovered area (SCA) from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor imagery, MOD10A2 images were collected with spatial and temporal resolutions of 500 meters and 8 days, respectively. The results show an average of 15% improvement in the model performance in the MSC strategy from the data period of 2008–2012. Also, an appropriate agreement with physical characteristics of the study area could be seen for the calibration parameters. The contribution of snowmelt in the river flow reaches its peak in April and May, then with increasing temperature, the contribution decreased gradually. Furthermore, analysis of parameters indicates that the SRM is sensitive to recession coefficient and runoff coefficients.


INTRODUCTION
In the mountainous regions, snowmelt in most conditions, particularly in warm seasons, is the primary resource to supply and guarantee freshwater (Han et al. 2019). In recent years, most studies have focused on precipitation-runoff models without considering snowmelt (Wu & Chau 2013;Fu et al. 2020), whereas snowmelt provides about 50% of the mountainous countries' surface water and groundwater (Vafakhah et al. 2015). Therefore, broadening our understanding of the process of snowmelt is the critical factor in such regions. The phenomenon of climate change and global warming could impact the time variation in snowmelt and the consequent amount of runoff. Thus, accurate simulation of snowmelt runoff is essential for planning and managing water resources. Among the natural disasters, floods are the most destructive, causing massive damage to human life, infrastructure, agriculture, and the socioeconomic system (Mosavi et al. 2018). Therefore, accurate estimation of snowmelt runoff can also decrease and control the risk of flooding due to rapid snowmelt (Abudu et al. 2012).
In past decades, hydrologists have searched for suitable ways to simulate snowmelt and its impact on runoff. Therefore, two principal approaches have been introduced, energy-balance and degree-day (Zhang et al. 2014). The energy-balance approach is the most comprehensive method for simulating and assessing the surface flow using energy flux between snow, soil, and air. A significant drawback of this approach is the fact it is data-intensive, as it cannot be implemented in basins that face a lack of This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/). sufficient data (Abudu et al. 2016). Unlike the energy-balance models, degree-day base models are more applicable, especially in data-sparse basins (Rulin et al. 2008). Notwithstanding the success of degree-day snowmelt runoff models, some researchers have argued that the temperature cannot be considered as the only practical factor in the process of melting, and have remarked on the necessity of using the solar radiation and surface albedo for providing the energy of snowmelt as well (Li & Wang 2008;Vafakhah et al. 2015). On the other hand, some other researchers have indicated via their studies that temperature is the most influential element in this process (Zeinivand & Smedt 2010;Saydi et al. 2019). Degree-day approach models are known for their simplicity and have been successfully applied in several studies (e.g., Tekeli et al. 2005;Kult et al. 2014;Zhang et al. 2014).
Snowmelt runoff model (SRM) is a degree-day based model, designed to simulate the impact of snowmelt, where it is a significant proportion of the water supply on watershed daily runoff (Martinec 1975). The SRM is a useful tool for high-precision simulation in mountain basins that provides a valuable opportunity to consider snowmelt as a new water resource. This model was previously used in small European basins, and increasing model recognition led to application of the model in vast basins (areas up to 917,444 km 2 and elevation of up to 8,840 m a.s.l.), even in basins where snowmelt does not play a dominant role in the river regime (Martinec et al. 2008). The SRM has been applied successfully in many basins of the world, more than 25 countries. It is one of the few models where the physical-based parameters and remote sensing data are applied to build model inputs. Due to the requirements of fewer data and remotely sensed images to calculate a snow-covered area (SCA) and hypsometric information, this model has been popular for use in data-sparse basins, mostly in inaccessible high mountainous areas.
The key feature of the SRM is the value of the SCA. Accurate assessment of the SCA is one of the basic operations in the field of water resources management in mountainous areas since most precipitation in those areas falls as snow. Lack or absence of snow monitoring stations in most mountainous basins has caused the use of remote sensing tools that have gained considerable popularity among researchers in recent years owing to being free and producing acceptable results. There has been invented a large number of earth observation satellites (e.g., Landsat, Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS)) for the purpose of environmental monitoring. MODIS free images have gained a high reputation among researchers owing to appropriate diversity in temporal resolution (Ringgs & Hall 2016). In recent years, several studies have been implemented on the performance of these images and compared with the results of studies that used different satellite images to estimate their accuracy (Tekeli et al. 2005;Zhang et al. 2014;Han et al. 2019;Hussainzada et al. 2021). Eight-day maximum MODIS snow products (MOD10A2) have been widely used in recent studies to calculate snowmelt runoff and have indicated good performance in the SRM model (Tahir et al. 2019;Sharma et al. 2020).
In the current study, the SRM, using a degree-day approach with utilization of both ground and satellite data, was used to simulate the Aji-Chay basin runoff contributed by snowmelt. The basin is located in the northwest of Iran. Therefore, two main strategies were applied, multi-station calibration (MSC), in which the basin is divided into several sub-basins that have flowmeters at the outlet, and single-station calibration (SSC). In the MSC strategy, the parameters are calibrated simultaneously using data from several observation stations. Investigations on hydrological models reveal that using the MSC reduces errors and increases modeling accuracy. Observation runoff data acquired from the whole basin outlet during the simulation period in large basins are incapable of representing real procesess of the basin due to heterogeneity of sub-basins. In other words, sub-basins in most large basins have different characteristics; hence, to obtain an accurate result, the behavior of all sub-basins must be reflected in the outlets (Andersen et al. 2001;Saeidifarzad et al. 2014). Most physically based models include a large number of parameters, so the SSC strategy cannot describe spatial variations in large basins. In general, MSC enables the use of all available data in the basin, reducing the uncertainty in results, especially in large basins with higher physical heterogeneities (Bekele & Nicklow 2007;Bai et al. 2017). The results of this assessment endeavor to provide further details about MSC and decrease the uncertainty of hydrological models, as well as helping planners and managers to know the surface water supply of the basin in order to make proper decisions against possible water tension.

Case study and data
The Aji-Chay river basin in the northwest of Iran (East Azerbaijan province) was selected as a pilot case study to perform the empirical phase. Aji-Chay River originates from the hillside of Sabalan Mountain and discharges into the Urmia Lake so that its annual inflow to Urmia Lake is 380 million cubic meters (MCM) (Fazel et al. 2017). Urmia Lake is a saline lake in the northwest of Iran. During two decades, the water level of Urmia Lake decreased from a historical maximum of 1,278 m by 6.5 m (1,271.5 m) in July 2020. During this time the lake lost almost 45% by area and 85% by volume (Sima et al. 2021). The Aji-Chay river basin is the largest basin in East Azerbaijan, and its area is 13.853 km 2 . In the present study, due to the acquisition of valid data, the selected case study basin is the upper area of Markid hydrometric station with a drainage area of 6,913.5 km 2 located between 37°68 0 N and 38°47 0 N and 46°47 0 E and 47°89 0 E ( Figure 1). Basin altitude ranges from 1,467 m to, 3,875 m a.s.l. One of the significant defects of the case study basin is the shortage of ground stations at high altitudes. The daily data of discharge and temperature as ground data are obtained from the East Azerbaijan Regional Water Department. Data from two hydrological stations, Markid and Arzanag, were used. They are the outlet for the whole case study and upper sub-basin, respectively. Temperature data from Sarab metrological station were also imposed into the model. Annual mean temperature and precipitation are 9°C and 230.7 mm, respectively (Barzegar et al. 2016). In this study, all data such as ground measurements and remote sensing data were acquired between 2008 and 2012, and are divided into calibration (2008)(2009)(2010)(2011) and validation (2012) periods, respectively.
Detection and determination of different snow properties using remote sensing data are widely accepted hydrology remote sensing tools since they have introduced novel approaches in acquiring required parameters such as depth and area of snow. In other words, spatial and temporal distributions can be detected by remote sensing data in all basins (Jain et al. 2010;Saydi et al. 2019). In order to calculate the SCA, MODIS images (MOD10A2) with the temporal and spatial resolutions of 8-day and 500 m spatial resolutions, respectively, were used. The normalized difference snow index (NDSI) was applied to distinguish snow from clouds and other types of areas (Hall et al. 1995). In this study, NDSI snow mapping algorithm was used to detect snow cover pixels from the cloud-like atmospheric blockage which has been executed automatically on the MOD10A2 images, because many previous studies have shown NDSI algorithm reduces the impact of cloud cover appropriately (Tekeli et al. 2005;Kult et al. 2014;Sharma et al. 2020). The accuracy of MODIS snow products has been already examined by several studies, all of which have illustrated that there is more than 90% agreement between these products (MOD10 A/Y 1/2) and in-situ measurements under clear sky with about 4 cm depth of snow (e.g., see, Hall et al. 2002;Huang et al. 2011). However, since most days of the year are cloudy, daily images (MOD10A1) might lose their accuracy. Hence, to overcome this issue, 8-day composite (MOD10A2) was developed (Abudu et al. 2016), and in this study, the values obtained from these images (maximum snow extent file) were used as the SCA inputs. The benefit of using MODIS images against other satellite images like Landsat is pixel size (500 Â 500 m 2 ) which is useful in large basins. Extracted SCA from the MODIS data is a more popular approach for use in the SRM, mainly because the MODIS data do not need to be processed to classify snow cover from other features such as lake, cloud, and no-snow area, and have been improved in terms of special resolution and geolocation accuracy by the user (Abudu et al. 2016).
As well, ASTER image is used as digital elevation model (DEM) for measuring the basin's physiographic characteristics. Figure 2 illustrates the daily time series of ground measurements (discharge, rain, and temperature) between 2008 and 2012.

Concept of SRM
Several snowmelt simulation models have been developed to suit particular necessities and hydrologic conditions. The SRM because of its simple concept (degree-day) is a feasible alternative to complex models which use massive data. Most models that can generally handle the various hydrologic conditions are data-intensive and complicated (Tekeli et al. 2005). The SRM, which utilizes the snow cover and meteorological data as inputs, has been the most broadly employed in mountainous basins covered by snow and ice (Abudu et al. 2012). The SRM is a degree-day and physically based model that uses temperature as an advocate of energy flux among snow, soil, and air to simulate the snowmelt process. According to Equation (1), the SRM estimates daily streamflow and discriminates runoff from snowmelt and precipitation as (Martinec & Rango 1986): |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} snowmelt runoff þ c rn P n : A : 0:116(1 À k nþ1 ) |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} precipitation Àrunoff þ (Qs n þ Qr n )k nþ1 |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} runoff contribution from the previous day (1) where Q is the average daily discharge (m 3 /s); c is the snow (c s ) or rain (c r ) runoff coefficient; a is the degree-day factor indicating the melt depth in each increment of the number of degree-days (cm°C À1 d À1 ); T is the number of the degree-day (°C d); DT is the temperature adjustment factor extrapolating the number of degree-day (temperature) from the station to the mean hypsometric elevation of each zone (°C d); S is the ratio of the SCA to the total study area; P is precipitation contributing to runoff (cm); A is the area of basin or elevation zone (km 2 ); n is the sequence of days during the discharge computation period; the coefficient of 0.116 converts data from cm km 2 /day to m 3 /s; and k is the recession coefficient. Generally, models are sensitive to some parameters in the area, thus, the estimation of these parameters' values requires consideration of basin characteristics (e.g., vegetation cover, soil condition, evapotranspiration, geographical location, etc.) and interactions between them. Since this process is time-consuming, developers of the SRM introduced several approaches according to hydrological judgment, empirical and mathematical relations.
Recession coefficient (k) represents the decline of discharge in a period without snowmelt or precipitation. This parameter is estimated with a diagram in which historical data of discharge (Q n and Q nþ1 ) are plotted against each other. X c and Y c are calculated according to Equations (2) and (3): The degree-day factor is estimated according to the empirical Equation (4): where r s is snow density and r w is water density. In order to calculate the value of the degree-day parameter according to Equation (4), limited observational data were recorded for 2 days over 2007-2012. An important feature of SRM is the critical temperature (T crit ) through which precipitation is considered as rain or snow. According to the suggestion of previous studies (Zhang et al. 2014;Adnan et al. 2017), the critical temperature of 0°C is used in the simulation in this study. As well, the temperature lapse rate (g) is the key factor in estimating temperature adjustment (DT). In general, the temperature lapse rate of À0.65°C/100 m is chosen for simulations, but Andaryani et al. (2019) found a lapse rate of À0.6°C/100 m for the Aji-Chay basin. The lag time parameter (L) indicates the time between snowmelt and the appearance of its effect at the runoff measuring station (Abudu et al. 2016). To get an accurate calculation of this parameter, hourly time series of runoff and temperature are needed. There are several proposed values to use as lag time, which determines the ratio of the snowmelt impact on runoff. Here, since the observed data in Aji-Chay basin are at daily scale, this parameter was selected as 18 h based on size, topography of the basin, and analogy with other comparable basins (Martinec et al. 2008).
A zone-wise simulation approach was used in this study. Martinec et al. (2008) recommended that the case studies should be divided into several elevation zones owing to extrapolating temperature and precipitation data. According to the previous studies (e.g., Zhang et al. 2014), the basin was divided into five elevation zones with 500 m increments. Daily temperature and precipitation were extrapolated from Sarab metrological station for five zones using lapse rate (À0.6°C/100 m) (Andaryani et al. 2019) and 2% increase per 100 m (Martinec et al. 2008), respectively. The SCA are extracted for five zones from MOD10A2 images. The most key parameters to have high impact and uncertainty compared to other parameters are runoff coefficients (c r and c s ). Runoff coefficients are related vegetation cover, topography, soil condition and the complex interaction between them, and thus appear to be the primary candidate for calibration (Martinec et al. 2008;Abudu et al. 2012). In this study, these parameters (runoff coefficients (c r and c s )) were adjusted during the process of calibration via both MSC and SSC strategies. Here, a version of SRM, WinSRM (see, Martinec et al. 2008), based on Windows platform was applied in this study.

MSC and SSC strategies
In this study, the MSC strategy was considered to calibrate the SRM. The results of this strategy were compared with SSC strategy results within a large and complex basin. As runoff coefficients are sensitive to vegetation cover and soil condition, they were selected to be calibrated. In hydrologic processes, the main aim of calibration is to achieve feasible values and not just to increase the accuracy. Therefore, MSC strategy can estimate parameters' values that are more rational to local characteristics (Nkiaka et al. 2018). In SSC strategy only data of one station (Markid hydrometric station) was used as observed data to assess performance of the model. However, in MSC strategy, runoff coefficients were examined simultaneously within subbasins, thus, the response of every change in the parameters are recorded in both sub-basins. This process continues until the model achieves an optimal performance. The aforementioned strategies were considered in two sub-basins, then the computed runoff was summed. Finally, to gain correct comparison between the performances of the two strategies, computed runoff was compared with the observed runoff data of the whole basin outlet (i.e., Markid station).

Performance criteria
Specifically, some criteria have been proposed to assess the accuracy of a hydrological model including Nash-Sutcliffe or determination coefficient (R 2 ) and volume difference (D v ) (Martinec et al. 2008). Here, quantitative assessment of model performance was carried out using the aforementioned indexes as Equations (5) and (6): where, in Equation (5), Q i is the measured daily discharge, Q 0 i is the computed daily discharge, and Q is the average measured discharge of the given period; in Equation (6), V R is the measured yearly (or seasonal) runoff volume, V 0 R is the computed yearly (or seasonal) runoff volume. The optimum values of R 2 and D v are 1 and 0, respectively. Positive D v means underestimated, whereas negative overestimated simulations.

Calculation of elevation zones and SCA values
In mountainous regions, most hydrological and geological phenomena such as temperature, hypsometry, and snowpack are affected by altitude (Zhang et al. 2014). The SRM is designed to use different elevation zones to improve the final result. Therefore, in most studies, basins are divided into several zones. In this study, the SRM for the case study was simulated considering five zones with 500 m increments (1,467 m to 3,874 m) (see Figure 3). Table 1 indicates the area of each zone for the case study. More than 3,700 km 2 of the total area (about 54% of the whole area) is in the elevation range of more than 1,800 m, which confirms that the given basin is a mountainous area. Among input variables, the SCA is considerably important, thus, to extract the SCA, MODIS snow images (MOD10A2) were used for the hydrologic years of 2008-2012. To obtain daily values as the input, linear interpolation was used. Figure 4 shows the time series of the proportion of SCA in the study basin, which indicates that the three months of December, January, and February have the highest snow cover during a year. However, this proportion plunges to a low of less than 10% in May and throughout the summer recorded SCA was approximately zero, not counting the highest altitudes.
Through the extracted snow images, the ratio of the probability of existing snow cover in each pixel of the study area was examined for every hydrologic year (2008)(2009)(2010)(2011)(2012). The spatiotemporal pattern in Figure 5 indicates that the high probability of snow cover throughout the year (30% and more) occurred at altitudes above 2,300 meters. However, most areas of the case study experienced less than 70 days (20%) cover of snow. Overall, it is readily apparent that the highest proportion of snow  cover occurred in 2012 with the mean SCA of 34.3%. The maximum snow-covered pixels belong to zone E in the northwest edge of the basin (see, Figure 3 and Table 1), covered by snow about 60% days of a year. Figure 6 shows the monthly mean values of SCA, precipitation, temperature, and measured discharge (Markid outlet) from 2008 to 2012. Except for the three months of December, January, and February, the average value of temperature is more than 0°C. Maximum and minimum values, 28.6°C and À17.9°C, are observed in July and January, respectively. The dominant   form of precipitation between March and November was precipitation and for other months it was snowfall. Average precipitation value within five years is 20.7 mm/month, and this value during April, May, and June is more than average (more than 22 mm). The most cover in snow pertains to the coldest month of the year (i.e., January) (see Figure 6). The lowest value for precipitation occurred in December as well. By increasing the temperature and precipitation, SCA is decreased. Therefore, the amount of runoff registers its peak in May (∼14 m 3 /s). From June to the end of October, except for the highest values (altitude more than 3,300 m), no significant snow cover is observed. It can be seen, according to the trends of precipitation and SCA, that precipitation on snow cover plays an accelerator role in the process of melting since the temperature of raindrops causes melting of shallow snowpack increasing the upper threshold which is not accounted for in the SRM (Singh et al. 1997;Zhang et al. 2014).

Sensitivity analysis results
The target of hydrological models is to simplify the different phases of the hydrological cycle. Therefore, models should be examined accurately. The sensitivity analysis of model parameters or inputs is a useful tool to understand their impact rate and determine those parameters which are very sensitive and have more effect on results. In this process, all parameters are kept constant except one, which is changed and its effect examined on the model results. Figure 7 indicates the deviation of R 2 from the optimum R 2 for different values of the parameters. Among the parameters, recession coefficient (k) and runoff coefficients (c r and c s ) show significant impact on R 2 value, respectively. Siemens et al. (2021) reported that runoff and recession coefficients have significant impacts on the SRM model. Also, Panday et al. (2013) found that simulated hydrograph was strongly sensitive to recession coefficients and temperature lapse rate. These parameters largely depend on seasonal variations, snow characteristics, and runoff behavior and are considered as primary elements having a great impact on the model performance so they should be carefully determined. The model did not share remarkable sensitivity with the other parameters, including degree day (a), critical temperature (T crit ), lapse rate (g), and lag time (T lag ). For instance, although different values from 9 to 27 h were devoted to the lag time, no significant change was observed in the simulation result.

Results of SRM simulation
After providing the SCA information for each elevation zone and other inputs (i.e., temperature and precipitation) of the model, a pre-simulation process was conducted. Table 2 shows the values of parameters that were estimated through theoretical and empirical formulas as well as recommendations from previous studies. Figure 7 | Sensitivity analysis of SRM parameters (i.e., recession coefficient (k), rain and snow runoff coefficients (c r and c s ), degree-day (a), critical temperature (T crit ), lag time (T lag ), and lapse rate (g)) based on R 2 criteria.
Finally, the SRM was calibrated and validated by both MSC and SSC strategies for five hydrologic years (data of 2008-2011 for calibration and 2012 for validation purposes). In the calibration process, parameters (runoff coefficients) were modified to have a better agreement between observation and computed hydrograph. In the MSC strategy, the SRM was run in both subbasins (see Figure 1) simultaneously with the same values to reach an optimum point according to the R 2 values. Average values of corresponding days of calibrated parameters for each zone of upper and lower sub-basins were calculated from 2008 to 2011, and then these values were used for the validation period (2012). According to the calibrated values obtained from the MSC and SSC strategies, the SRM was performed, and the results are shown in Table 3. Figure 8 shows the trend of R 2 criterion variations of lower and upper sub-basins in the calibration period acquired from modifying runoff coefficients (c s and c r ), the purpose of which indicates how variations in one sub-basin affect another. In this way, the optimum points (red points), which are maximum turning points, were found as Pareto front. In each attempt, the intensity of variation in the lower sub-basin is larger than the upper sub-basin due to the difference in the area value, since the runoff coefficient from rain (c r ) for the same amount of precipitation leads to more runoff in bigger sub-basins (see Figure 8(b)). Therefore, the impact of snowmelt on runoff is smoother than precipitation. The snow runoff coefficient attains higher values in higher altitudes because of the high proportion of SCA in three upper zones (altitudes above 2,300 meters) in comparison to the lower zones. Furthermore, lower altitudes have experienced a small proportion of snow cover (see Figure 5), thus, in days with no precipitation, snow runoff coefficient should be given a bigger value in order to achieve good model performance. The number of tests to find the best values of calibration parameters with the maximum value of the R 2 criterion was different. In this method, all attempts were performed by variation of runoff coefficients due to the uncertainty of these parameters in large basins with high spatial and altitude heterogeneities. In the MSC strategy, simultaneous modeling and calibration of parameters for two sub-basins lead to dissemination of the effects and conditions of the upstream sub-basin parameters to downstream and vice versa. Therefore, accurate estimations of parameters for each sub-basin could lead to better understanding of the physical interpretation of the model parameters according to the area and altitude.
Values of runoff coefficients (c s and c r ) that were calculated from the MSC calibration strategy are shown in Figure 9. According to the time series, snow runoff coefficient (c s ) had larger numbers in zones D and E (see Table 1 and Figure 3),  which have a high impact on amount of runoff. This may be because high altitudes are covered with rock and, as a result, less penetration occurs. On the other hand, the runoff coefficient from rain (c r ) plays a dominant role in zones A and B because of lack of snow cover in such areas which also have greater areas compared to the rest of the highlands. Also, there is remarkable variation in values from February to the end of May and its impact on runoff increases in March and May when runoff reaches the peak values (see Figure 10). According to the obtained values of R 2 and D v in Table 3, the SRM has acceptable accuracy in simulation of runoff with the contribution of snowmelt. The MSC strategy leads to better performance than the SSC strategy. As the basin has a large area, in whole basin modeling (SSC), the difference in characteristics and features of each sub-basin was not taken into account. However, by dividing the area into several smaller sub-basins, the final accuracy was increased. In 2009, the R 2 index shows a low accuracy in both strategies, which might be due to the low amount of annual runoff in comparison to others. The total volume of runoff in 2009 was 36.41 MCM, which is very low compared to the measured average volume of five years (112.11 MCM). Moreover, another error source and reduction of modeling accuracy is related to recession coefficient. The runoff value is increased and decreased sharply on some days having different behavior in comparison to most days having almost smooth variation; therefore, due to such different behaviors, SRM was not able to simulate such rapid changes for those days appropriately. Furthermore, having data from only one meteorological station also has a negative impact on the modeling efficiency; since in a large basin, temperature spatially deviating during the day compensates for this shortcoming. Figure 10 shows measured and computed daily discharge values. It is evident that the MSC strategy has better performance to catch peak discharge value than SSC. However, by the SSC strategy, the average computed runoff value of five years is 97.75 MCM, but is 92.61 MCM for the MSC strategy in the same period. The reason for this difference is overestimation of the SSC strategy in the validation period where the computed volume is 143.87 MCM. Looking closer at Figure 10(c) shows that the computed volume of runoff for May until the end of June of 2012 in both MSC and SSC was 34.4 and 60.8 MCM, respectively, whereas the measured volume was 35.8 MCM. A small change in the parameters' values, especially runoff coefficients which depend on vegetation cover and soil condition, can affect results significantly. In order to deal with this issue, by categorizing the basin into several small sub-basins, heterogeneity is considered and reliability of calibration parameters is also improved. Furthermore, simulation of peak flow is a major weakness of most hydrological models such as the Soil and Water Assessment Tool (SWAT) (Andaryani et al. 2019). However, according to Figure 11, the SRM can solve this issue and present a good performance.
Hydrographs of snowmelt runoff without the impact of precipitation and measured runoff are presented in Figure 11. The average volumetric values of snowmelt on runoff of five hydrologic years through MSC and SCC strategies are 36.73 MCM and 33 MCM, which accounted for 43 and 36% of the computed runoff volume, respectively. The highest amounts of snowmelt impact in a year occurred in 2008 (about 60%) and 2012 (about 40%), respectively.
The SRM is very practical for basins that are mainly covered by snow. In such basins, a significant amount of annual discharge belongs to snowmelt. According to Piracha & Majeed (2011) and Sharma et al. (2020), the contribution of snowmelt in total annual discharge is more than 60% in the Himalaya region. In this study, however, about 40% of annual discharge was accounted for by snowmelt. Given that Iran is a semi-arid region, a comprehensive study and modeling of available resources is crucial to supply water for the agricultural sector, generate electricity from hydropower (dam), as well as prevent floods caused by rapid melting. The model results represent that the proposed strategy (MSC) leads to better performance through reducing uncertainties and improving the accuracy. Calibration of SRM parameters through the MSC strategy  tends to improve the model accuracy. Sharma et al. (2020) reported that the SRM cannot simulate peak flow. However, the MSC strategy used in this study overcomes this inability, and also increases the reliability of calibrated parameters in large basins. Studying snowmelt in mountainous regions is an important issue in sustainable water planning for the future. Since Aji-Chay is one of the largest basins in Urmia Lake basin, more attention should be paid to the potential of snowmelt runoff.

CONCLUSIONS
Overall, it can be concluded that the 8-day MODIS images (MOD10A2) are reliable data to determine SCA of the Aji-Chay basin regarding the performance of the model in the calibration and validation steps. In general, the quality of input data, the accuracy of calibrated parameters and non-calibrated parameters are very effective in the performance of hydrological models, and similarly in modeling by SRM, the quality of SCA and temperature data as the two principal inputs as well as the number of altitude classes can affect the model performance significantly. Accurate simulation and forecasting of snowmelt process and its impact on runoff are one of the most important issues of water resources management, which should be considered in various fields such as irrigation, flood control, and hydropower generation. Since most of the northwestern and western basins of Iran, including Aji-Chay basin, have a high potential of experiencing snowy days, the use of several models (such as SRM) to predict the amount of snowmelt runoff is an essential step. Furthermore, estimation of snow cover as an important resource is important for sustainable water management, particularly in this basin which is the largest sub-basin feeding Urmia Lake.
In the study basin, there is no snow monitoring station, hence remote sensing data were used. MOD10A2 snow images with the 8-day temporal resolution were used to extract the SCA information. Also, MSC and SCC strategies were used to calibrate the model. According to the results, the SRM was capable of simulating runoff contributing snowmelt in mountainous basins appropriately. To improve performance of the model, the MSC strategy was used. The MSC strategy considers data and physical information of sub-basins to report more reliable results than the SSC strategy. Modeling accuracy in terms of R 2 values was improved up to an average of 12 and 22% when using the MSC strategy with regard to the SSC strategy in the calibration and validation steps, respectively. Analyzing sensitivity of parameters revealed that recession coefficients and runoff coefficients are more sensitive than other parameters.
Due to the vast area of the basin, data from several meteorological stations to calculate the temperature and lapse rate are recommended to consider in future study. Also, recession coefficient is highly sensitive to the sharp fluctuations of the runoff over short time intervals and it is seen to be estimated using data with fine temporal resolution. These two issues can be considered in future studies to increase the accuracy of modeling, of course, in the case of availability of such data. The degreeday (a) and lag time (T lag ) parameters can also be calibrated using appropriate optimization methods like genetic algorithm, as well as using more accurate experimental methods such as in-situ data collection to increase the parameters' accuracy. To overcome the limitations of optical images in terms of cloud cover, unmanned aerial vehicle-based radar data such as synthetic aperture radar could be useful in future studies.

CONFLICTS OF INTEREST
There is no conflict of interest.

AVAILABILITY OF DATA AND MATERIAL
The data will be available upon request.

CODE AVAILABILITY
The code will be available upon request.

DATA AVAILABILITY STATEMENT
Data cannot be made publicly available; readers should contact the corresponding author for details.