Using Hydrological Models and Remote Sensing to Characterize Snow Pack Behaviour in High Mountain Environments

1. Introduction

Sparse meteorological and snowpack data pose challenges for hydrological modelling in the Upper Indus Basin. The amount of water stored in snowpacks in the Upper Indus Basin (UIB), known as snow water equivalent (SWE), and its utilization in hydrological applications is constrained due to significant uncertainty for water resource management and downstream floodings [1,2,3]. Snow plays a crucial role in high altitude settings regulating resources such as water supply, hydropower and flood hazard [4,5]. The topographic complexity of montane terrain increases the spatial variability in snow mass and snow extent. This is influenced by a myriad of factors including elevation, slope, aspect, topography, orographic effects, vegetation coverage, wind redistribution, aerodynamic roughness atmospheric circulation [6,7]. Given the importance and the dynamic nature of mountain snowpacks, accurately estimating the SWE stored within these regions is an important scientific and societal [8,9].

In recent decades hydrological models have risen as potent tools for investigating SWE [10,11]. They are versatile because they primarily use widely available meteorological data requiring site-specific data for calibration and validation. Hydrological modelling has its greatest challenges where data, including meteorological inputs are scarce, and especially in remote montane regions where snowpack data are exceptionally rare. Utilizing global satellite and re-analysis datasets that furnish gridded meteorological and SWE information might be a promising avenue to enhance the assessment of hydrological processes and simulate snowpack behaviour [1]. Calibrating hydrological models by satellite and reanalysis measurements also offers scope to improve the limitations of traditional calibration based solely on observed streamflow [12,13]. Incorporating various variables, both individually and in combination with streamflow measurements, offers potential to improve model performance and enrich understanding of hydrological processes [14]. However, the challenge in data-scarce regions is the absence of quality ground observations which instead is typified by low-quality, including missing, irregularly sampled, or poorly validated data. One promising avenue involves leveraging emerging earth observation technologies [15].

Hydrological models span a range from data-intensive energy balance models to temperature index models [16]. Previous investigations that have compared temperature-index and energy-balance models have yielded inconclusive results regarding the superiority of one approach over the other [17]. The energy balance model presents a favourable alternative due to its reduced calibration demands. However, the model performance is contingent upon the availability and accuracy of the supplementary climate input data it requires [16]. The temperature-index model benefits from readily available input data, although its calibration for individual catchments demands substantial labour [18]. In high altitude and high latitude settings, the Soil and Water Assessment Tool (SWAT) operates as a temperature-index model, relying on air temperature and precipitation as its primary inputs to replicate the dynamics of snow accumulation and snowmelt processes.

In mountain environments securing reliable observed climate data is challenging due to limited spatial coverage of field stations and because global gridded precipitation products are prone to significant local errors. This creates problems for monitoring or forecasting the risk of hydrometeorological hazards like flooding and landslides [19]. However, remote sensing and reanalysis precipitation datasets hold significance in these regions where data availability is limited [20]. Precipitation stands as the primary contributor to uncertainty which leads to varying optimal ranges for the fine-tuned parameters [21]. The elevation band method employed in SWAT addresses the challenge of having limited precipitation stations within a catchment [21]. However, this method can result in errors when estimating precipitation data for specific sub-catchments, especially in those with varied topography. The most promising alternative involves incorporating worldwide data, such as high-resolution daily precipitation datasets like ERA5-Land (European Center for Medium-Range Weather Forecasts fifth generation), satellite precipitation data from CMORPH (Climate Prediction Center Morphing Technique), CPC UPP (Climate Prediction Center Unified Precipitation Project), APHRODITE (Asian Precipitation Highly Resolved Observational Data), and TRMM (Tropical Rainfall Measuring Mission). Satellite and reanalysis precipitation products require evaluation against in-situ data and calibration to hydrological models before use, with varying regional impacts on product performance [22,23].

The available tools for SWAT parameter calibration, sensitivity analysis, and uncertainty analysis are either commercial products or free tools with limited functionalities. SWAT-CUP is one of the most widely used software to auto-calibrate SWAT model. SWAT-CUP is commercial software, and its free version has certain limitations, such as restrictions on the maximum number of simulations and parallel simulation options. ArcSWAT and SWATEditor also provide options for running SWAT with modified parameter values, but they lack automatic parameter calibration, sensitivity analysis, and uncertainty analysis techniques. In contrast, several free or open-source tools have been developed specifically for SWAT+, which is a restructured version of SWAT [24]. These tools include the SWAT+ toolbox and the Integrated Parameter Estimation and Uncertainty Analysis Tool Plus (IPEAT+) [25]. The development of such tools for SWAT model is essential given the widespread use of this modelling platform. Nevertheless, tasks such as parameter calibration, sensitivity analysis, and uncertainty analysis continue to pose significant challenges. R-SWAT provides a valuable and versatile solution for overcoming the challenges in SWAT parameter calibration, sensitivity analysis, and uncertainty analysis. Its interactive interface and compatibility with R’s wide range of functions and packages make it a powerful tool for advancing hydrological research and modelling [26]. R-SWAT facilitates a deeper understanding of hydrological processes by leveraging the power of open-source SWAT and R. One valuable aspect of utilizing R-SWAT is its capability to provide real-time visualizations of adjustments in simulated and observed flow, snowfall, snowmelt, and other water balance elements at the catchment level when modifying hydrological and snow parameters. This specific feature of R-SWAT enhances its effectiveness for the manual calibration of snow parameters, substantially reducing the labour-intensive nature of the manual calibration process in the SWAT model.

Calibrating hydrological models involves various considerations, inclusive of the algorithms used, the parameters requiring calibration and whether they need calibrating individually or collectively [27,28]. Any model calibrated solely on streamflow data is unlikely to yield dependable spatio-temporal portrayal of other hydrological fluxes and states [29]. The calibration procedure can differ in multiple aspects, including the time step (monthly, daily, etc.), calibration using single vs multiple variables, spatial scale (single point or multi-point calibration), auto- and manual calibration, the choice of objective function employed for optimizing the variable of interest. Multivariate parameter estimation can lead to more realistic internal model dynamics and associated hydrological characteristics [30,31], ultimately enhancing the overall representation of catchment behaviour. Prior research has demonstrated the effectiveness of using remotely sensed data to implement and calibrate hydrological models [32,33,34]. The issue of model equifinality, where multiple parameter sets can produce equally good model performance, introduces uncertainty that can undermine the reliability of hydrological models in forecasting future water resources, especially in mountainous regions where snow processes significantly influence the hydrological cycle [19]. Multi-variable calibration schemes have been proposed as useful means to address model equifinality [35,36] and reduce prediction uncertainty [37]. However, results from multivariate calibration approach have enhanced the power of model simulation in some cases [38], while others have reported no significant improvements or even substantial deteriorations [39,40,41]. There could be multiple reasons to these contradictory findings such as differences in model structures, parameterization, wide variety of remote sensing datasets or case study specific. Nonetheless, multivariate parameter estimation using remote sensing data remains a promising approach, especially in situations where streamflow data availability is limited, or data quality is questionable.

The choice of calibration scheme can significantly influence the outcomes of hydrological models especially in mountainous catchments with sparse input data [27]. Multi-variable (in addition to observed streamflow) calibration is commonly used to improve streamflow simulations in SWAT model, but snow reference data are rarely available due to collection challenges. This was the case in the study catchment necessitating examination of different calibration methods, encompassing single-, multi-variable, manual, and auto-calibration of both lumped and spatially-varied snow parameters. The primary focus of this study is to comprehend the estimation of SWE during the accumulation season. This emphasis is driven by the critical role of SWE and an accurate characterization of SWE is a fundamental requirement for a precise representation of the available water supply. The novelty of this research lies in incorporating snow reference data to manually calibrate snow parameters within R-SWAT to estimate SWE at basin and sub-catchment scales in the Chenab River Basin (CRB). In previous studies, auto-calibration using SWAT-CUP has been widely adopted foregoing arduous manual calibration. It is worth highlighting that the incorporation of SWE data into the calibration process offers a direct means to achieve optimal parameter values with efficiency.

2. Materials and Methods

2.1. Study area

The Chenab River Basin (CRB) is an eastern tributary of the Indus River (Figure 1), originating from the Kulu and Kangra districts of Himachal Pradesh, India. It is formed by the confluence of two major tributaries, the Bhaga, and the Chandra and flows through Indian-administered Jammu and Kashmir before entering the plains of the Punjab province in Pakistan, covering a length of 974 km [42]. The south and north of the CRB are located in the foothills and very high mountains, respectively of the western Himalayas. Approximately half of the total water supply comes from the eastern Hindukush, Karakoram, and western Himalayas [43]. The catchment area of the CRB, upstream of the Marala Barrage is approximately 26258 km². The elevation in the upper snowy catchment of the Chenab River peaks at 7,085 m.

The CRB is situated in an active monsoon belt, experiencing monsoon rainfall between mid-June to mid-October. The character and timing of precipitation vary according to location and altitude. Singh et al. (1995) observed that precipitation mainly occurs during the monsoon (about 75%) and pre-monsoon, while 15-26% occurs during winter months as snowfall in both the Middle and Greater Himalayan ranges. In the outer Himalayan ranges, winter precipitation falls as rainfall due to lower altitudes.

A significant proportion of the flow in the Chenab River derives from snowmelt during the mid to later summer season, augmented by pre-monsoon and monsoon rainfall. Approximately 84% of the annual discharge occurs between April and September, especially from June onwards due to synchronized flow from snowmelt and precipitation [44].

2.2. Data Sets Employed

Data used to configure, calibrate and run SWAT across the CRB spanned the years 1998 to 2015. The model requires topographic, soils and land cover inputs along with climate data, inclusive of precipitation patterns along with minimum and maximum temperatures, relative humidity, solar radiation, and wind speed. Streamflow records from the basin outlet and 19 sub-catchments offers scope to validate simulated runoff and elucidate the role of the snowpack in hydrological response. Accessing such data at scale and over extended timescales is always challenging, especially in mountainous environments [19]. Table 1 provides a summary of the data used in this study, encompassing their attributes and sources.

In hydrological modelling, land use and land cover (LULC) play a crucial role in influencing runoff generation. The LULC map of the CRB reveals that the northern higher-altitude region of the catchment is predominantly characterized by rangeland (covering 10,619 km², which is approximately 40.4% of the area) and permanent snow/ice (constituting 10.70% of the total area, equivalent to 2,827 km², as depicted in Figure 2). Meanwhile, the middle and lower sections of the catchment comprise crops (covering 1,082 km², or 4.12%), trees (5,889 km², approximately 22.4%), barren land (5,130 km², or about 19.5%), and built-up areas (occupying 547 km², approximately 2.1%).

Herein the reference SWE datasets such as High Mountain Asia (HMA) Daily Snow Reanalysis by University of California (UCLA) [45] and modelled snow data [46], hereafter referred to them as UCLA-SWE and KRA-SWE, respectively (Figure 2). These SWE products were used to calibrate snow parameters in the SWAT model and whilst they do not constitute direct observations, they are the best available alternative for in-situ SWE measurements otherwise absent in the basin.

A substantial portion of the CRB is enveloped by glaciers or semi-permanent snow (Figure 2). As a result, it becomes crucial to differentiate between snow that is seasonal and snow that persists year-round on land or glacier surfaces. Potentially erroneous pixels (non-seasonal snow and ice pixels) are masked out including pixels where complete snow melt has not occurred in the UCLA-SWE product (Figure 3b). UCLA-SWE suffers from the absence of validation data across the HMA but has been shown to be effective in tests elsewhere. In both the Sierra Nevada (USA) and the central Andes (South America) satellite-derived snow cover and SWE reanalysis methods yielded positive outcomes relative to observed data [45].

The KRA-SWE product across HMA was developed by Kraaijenbrink et al. [46](Figure 3c). They conducted a comprehensive assessment of SWE and snowmelt at high spatial resolution in all major river basins within the HMA region by integrating a regional-scale snow model with gridded climate data and satellite snow cover observations. The snow model operates at a resolution of 0.05° (~5.7 km), allowing for detailed SWE quantification. To validate the model’s accuracy, they compared the modelled SWE with in-situ measurements of SWE in the Nepalese Himalaya and daily snow-depth maps derived from spaceborne Sentinel-1 radar data.

2.3. Hydrological Model

In mountainous regions precipitation takes the form of both snow and rain as well as in combination. The SWAT model operates as a temperature-based degree-day index model, and it effectively segregates precipitation into rain and snow categories based on temperature thresholds and melting rates. This segregation is fundamental for modelling the intricate interplay between precipitation types and their contribution to streamflow dynamics. Capturing the accretion/depletion behaviour of seasonal snowpacks is also pivotal in determining the snow melt contribution to streamflow.

The setup of the SWAT model involved the utilization of open-source data, including Digital Elevation Models (DEM), Land Use and Land Cover (LULC) information, soil data, and precipitation datasets (Table 1). To consider the spatial variations within the catchment, the basin was divided into sub-basins by the SWAT model, which are further disaggregated into Hydrological Response Units (HRUs). Within the CRB, the subdivision was carried out into 19 sub-basins, each delineated based on a drainage area threshold of 1,000 ha and categorized into three slope classes: 0–3%, 3–10%, and 10–100%. Further division of these sub-basins resulted in a total of 92 HRUs using the SWAT 2012 interface within the ArcGIS 10.5 environment. Each HRU is characterized by unique combinations of soil, land use, and slope. In mountainous basins, where precipitation and temperature exhibit significant orographic effects, SWAT subdivides HRUs into various elevation bands to extrapolate the meteorological forcing point. These elevation bands in SWAT help in accounting for the topographical influences of temperature and precipitation on the simulation of snowpack and snowmelt dynamics. Model simulations were executed at daily time scale, spanning from 1998 to 2015. Notably, the initial 2 years (1998–1999) were designated as a warm-up period, primarily aimed at mitigating the influence of user-estimated initial state variables on the outcomes.

2.4. SWAT snow module

The SWAT model can be implemented as a distributed model with parameters defined at various spatial levels, such as HRU and sub-basins. SWAT calculates the SWE at HRU scale and averaged at each sub-basin. The SWAT inputs are stored in ASCII file format with different file extensions like “. snw", ".hru," ".sub," and ".bsn," representing distinct spatial modelling units and/or properties. The SWE is stored in the ‘.hru’ and in ‘.snw’ file with specific elevation bands.

The model categorizes precipitation as rain or snow based on a user-specified threshold value of the mean daily air temperature. Snowfall accumulates as a snowpack at the ground surface, and the water stored in the snowpack is represented as SWE. The snowpack can increase with snowfall and accumulation or decrease through snowmelt or sublimation. Snowmelt and rainfall contribute to the computation of runoff and percolation. The snow mass balance simulations in the SWAT model are significantly influenced by the snow parameters SMFMX and SMFMN, representing the maximum and minimum snowmelt rates, respectively. As elucidated in the SWAT Input/Output files [47], SMFMX is the maximum melt factor for regions in the Northern Hemisphere and functions as the minimum melt factor for regions in the Southern Hemisphere. In the geographical context of this study SMFMX is anticipated to be greater than SMFMN.

The snowmelt is computed as a linear function of the difference between the average snowpack maximum air temperature and the base temperature, or threshold temperature, for snowmelt. Sublimation from the snow surface is determined based on the function of potential evapotranspiration. Additionally, a temperature lag factor is incorporated to account for the thermal characteristics of the snowpack, which influences the snowpack density, depth, and exposure. As this temperature lag factor approaches 1.0, the mean air temperature on the current day has an increasingly greater influence on the snowpack temperature, while the snowpack temperature from the previous day becomes less influential.

The snowmelt algorithm employed in SWAT is based on a simple temperature index method. The primary climate data required for this algorithm includes maximum and minimum temperature, as well as precipitation. When the mean daily air temperature falls below the snowfall temperature (SFTMP), the precipitation is considered as snow, and its liquid water equivalent is added to the snowpack. The snowpack’s size increases with additional snowfall but decreases due to snowmelt or sublimation.

$${SNO}_i={SNO}_{i-1}+R_{sfi}-E_{subi}-{SNO}_{mlti-1}$$

(1)

Where:

${SNO}_i$ and ${SNO}_{i-1}$ represents the snow water equivalent of snowpack on current day $i$ (mm H₂O) and on the previous day $i-1$ (mm H₂O), respectively.

$R_{sfi}$ represent the amount of precipitation in terms of snowfall on the day $i$ (mm H₂O)

$E_{subi}$ represents sublimation loss during the day $i$ (mm H₂O)

${SNO}_{mlti-1}$ represents the snowmelt loss from the previous day $i-1$ (mm H₂O)

The snowfall temperature (SFTMP) serves as a benchmark for distinguishing between snowfall and rainfall. If the average daily air temperature falls below the SFTMP, precipitation labelled as snowfall, and it will contribute to the accumulation of the snowpack.

The calculation of snow sublimation depends on the potential evapotranspiration within the evapotranspiration module (Neitsch et al., 2011). The snowpack remains in its solid state until its temperature surpasses a specific threshold referred to as the snowmelt base temperature (SMTMP) in degrees Celsius. The computation of the snowpack temperature is carried out according to the following equation:

$$T_{snowi}=T_{snowi-1}\left(1-TIMP\right)+T_{avi}.TIMP$$

(2)

where:

$T_{snowi}$ and $T_{snowi-1}$ are the snowpack temperature on the current day $i$ and the previous day $i-1$, respectively, ${TIMP and T}_{avi}$ are the snow temperature lag factor and mean air temperature on the current day $i$, respectively.

Once the snowpack’s temperature rises above the designated snowmelt base temperature (SMTMP), the process of melting initiates. The amount of snowmelt can be calculated by following equation:

$${SNO}_{mlti}=b_{mlti}.{sno}_{covi}.[\frac{T_{snowi}+T_{maxi}}{2}-\mathrm{SMTMP}]$$

(3)

where:

$$b_{mlti}=\frac{(SMFMX+SMFMN)}{2}+\frac{(SMFMX-SMFMN)}{2}.\sin \left[\frac{2\pi }{365}\right(i-81)]$$

(4)

Here, ${SNO}_{mlti}$ represents the quantity of snowmelt on day $i$ in millimeters of water (mm H2O). The parameter $b_{mlti}$ represents the melt factor for that specific day $i$ in units of millimeters of water per degree Celsius per day (mm H2O °C⁻¹ d⁻¹). The variable ${sno}_{covi}$ corresponds to the proportion of the Hydrologic Response Unit (HRU) area covered by snow. $T_{maxi}$ denotes the highest air temperature on the given day. SMFMX represents the melt factor for June 21st in millimeters of water per degree Celsius per day (mm H2O °C⁻¹ d⁻¹), while SMFMN stands for the melt factor for December 21st in the same units.

2.5. Model evaluation and calibration

The Root Mean Square Error (RMSE), serving as a measure of dissimilarity, along with the Nash–Sutcliffe Efficiency (NSE) and R², recognized as prominent similarity metrics, percent bias (PBIAS) measuring the average tendency of the simulated values larger or smaller than measured data and Kalling-Gupta Efficiency (KGE), which is a hybrid statistical metric, were employed to compare the observed and simulated streamflow within the suggested scenarios. The mathematical formulations for these metrics are expressed as follows:

$$PBIAS=100\mathrm{*}\left[sum\right(sim-obs)/sum(obs\left)\right]$$

(5)

$$RMSE=\sqrt{\frac{\sum _{i=1}^n(Q_i^{obs}-Q_i^{sim})2}{n}}$$

(6)

$$\mathrm{KGE}={\left(r_{-1}\right)}^2+{\left(1-\frac{u_g}{u_{ObS}}\right)}^2+\left(\frac{{\sigma }_g∕u_g}{{\sigma }_{Obs}∕u_{obs}}\right)$$

(7)

$$NSE=1-\frac{⟨{\left(P_1-P_2\right)}^2⟩}{⟨{\left(P_1-P_2\right)}^2⟩}$$

(8)

$$R^2=\frac{{\left[{\sum }_{i=1}^n\left(O_i-\overline{O}\right)\left(s_i-\overline{s}\right)\right]}^2}{{\sum }_{i=1}^n\left(O_{i=0}^{-}\right){\Sigma }_{i=1}^n{\left(s_i-\overline{s}\right)}^2}$$

(9)

2.6. Calibration schemes

Commonly, hydrological models undergo calibration through the utilization of recorded streamflow data collected at hydrometric gauging stations. This approach is effective for refining the model’s performance, but because complex models like SWAT have many parameters that can be optimised depending on the resources and expertise of any given user group – it is possible to generate plausible results but with incorrect behavioural representation. It is thus essential to confirm that calibration through parameterising the model is scientifically robust and gives meaningful and representative model outputs. Additionally, traditional calibration methods could lead to a higher dependence of internal processes, such as the snow module, on the outlet’s streamflow. To address these limitations, a two-stage calibration strategy was implemented. This approach, coupled with the integration of multi-source information, aimed to diminish the reliance of the snow process on the outlet’s streamflow, thus improving the overall constraint of the hydrological model. In this context, SWE products were introduced as novel references for calibrating the SWAT model under the proposed scenarios.

We employed two calibration methods to evaluate the performance of the SWAT model in simulating streamflow and SWE within the mountainous CRB. In the initial approach, the SWAT was calibrated solely using observed streamflow data. In the second calibration approach, SWE reference data were incorporated alongside streamflow data for calibration.

The snow temperature index method in the SWAT model involves several parameters, such as snowpack/snowmelt factors (SMFMX, SMFMN), snowfall/snowmelt temperature thresholds (SFTMP, SMTMP), and temperature lag factor (TIMP). These parameters were manually fine-tuned until achieving acceptable evaluation statistics (NSE, KGE, R² and PBIAS), aiming to optimize the best match between simulated and observed streamflow. Following calibration strategies were applied:

• Initially, a catchment-scale automatic calibration approach was employed with default snow parameters using single variable (streamflow) to optimize all model parameters. This scenario, often employed as the conventional method for calibrating hydrological models in numerous studies.
• Following the process of auto-calibration, manual calibration of snow parameters was performed by leveraging SWE and streamflow data together to attain the best possible outcomes.
• Subsequently, the calibration was performed at a sub-basin-scale, recognizing the need to consider the spatially varied snow parameters.

Drawing on guidance from the literature a total 16 parameters were selected (Table 2) for calibrating the SWAT model in the study basin [20,48]. The selected parameters were divided into three groups comprising snow parameters, catchment parameters and elevation lapse rate parameters. It is important to emphasize that the elevation lapse rate parameters - temperature and precipitation lapse rates (Tlapse and Plapse), can influence the partition of precipitation into rain and snow, as well as the distribution of snow in high-elevation areas. To address orographic controls in mountainous catchments, which affect both temperature and precipitation, SWAT permits the definition of up to 10 elevation bands within each sub-basin. Prior research has indicated that, typically between three and five elevation bands can achieve the accuracy required for most model predictions while maintaining low computational demands, such as running time [49]. Here we followed the upper limit by using five elevation bands to simulate streamflow and SWE in the CRB.

3. Results

3.1. Data evaluation

The area weighted averaged satellite, gridded and reanalysis precipitation datasets in the CRB were evaluated against the available observed precipitation station data (Muzaffarabad and Garhi Dupatta) located just outside the catchment boundary before using it as forcing data to derive SWAT model (Table 3).

All satellite and reanalysis precipitation datasets demonstrate very high mean squared error (MSE) values when compared with observed precipitation data at Muzaffarabad (station-01) and Garhi Dupatta (station-02). Monthly precipitation range (310 to 1189 mm) of CMORPH is much higher than all other precipitation datasets. MSE values of CMORPH are the highest while, ERA5-Land showed the lowest values against observed precipitation. ERA5-Land outperformed all other datasets on the bases of MSE. On the basis of the PBIAS metric, all datasets showed worst performance except ERA5-Land. TRMM and CMORPH shown very high negative and positive bias respectively against both observed precipitation stations as shown in the Figure 4. CMORPH was highly overestimated (Figure 4b) while all remaining precipitation products underestimated precipitation when compared with observed station data (Figure 4a). However, all precipitation products are negatively biased against observed except CMORPH. A cross-validation procedure involving five precipitation datasets over the CRB demonstrated that ERA5-Land provides better estimates of precipitation amounts and variability at monthly to annual scales.

3.2. Global and regional SWE products

Global atmospheric reanalysis products provide an alternative approach for estimating large-scale SWE estimates as by-products generated by their land surface models (Figure 5). These products encompass the Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004) and reanalysis datasets from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5-Land (Hersbach et al., 2020). There remains significant uncertainty regarding how SWE is represented within these datasets, as evident from the notable variations in simulated SWE. Figure 5a depict the comparison of global and regional (e.g., UCLA and KRA) SWE products and Figure 5b illustrate the snow cover area (SCA) at monthly scale in the CRB. The SCA is relevant because the utilization of satellite data to estimate SCA through visible and near-infrared bands indirectly relates to snow mass (Painter et al., 2009).

There is a notable difference in SWE amounts at the catchment level between satellite-derived, modelled and reanalysis products. Such disparities could stem from biases present in atmospheric input data or disparities in model resolution. The ERA5-Land reanalysis product notably underestimated SWE when compared to the higher resolution satellite-derived UCLA-SWE product. GLDAS, a reanalysis SWE product exhibits some similarity with UCLA-SWE between June and December while, it underestimates SWE during the accumulation period (January to April). Similarly, both modelled and satellite derived products also demonstrate similarities during snow accumulation, but they diverge considerably during the snowmelt period. The KRA-SWE product diverges notably from other datasets both in terms of quantity and pattern during the snow accumulation and snowmelt periods.

The KRA-SWE dataset distinguishes itself from other datasets, displaying remarkable disparities in both quantity and pattern during both snow accumulation and snowmelt periods. When examining SWE products in the CRB alongside SCA on a monthly time scale, distinct differences emerge. For instance, the KRA-SWE dataset indicates the accumulation of snow during June and July, contrasting with the decreasing trend observed in the SCA dataset. In contrast, the UCLA-SWE product shows snow accumulation occurring after September, while the SCA product depicts snow accumulation commencing after August. Similarly, both the GLDAS and ERA5 products exhibit even greater discrepancies when compared with the SCA product at the monthly scale.

Figure 6 shows the spatial distributions of the muti-year mean daily KRA-SWE and UCLA-SWE products. Both SWE products substantially differ in spatial distribution of SWE at sub-basin level in the CRB. In the KRA-SWE and UCLA-SWE dataset, sub-basins 14, 15, and 16 exhibited no discernible SWE. In other sub-basins, significant disparities emerge when examining the SWE representation at the sub-basin scale. Over the period from 2000-2015, the KRA (152 mm) and UCLA (109 mm) SWE products exhibit considerable discrepancies at the basin scale which extended to the sub-basin level also.

In this regard, the most notable differences occurring between Sub-basins 10 and 17, where KRA reports overestimation of SWE values of 198.6 mm and 238.5 mm respectively, while UCLA reports underestimation of SWE only 1.3 mm and 62.3 mm, respectively. Similar substantial differences are observed in Sub-basins 2, 5, 18, and 19, where KRA-SWE records values of 296 mm, 271 mm, 275 mm, and 293.3 mm, contrasting with UCLA-SWE values of 182.3 mm, 120.2 mm, 188 mm, and 220.2 mm, respectively. However, it’s worth noting that in sub-basins 3, 4, and 11, both KRA-SWE (180.8 mm, 117.8 mm, and 204.5 mm, respectively) and UCLA-SWE (182 mm, 122.3 mm, and 208.1 mm, respectively) depict SWE with a higher degree of similarity.

3.3. Pre-calibrated SWAT model simulations

The SWAT model was initially run with default parameters and input data from precipitation products as depicted in Figure 6. It is evident that CMORPH and TRMM precipitation inputs did not yield adequate discharge records in the CRB (Figure 6a and 6b, respectively). The former overestimated streamflow while the latter underestimated streamflow. Baseflows and peak flows are not well captured by simulated streamflow with input data from CPCUPP (Figure 6c). However, SWAT model with input data from APHRODITE captures baseflow very well. Only simulated streamflow with input data from ERA5-Land captures both baseflow and peakflow satisfactorily. In general, the ERA5-Land precipitation data demonstrate superior performance when it comes to driving the SWAT model with default parameters. SWAT model simulations derived with other than ERA5-Land precipitation products achieved unsatisfactory results (NS ≤ 0.50).

SWAT model run with APHRODITE, CPCUPP and ERA5-Land precipitation data was further evaluated for its performance in simulation of SWE in the CRB as depicted in Figure 7. These simulations were then compared with two existing SWE products (KRA-SWE and UCLA-SWE). SWE output from SWAT model run with ERA5-Land precipitation data exhibited a better performance when compared with reference snow datasets (KRA-SWE and UCLA-SWE). The simulated SWE showed higher peaks than UCLA-SWE but lower than the KRA-SWE product. However, the SWE simulated by the SWAT model using APHRODITE and CPCUPP precipitation data did not yield satisfactory results. As a result, we chose to include only the SWAT model run with ERA5-Land precipitation data in the subsequent analysis, as the other precipitation products were excluded due to their inadequate performance in simulating streamflow and SWE in the CRB.

Figure 7. Model performance with various precipitation datasets for daily streamflow during 2000-2015. Streamflow simulation is conducted with default SWAT parameters by input data from the satellite and reanalysis precipitation products such as (a) ERA5-Land (b) CPCUPP (c) APHRODITE (d) TRMM (e) CMORPH.

Figure 7. Model performance with various precipitation datasets for daily streamflow during 2000-2015. Streamflow simulation is conducted with default SWAT parameters by input data from the satellite and reanalysis precipitation products such as (a) ERA5-Land (b) CPCUPP (c) APHRODITE (d) TRMM (e) CMORPH.

Figure 8. Daily snow water equivalent (SWE) simulation with SWAT model run with precipitation data from APHRODITE, CPCUPP and ERA5-Land and comparison with SWE products (UCLA and KRA) in the CRB. TRMM and CMORPH data is not shown here due to their poor performance.

Figure 8. Daily snow water equivalent (SWE) simulation with SWAT model run with precipitation data from APHRODITE, CPCUPP and ERA5-Land and comparison with SWE products (UCLA and KRA) in the CRB. TRMM and CMORPH data is not shown here due to their poor performance.

3.4. Single Variable Calibration Approach

After thoroughly investigating the impact of precipitation data on SWAT model simulations (e.g., streamflow and SWE), the best performing precipitation data (ERA5-Land) was used to calibrate and validate the SWAT model. First, all hydrological parameters were calibrated with autocalibration in R-SWAT using observed streamflow in the CRB. After getting the satisfactory ranges of hydrological parameter, SWAT snow parameters were calibrated manually. Table 4 presents a synopsis of the goodness-of-fit statistical outcomes.

The model performance is evaluated with NSE, R² and PBIAS. For calibration period, NSE, R² and PBIAS were 0.63, 0.78, and 3.2 % while during validation period NSE, R² and PBIAS were 0.71, 0.83 and 3.9 % respectively. The match-up between the simulated and observed data suggests that the model is well-calibrated and that peak flows are well captured during both wet and dry years. There is also promising evidence in the validation period with respect to base flows.

Figure 9. Daily streamflow simulation using single variable calibration approach during calibration (2000 to 2011) and validation (2011-2015) in the CRB.

Figure 9. Daily streamflow simulation using single variable calibration approach during calibration (2000 to 2011) and validation (2011-2015) in the CRB.

The simulated SWE (SWAT-SWE) output from calibrated SWAT model was compared with snow reference data (KRA-SWE and UCLA-SWE) in the CRB as shown in Figure 10. SWAT-SWE consistently underestimated peaks in SWE when compared to KRA-SWE during calibration and validation period. However, this trend was irregular when compared to UCLA-SWE. For instance, SWAT-SWE peaks were overestimated during 2000, 2003, 2005, 2007, 2009, 2011, 2012 while, underestimated during 2001 and 2004. The SWAT-SWE showed late accumulation and early snow melt when compared with KRA-SWE. In contrast, SWAT-SWE consistently showed early accumulation and late snowmelt when compared with UCLA-SWE. SWAT-SWE aligned more closely with KRA-SWE during some years (2000, 2007) while, with UCLA-SWE product during certain years (2002, 2008, and 2010). In general, SWAT-SWE peaks remained lower compared to KRA-SWE while, higher than UCLA-SWE in the CRB.

3.5. Multi-Variable Calibration Approach

Figure 11 illustrates the streamflow simulations conducted with reference SWE data (KRA-SWE and UCLA-SWE). The simulated streamflow results using KRA-SWE product to calibrate SWAT model closely resemble those results from a single-variable calibration approach, whereas those utilizing the UCLA-SWE product perform less satisfactorily. Calibrating snow parameters with UCLA-SWE negatively impacted both the timing and magnitude of peak flows (Figure 11a) weakening streamflow simulations. However, the multi-calibration techniques improved the simulated streamflow performance using KR-SWE as opposed to the UCLA-SWE as indicated by the statistical indices provided in Table 4.

Figure 11 illustrates the SWAT-SWE outputs conducted with reference SWE data (KRA-SWE and UCLA-SWE) showing that the SWAT model performed better when calibrated with KRA-SWE as opposed to UCLA-SWE data. The utilization of KRA-SWE for snow parameters calibration significantly improved the simulation of SWAT-SWE. Specifically, NSE (0.94), R² (0.96), and KGE (0.86) values confirm that SWAT model performance was enhanced during calibration and the validation (NSE=0.94, KGE=0.87, and R²=0.97). Using KRA-SWE data, SWAT-SWE adeptly captures the peak values does less well with timings, especially at the onset of snow accumulation and end of the snow melting period between 2000 and 2009. Thereafter SWAT-SWE more closely matches the KRA-SWE and overall the performance of the SWAT model is with the KRA-SWE product as compared to the UCLA-SWE data.

The SWAT-SWE with UCLA-SWE data did not produce as good results as with KRA-SWE data. However, performance of SWAT model to simulate SWE was enhanced during validation (NSE=0.87, KGE=0.91, and R²=0.88) period compared to calibration (NSE=0.64, KGE=0.78, and R²=0.65). Except some years (2008 & 2010), peaks of SWAT-SWE remained higher during calibration period. However, SWAT-SWE captures peaks and base in UCLA-SWE quite accurately.

Figure 12. Comparison of simulated SWE (SWAT-SWE) with UCLA-SWE and KRA-SWE in the CRB. In panel (a), simulated SWE using UCLA-SWE product during calibration and validation period while panel (b), simulated SWE using KRA SWE product during calibration and validation period are shown.

Figure 12. Comparison of simulated SWE (SWAT-SWE) with UCLA-SWE and KRA-SWE in the CRB. In panel (a), simulated SWE using UCLA-SWE product during calibration and validation period while panel (b), simulated SWE using KRA SWE product during calibration and validation period are shown.

3.7. Snow Water Equivalent simulations (Sub-basins)

The comparison of SWAT simulated SWE (SWAT-SWE) with KRA-SWE and UCLA-SWE is presented in Table 5 through various statistical indices. The SWAT model demonstrated exceptional performance to simulate SWE using KRA-SWE as reference snow data in sub-basins BSN01, BSN02, BSN05, BSN09, and BSN12, as indicated by PBIAS<±10. Likewise, sub-basins BSN01, BSN02, BSN05, BSN12, BSN13, and BSN17 exhibited commendable agreement between KRA-SWE and simulated SWE (NSE>0.5). Based on R² (R²>0.60), SWAT simulated streamflow exhibited good agreement with KRA-SWE in all sub-basins except sub-basins BSN06 and BSN11. In fact, the SWAT performance was deemed satisfactory in the majority of sub-basins (BSN01, BSN02, BSN09, BSN12, BSN13, BSN17, BSN18, and BSN19) based on evaluation statistics.

Using UCLA-SWE as reference snow data, the SWAT-SWE simulations achieved satisfactory results in nearly all sub-basins, except for sub-basin BSN13 and BSN10, as evidenced by R² (>0.5). Additionally, based on NSE (>0.5) and KGE (>0.4), the SWAT model displayed a good performance in simulating SWE across sub-basins (BSN01, BSN02, BSN03, BSN04, BSN07, BSN08, and BSN11. However, when assessed using PBIAS (>±15%) evaluation metric, the SWAT-SWE consistently displayed a considerable bias against UCLA-SWE across all sub-basins.

Figure 13 presents a visual depiction of SWAT-SWE with KRA-SWE and UCLA-SWE in selected sub-basins (1a, 2b, 5c, 12d, 13f, and 17f). This graphical representation reveals that, during certain years, the SWAT-SWE closely aligns with KRA-SWE and UCLA-SWE in these sub-basins. However, for most sub-basins, the SWAT-SWE consistently exceeded the KRA-SWE while, underestimated UCLA-SWE. The SWAT model tended to overestimate both snow accumulation and melting, except for sub-basin 17f, where close agreement in SWE peaks was observed. The timing of snow accumulation and melting was also well-replicated in SWAT-SWE with KRA-SWE across various sub-basins, demonstrating a good match at sub-basin level.

4. Discussion

Modelling hydrology in snowy mountain regions requires knowledge of snowpack behaviour and its connectivity to streamflow to validate models with wider spatial applications. However, the case remains that such data are either extremely rare or absent all together placing emphasis on conceptually sound models and making virtue of new data sources including proxies for calibration and validation. Here we used the well described and widely employed SWAT model to estimate SWE incorporating information from remotely sensed and modelled SWE products. The results show that SWE predictions can be improved in SWAT model using reference snow data in general and enhance the overall prediction reliability.

The current discourse in hydrological modelling and calibration involves two opposing approaches: simultaneous calibration and sequential calibration. The simultaneous approach (Brighenti et al., 2019), emphasizes calibrating discharge measurements and distributed snow observations together to complement each other and improve overall model accuracy. In contrast, the sequential approach [28,47], suggests calibrating these components step by step to achieve a more detailed calibration avoiding hidden errors. We adopted the latter - optimizing snow parameters independently before calibrating other hydrological parameters. This approach aims for a systematic and thorough calibration process, enhancing the understanding of the model’s behaviour and potentially improving accuracy by addressing errors in a controlled manner.

This study revealed that satellite-derived UCLA-SWE performed less well than the modelled snow KRA-SWE product in estimating SWE using SWAT in the CRB. Calibrating snow parameters using UCLA-SWE weakened the simulation of streamflow. A key reason for the lack of accuracy of this data lies in the fact that the dataset doesn’t account for the substantial coverage of permanent snow and ice, which makes up approximately 10.70% of the CRB and is entirely omitted in the UCLA-SWE dataset. The accuracy of satellite derived SWE products can vary due to multiple factors, including large biases in the satellite derived products in mountainous regions, low reliability of meteorological data, the geographic and temporal scales considered, cloud cover, masking areas of permanent snow/ice and a lack of ground-truth data for validation in the studied region. Consequently, this permanent snow and ice cover area doesn’t contribute to the catchment’s overall flow, leading to discrepancies in flow simulations using the SWAT model when calibration was performed with UCLA-SWE data.

When employing a single variable calibration approach that relies solely on observed streamflow, the SWAT model demonstrated good agreement between simulated and observed streamflow, but it faced challenges in accurately modeling snowpack. The implementation of a multivariate calibration method, which incorporated snow reference data from KRA-SWE, resulted in enhanced predictions for both SWE and streamflow. However, it’s important to note that UCLA-SWE performed well in simulating SWE in certain individual sub-basins (Figure 13) as compared to the performance at the basin-wide level. This observation suggests that UCLA-SWE data could be valuable in sub-basins where a relatively smaller area is masked out due to permanent snow/ice in the product. This study underscores the effectiveness of the SWAT model in estimating SWE when utilizing a multivariate calibration approach.

The calibration methods detailed in this study have improved streamflow and SWE prediction by incorporating snow reference data, offering potential for more effective operational strategies in forecasting. However, the observation data about snowpack and streamflow may have inherent errors or uncertainties in the measurements. These uncertainties can make it challenging to validate hydrological models accurately. Additionally, when you want to demonstrate that improving SWE simulation results in better predictions of streamflow, these uncertainties in the data can hinder the ability to establish a clear cause-and-effect relationship between the two. In other words, it may be challenging to prove that SWE improvements directly lead to better streamflow predictions because the initial data used for validation contains uncertainties. The results of employing a multivariate calibration approach have strengthened modelling accuracy in specific cases, as noted by [38]. However, in other instances, researchers have reported either minimal improvements or even substantial declines [12,39,41]. These models may not fully capture the complexities of snowpack dynamics, such as the spatial distribution of snow, variations in snow density, and the effects of forest cover. Even if the SWE simulation is improved, the overall model may still have limitations in representing the hydrological processes [50]. The simulation of streamflow is influenced by various hydrological processes beyond solely snowmelt. These processes include groundwater inputs, soil properties, evapotranspiration, alterations in land use, and more. Even with enhancements in SWE simulation, the interplay with these additional processes may introduce uncertainties in streamflow forecasts. Similarly, external factors such as climate change, land use changes, and infrastructure development have a substantial impact on streamflow, but these factors are frequently underrepresented in hydrological models, making it difficult to enhance streamflow simulation by focusing solely on improving SWE simulation. The calibrated parameters used for streamflow simulation may help compensate for SWE discrepancies and potential deteriorations, such as earlier peak flows, if SWE is enhanced by the inclusion of reference snow data [51].

Various approaches to measuring snowpack come with their own advantages and limitations. In the case of in situ monitoring for SWE in mountainous regions, it can yield reliable data but may fall short in capturing local variations due to limited network coverage or placement primarily in valley bottoms. Remote sensing, on the other hand, provides extensive spatial coverage but lacks certain snow parameters, contains errors, and offers coarse resolutions. Hydrological models, while capable of compensating for these shortcomings, introduce errors due to numerical simplifications and a limited understanding of snow physics. To address these challenges, it is recommended to combine these three techniques, enhancing our ability to monitor, estimate, and predict SWE and associated hydrological processes. This integrated approach leverages the strengths of each method, thereby improving the accuracy of hydrological models. Remote sensing, with its capability for large-scale spatial observations, presents opportunities for more comprehensive land surface modeling, especially when integrated with ground observations.

This study brought attention to the challenges associated with using SWE data retrieved from regional products to calibrate snow parameters in a semi-distributed hydrological model to simulate snowpack and discharge. In particular, forcing data from reanalysis product resulted in more realistic simulation of SWE (even with uncalibrated SWAT model) which, in turn, supports the model’s utility for simulating SWE. This underscores the significance of accurate precipitation data and modelled SWE data, especially in regions where snowpack dynamics play a crucial role in hydrological processes. The ongoing advancement of these datasets and modelling frameworks holds great promise. It contributes to enhancing our comprehension of water resources in regions with limited data availability such as the CRB, especially as they grapple with the challenges posed by climate change.

5. Conclusions

Accurate understanding of the Snow Water Equivalent (SWE) is an essential prerequisite for predicting the availability of water in a catchment during the warm season. This study shows that SWAT can accurately simulate SWE and streamflow in snow dominated catchments using snow reference data to calibrate the snow parameters especially in high-altitude catchments. The capability of SWAT model was primarily reliant on the use of precipitation as primary input for simulating streamflow. In this research, the SWAT model underwent calibration using two different methods: single-variable calibration based on streamflow data and multi-variable calibration using streamflow and SWE data. In the single-variable calibration approach, the SWAT model successfully replicated streamflow according to statistical criteria, but the simulated SWE did not match with any of the available reference SWE products. Conversely, when employing a multi-variable calibration approach with observed streamflow and KRA-SWE as the reference snow data, both simulated streamflow and SWE demonstrated good agreement with observed streamflow and reference KRA-SWE. Significant enhancements in the accuracy of simulated SWE were achieved, and there were also improvements, to a certain extent, in the simulation of streamflow when fine-tuning the snow parameters of the SWAT model with KRA-SWE data. However, incorporating UCLA-SWE for calibrating snow parameters resulted in improved SWE simulation but led to a decline in the accuracy of simulated streamflow. The poor performance of the UCLA-SWE in accurately representing streamflow can be ascribed to the fact that UCLA-SWE obscures a substantial portion of permanent snow/ice in the CRB. These results imply that using modelled snow data to fine-tune the snow parameters within the SWAT model, for simulating both SWE and streamflow, may be the preferable method in non-seasonal snow/ice-dominated catchments, especially in cases where remote sensing products might not provide adequate information regarding spatial distributions of SWE estimates. Furthermore, these findings suggest that integrating data from reliable snow reference products to calibrate the snow parameters in SWAT has the potential to accurately simulate SWE in terms of its timing and magnitude. The optimization of parameters notably improved the performance of the SWAT model in estimating SWE.

The primary aim in refining our snow modelling approach was to enhance the accuracy of SWE representation by integrating each SWE product individually. By adjusting snow parameters, we were able to achieve a good result in estimating SWE during accumulation and melting phases, along with slight improvements in model performance concerning discharge. However, considerable differences were observed at the sub-basin level during the snow accumulation and depletion periods. These disparities at the sub-basin level are probably attributable to the intricate topographical and climatic factors characteristic of this region.

Multivariate approach permitted a comprehensive evaluation of snow modelling across the catchment and throughout the entire period. This suggests that model showed its potential for application in the SWE and streamflow simulation to advance water supply forecasting for the future changing conditions. Furthermore, it emphasized the significance of incorporating satellite, modelled, and reanalysis data, for estimating model parameters necessary for simulating streamflow and SWE in a mountainous watershed. The good agreement between the SWAT simulated SWE and KRA-SWE underscores the importance of a multivariate calibration approach, especially in alpine regions with limited ground observations. It holds relevance for ongoing scientific research focused on estimating SWE using remote sensing and re-analysis data. This work also assesses the reliability of recently established reference dataset (UCLA-SWE) for snowpack in a snow-dominated catchment within the Himalayan region. Overall, the integration of remote sensing data into computational model enhances the accuracy of SWE simulations, particularly at the watershed level. This approach can provide valuable information for water resource management, flood prediction, and understanding the hydrological behaviour of a given region. The methodologies detailed in this study can be readily applied to other catchments which often lack comprehensive meteorological data for the assessment of hydrological processes. The implications of this study extend to operational water management in regions with limited meteorological data in high-altitude catchments.