Combining Hydrological Models and Remote Sensing to Characterize Snowpack Dynamics in High Mountains

Abstract

Seasonal snowpacks, characterized by their snow water equivalent (SWE), can play a major role in the hydrological cycle of montane environments with months of snow accretion followed by episodes of melt controlling flood risk and water resource availability downstream. Quantifying the temporal and spatial patterns of snowpack accumulation and its subsequent melt and runoff is an internationally significant challenge, particularly within mountainous regions featuring complex terrain with limited or absent observational data. Here we report a new approach to snowpack characterization using open-source global satellite and modelled data products (precipitation and SWE) greatly enhancing the utility of the widely used Soil and Water Assessment Tool (SWAT). The paper focusses on the c. 23,000 km² Chenab river basin (CRB) in the headwaters of the Indus Basin, globally important because of its large and growing population and increasing water insecurity due to climate change. We used five area-weighted averaged satellite, gridded and reanalysis precipitation datasets: ERA5-Land, CMORPH, TRMM, APHRODITE and CPC UPP. As well as comparison to local weather station data, these were used in SWAT to model streamflow for evaluation against observed streamflow at the basin outlet. ERA5-Land data provided the best streamflow match-ups and was used to infer snowpack (SWE) dynamics at basin and sub-basin scales. Snow reference data were derived from remote sensing and modelled SWE re-analysis products: ULCA-SWE and KRA-SWE, respectively. Beyond conventional auto-calibration and single-variable approaches we undertook multi-variable calibration using R-SWAT to manually adjust snow parameters alongside observed streamflow data. Characterization of basin-wide patterns of snowpack build-up and melt (SWE dynamics) were greatly strengthened using KRA-SWE data accompanied by improved streamflow simulation in sub-basins dominated by seasonal snow cover. UCLA-SWE data also improved SWE estimations using R-SWAT but weakened the performance of simulated streamflow due to under capture of seasonal runoff from permanent snow/ice fields in the CRB. This research highlights the utility and value of remote sensing and modelling data to drive better understanding of snowpack dynamics and their contribution to runoff in the absence of in situ snowpack data in high-altitude environments. An improved understanding of snow-bound water is vital in natural hazard risk assessment and in better managing worldwide water resources in the populous downstream regions of mountain-fed large rivers under threat from climate change.

1. Introduction

Sparse meteorological and snowpack data pose challenges for hydrological modelling in the Upper Indus Basin (UIB). The amount of water stored in snowpacks in the 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 challenge [8,9].

In recent decades, hydrological models have risen as potent tools for investigating the 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 the SWE information is 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 Centre 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 poor model performance or even substantial deteriorations [39,40,41]. There could be multiple reasons for 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 26,258 km². The elevation in the upper snowy catchment of the Chenab River peaks at 7085 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.

Table 1. Data sets used in this hydrological modelling study.

Data Type	Sources	Spatial Resolution	Time Span	Region	Reference
Precipitation	CMORPH TRMM CPCUPP APHRODITE ERA5-Land	~8 km ~25 km ~55 km ~25 km ~9 km	1998–2015	Global	https://www.climateengine.org/
Streamflow	Observed	Basin and sub-catchment	2000–2015	Local	Pakistan Water and Power Development Authority
Min and max temperature, relative humidity, solar radiation, wind speed	Climate Forecast Reanalysis	~38 km	2000–2015	Global	https://climatedataguide.ucar.edu/climate-data/climate-forecast-system-reanalysis-cfsr Accessed on 22 May 2023
Digital Elevation Model	Open Topography SRTM	30 m	-	Global Global	https://portal.opentopography.org https://srtm.csi.cgiar.org
Land use and land cover	ESA CCI LC	300 m	-	Global	https://www.esa-landcover-cci.org/ Accessed on 12 July 2023
Soil	FAO Digital Soil Map	1 km	-	Global	http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/faounesco-soil-map-of-the-world/en/ Accessed on 12 July 2023
Snow Water Equivalent	High Mountain Asia Snow Reanalysis (UCLA SWE)	500 m	2000–2014	High Mountain Asia	https://nsidc.org/data/hma_sr_d/versions/1 Accessed on 15 July 2023
	Daily Snow Modelling (KRA SWE)	5 km	2000–2015	High Mountain Asia	https://zenodo.org/record/4715953 Accessed on 15 July 2023
	ERA5-Land	25 km	2000–2014	Global	https://www.climateengine.org/ Accessed on 15 July 2023
	GLDAS	25 km	2000–2010	Global	https://www.climateengine.org/ Accessed on 15 July 2023

provides a summary of the data used in this study, encompassing their attributes and sources. 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].

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 2827 km², as depicted in Figure 2). Meanwhile, the middle and lower sections of the catchment comprise crops (covering 1082 km², or 4.12%), trees (5889 km², approximately 22.4%), barren land (5130 km², or about 19.5%), and built-up areas (occupying 547 km², approximately 2.1%).

For reference SWE datasets we employed the High Mountain Asia (HMA) daily snow reanalysis by University of California (UCLA) [45] and modelled snow data (KRA) [46], hereafter referred to as UCLA-SWE and KRA-SWE, respectively (Figure 3). The KRA-SWE product across HMA was developed by Kraaijenbrink et al. [46] (Figure 3c). They conducted a comprehensive assessment of the 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.

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].

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.

2.3. Methods

2.3.1. Hydrological Model

The implementation of SWAT hydrological model to simulate streamflow and the SWE is summarized in Figure 4. 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; Figure 4). 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 1000 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.3.2. SWAT Snow Module

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 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 was 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 was incorporated to account for the thermal characteristics of the snowpack, which influences the snowpack density, depth, and exposure. As this temperature lag factor approached 1.0, the mean air temperature had an increasingly greater influence on the snowpack temperature, while the snowpack temperature from the previous day becomes less influential.

The snowmelt algorithm employed in the 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 is 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. 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 \mathrm{a}\mathrm{n}\mathrm{d} 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{S}\mathrm{M}\mathrm{T}\mathrm{M}\mathrm{P}]$$

(3)

where:

$$b_{mlti}=\frac{(\mathrm{S}\mathrm{M}\mathrm{F}\mathrm{M}\mathrm{X}+\mathrm{S}\mathrm{M}\mathrm{F}\mathrm{M}\mathrm{N})}{2}+\frac{(\mathrm{S}\mathrm{M}\mathrm{F}\mathrm{M}\mathrm{X}-\mathrm{S}\mathrm{M}\mathrm{F}\mathrm{M}\mathrm{N})}{2}.\mathrm{s}\mathrm{i}\mathrm{n}\left[\frac{2\pi }{365}\right(i-81\left)\right]$$

(4)

Here, ${SNO}_{mlti}$ represents the quantity of snowmelt on day $i$ in millimeters of water (mm H₂O). 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 H₂O °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 21 June in millimetres of water per degree Celsius per day (mm H₂O °C⁻¹ d⁻¹), while SMFMN stands for the melt factor for 21 December in the same units.

2.3.3. 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:

$$\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}=100\ast \left[sum\right(sim-obs)/sum(obs\left)\right]$$

(5)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{\sum _{i=1}^n(Q_i^{obs}-Q_i^{sim})2}{n}}$$

(6)

$$\mathrm{K}\mathrm{G}\mathrm{E}={\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)

$$\mathrm{N}\mathrm{S}\mathrm{E}=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.4. 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, the 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 the SWE within the mountainous CRB. In the initial approach, the SWAT was calibrated solely using observed streamflow data. In the second (multi-variable) calibration approach to the 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 the 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.

The selected parameters were divided into three groups comprising snow parameters, basin parameters and elevation lapse rate parameters. In earlier research endeavours [10,11,34,48] aimed at stimulating the SWE, SNOCOVMX and SNO50COV were deemed insensitive and consequently omitted. Nevertheless, in our investigation, we opted to incorporate these parameters (SNOCOVMX and SNO50COV) to assess their influence on streamflow simulation. Interestingly, it was observed that the inclusion of SNOCOVMX and SNO50COV, along with five additional snow parameters (SFTMP, SMTMP, SMFMX, SMFMN, TIMP) had no discernible impact on the simulation of streamflow. We excluded SNOCOVMX and SNO50COV, retaining only the five snow parameters for further analysis. The calibration process then proceeded with the remaining 16 parameters, as outlined in

Table 2. The SWAT model parameters including 5 snow parameters, 7 basin parameters and 2 elevation band parameters used in this study to calibrate model in the CRB, their initial values and final calibrated values.

ID	Name	Description	Range
Snow Parameters	V_SFTMP	Snowfall temperature (_C)	0–5
	V_SMTMP	Snowmelt base temperature (_C)	−5–5
	V_SMFMX	Melt factor for snow on 21 June (mm H2O/_C-day)	−5–5
	V_SMFMN	Melt factor for snow on 21 December (mm H2O/_C-day)	0–5
	V_TIMP	Snowpack temperature lag factor	0–1
Basin Parameters	R_CN2	SCS curve number	−0.2–0.2
	V_GW_DELAY	Groundwater Delay (days)	0–500
	V_ALPHA_BF	Baseflow alpha factor (days)	0–1
	R_SOL_AWC	Available water capacity of the soil layer (mm H2O/mm soil)	−0.5–0.5
	V_GWQMN	shallow aquifer required for return flow to occur (mm)	0–5000
	V_GW_REVAP	Groundwater “revap” coefficient	0.02–0.2
	V_RCHRG_DP	Deep aquifer percolation fraction	0–1
	R_ESCO	Soil evaporation compensation factor	0–1
	V_CH_K2	Effective hydraulic conductivity in the main channel (mm/h)	0–500
Elevation Lapse rate Parameters	PLAPS	Precipitation lapse rate (mm H2O/km)	0–1000
	TLAPS	Temperature lapse rate (mm H2O/km)	−8–−4

. 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 basins, which affect both temperature and precipitation, the SWAT model 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 the 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 basin before using it as forcing data to derive the SWAT model (

Table 3. Validation of satellite gridded and reanalysis precipitation (SGRP) data with station precipitation at Muzaffarabad (station-01) and Garhi Dupatta (station-02) in the CRB.

Satellite-Reanalysis Precipitation	Observed Precipitation	RMSE	PBIAS (%)	NSE	KGE	R2
CMORPH	Station-01	617.56	419.7	−34.20	−4.10	0.27
	Station-02	634.51	449.9	−42.80	−4.81	0.13
TRMM	Station-01	128.02	−79.7	−0.53	−0.12	0.51
	Station-02	118.24	−78.5	−0.52	−0.16	0.48
CPCUPP	Station-01	108.87	−52.0	−0.09	0.11	0.26
	Station-02	105.26	−54.0	−0.20	0.08	0.15
APHRODITE	Station-01	91.10	−32.0	0.23	0.32	0.36
	Station-02	81.96	−28.7	0.27	0.38	0.37
ERA5-Land	Station-01	74.26	−0.7	0.49	0.56	0.49
	Station-02	62.12	5.1	0.58	0.64	0.58

).

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 showed very high negative and positive bias, respectively against both observed precipitation stations as shown in Figure 5. CMORPH was highly overestimated (Figure 5) while all remaining precipitation products underestimated precipitation when compared with observed station data (Figure 5). 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

Figure 6a depict the comparison of global and regional (e.g., UCLA and KRA) SWEproducts and Figure 6b illustrate the snow cover area (SCA) at monthly scale in the CRB. There is a notable difference in SWE amounts over basin scales 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 underestimate SWE when compared to the higher resolution satellite-derived UCLA-SWE product. GLDAS, a reanalysis of the SWE product exhibits some similarity with UCLA-SWE between June and December while, it underestimates the 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 in terms of quantity and pattern during the snow accumulation and snowmelt periods.

TWhen examining the 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 7 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 the 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 the SWE values of 199 mm and 239 mm, respectively, while UCLA reports underestimation of the 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 mm, contrasting with UCLA-SWE values of 182 mm, 120 mm, 188 mm, and 220 mm, respectively. However, it’s worth noting that in sub-basins 3, 4, and 11, both KRA-SWE (181 mm, 118 mm, and 205 mm, respectively) and UCLA-SWE (182 mm, 122 mm, and 208 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 8. It is evident that CMORPH and TRMM precipitation inputs did not yield adequate discharge records in the CRB (Figure 8a,b, 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 8c). However, the 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. The SWAT model simulations derived with other than ERA5-Land precipitation products achieved unsatisfactory results (NSE ≤ 0.50).

The SWAT model ran with the APHRODITE, CPCUPP and ERA5-Land precipitation data was further evaluated for its performance in simulation of the SWE in the CRB as depicted in Figure 9. These simulations were then compared with two existing SWE products (KRA-SWE and UCLA-SWE). The SWE output from the 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 that used 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 the streamflow and SWE in the CRB.

3.4. Single Variable Calibration Approach

After thoroughly investigating the impact of precipitation data on the SWAT model simulations (e.g., streamflow and the 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, the SWAT snow parameters were calibrated manually.

Table 4. Statistical evaluation of the simulated SWE and streamflow using the SWAT model with UCLA and KRA SWE products during calibration and validation period.

Metrics	Calibration	Validation	Calibration	Validation
	(2000–2011)	(2012–2015)	(2000–2011)	(2012–2015)
	UCLA Streamflow		KRA Streamflow
PBIAS	−13.9	−8.1	−8.2	−5.2
NSE	0.33	0.42	0.71	0.77
R2	0.36	0.43	0.66	0.78
KGE	0.46	0.5	0.77	0.87
	UCLA SWE		KRA SWE
PBIAS	3.9	5.4	11.5	−10.7
NSE	0.64	0.87	0.94	0.94
R2	0.65	0.88	0.96	0.97
KGE	0.78	0.91	0.86	0.87

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 (Figure 10). There is also promising evidence in the validation period with respect to base flows.

The simulated SWE from the SWAT model (SWAT-SWE) was compared with snow reference data (KRA-SWE and UCLA-SWE) in the CRB as shown in Figure 11. SWAT-SWE consistently underestimated peaks in the 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 12 illustrates the streamflow simulations conducted with reference SWE data (KRA-SWE and UCLA-SWE). The simulated streamflow results using KRA-SWE product to calibrate the 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 12a) weakening streamflow simulations. However, the multi-calibration techniques improved the simulated streamflow performance using KRA-SWE as indicated by the statistical indices provided in Table 4.

Figure 13 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 the 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 matched 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 the 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 for 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.

3.6. Snow Water Equivalent Simulations (Sub-Basins)

The comparison of the SWAT-simulated SWE (SWAT-SWE) with KRA-SWE and UCLA-SWE is presented in

Table 5. Evaluation of simulated SWAT-SWE in each sub-basin of the Chenab River Basin (CRB) using UCLA-SWE and KRA-SWE data. NSE ≥ 0.5 are shown bold.

	PBIAS%	NSE	R2	KGE	PBIAS%	NSE	R2	KGE
Sub-Basins	KRA-SWE				UCLA-SWE
BSN01	10.4	0.82	0.95	0.66	17.5	0.76	0.81	0.79
BSN02	5.4	0.71	0.96	0.52	24.7	0.6	0.64	0.66
BSN03	72.2	−2.29	0.89	−0.59	−31	0.62	0.76	0.46
BSN04	60.9	−2.05	0.92	−0.63	−31.9	0.55	0.68	0.42
BSN05	−9.2	0.6	0.92	0.48	21.8	0.39	0.44	0.55
BSN06	−77.7	−0.24	0.21	−0.01	258.8	−8.9	0.7	−2.46
BSN07	330.6	−17.82	0.83	−3.69	−34.5	0.52	0.62	0.39
BSN08	143	−3.14	0.88	−1.11	−14.7	0.51	0.54	0.5
BSN09	−10.9	0.22	0.63	0.51	138.1	−2.11	0.61	−0.85
BSN10	57.3	−0.3	0.91	−0.01	9239	−761.45	0	−93.22
BSN11	348.9	−30.06	0.1	−3.5	−43	0.51	0.82	0.37
BSN12	0.4	0.72	0.93	0.57	−66.8	−0.05	0.51	0.04
BSN13	−24.8	0.52	0.82	0.59	9.4	0.32	0.33	0.46
BSN14	NA	NA	NA	NA	NA	NA	NA	NA
BSN15	NA	NA	NA	NA	NA	NA	NA	NA
BSN16	NA	NA	NA	NA	NA	NA	NA	NA
BSN17	−28.7	0.65	0.84	0.68	46.5	0.41	0.57	0.47
BSN18	31.4	−0.21	0.63	0.33	−45.1	0.43	0.75	0.34
BSN19	23.4	−0.01	0.61	0.43	−40.7	0.43	0.64	0.4

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 SWESWE (NSE > 0.5). Based on R² (R² > 0.60), the 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 the 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 14 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 the 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 [50] 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.

In snow-dominated basins, conventional calibration of snow parameters often involves simultaneous adjustment with hydrological parameters against observed discharge, disregarding snow conditions [51]. However, some studies highlight the benefit of integrating snow observations in a multi-objective calibration to enhance SWE simulation [11,52,53]. In previous studies, satellite-derived snow cover area from MODIS and Landsat, is commonly used for calibration, despite its potential limitations such as cloud cover and forest interference [52]. Although in situ snow data are less frequently employed due to challenges in heterogeneous snow conditions, they have proven valuable in calibrating hydrological models, particularly in glacierized basins [54]. This study addresses the absence of in situ snow data in the study area by using modelled and remote sensing SWE data to calibrate and validate the SWAT model, building on the success of previous studies that improved the SWE simulations with the SWAT model using in situ measurements [45,53]. This study revealed that satellite-derived UCLA-SWE performed less well than the modelled snow KRA-SWE product in estimating the 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 does not 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 does not contribute to the basin’s overall flow, leading to discrepancies in flow simulations using the SWAT model when calibration was performed with UCLA-SWE data.

We applied a methodology akin to that used by [10,11,34,49] for calibrating snow parameters to simulate SWE and streamflow in the study basin. However, in the absence of in situ data in the study basin, we utilized two newly developed datasets as reference data for snow parameters calibration. Our approach surpasses previous methods, as we manually adjusted the unrealistic combination of SMFMX and SMFMN in R-SWAT resulting from autocalibration in SWAT-CUP, as highlighted by [10]. 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 modelling snowpack. The implementation of a multivariate calibration method, which incorporated snow reference data from KRA-SWE product, resulted in enhanced predictions for both SWE and streamflow. However, it’s important to note that UCLA-SWE product performed well in simulating the SWE in certain individual sub-basins (Figure 14) 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 the SWE when utilizing a multivariate calibration approach.

The calibration methods detailed in this study have improved streamflow and the 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 the 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 [55]. 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 the 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 the SWE simulation. The calibrated parameters used for streamflow simulation may help compensate for the SWE discrepancies and potential deteriorations, such as earlier peak flows, if the SWE is enhanced by the inclusion of reference snow data [56].

Various approaches to measuring snowpack come with their own advantages and limitations. In the case of in situ monitoring for SWEs 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 modelling, especially when integrated with ground observations.

This study brought attention to the challenges associated with using the 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 the SWE (even with uncalibrated SWAT model) which, in turn, supports the model’s utility for simulating the 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 basin during the warm season. This study shows that SWAT can accurately simulate SWE and streamflow in snow dominated basins using snow reference data to calibrate the snow parameters especially in high-altitude systems. The capability of the SWAT model is traditionally 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 product as the reference snow data, both the simulated streamflow and SWE demonstrated good agreement with observed streamflow and reference KRA-SWE product. Significant enhancements in the accuracy of simulated SWEs 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 basins, 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 the SWE.

The primary aim in refining our snow modelling approach was to enhance the accuracy of the basin-wide SWE characterization by integrating different SWE products individually. By adjusting snow parameters, we were able to achieve a good result in estimating the 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.

Our multivariate calibration/validation approach permitted a comprehensive evaluation of snow modelling across the basin 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 the 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. Multivariate calibration, demanding ample meteorological and hydrological data, faces challenges in regions with limited data availability. Uncertainties in input data and parameter estimates, as well as correlations among parameters, may not be adequately accounted for, impacting model reliability. Additionally, errors in observed data can propagate into the calibrated model, affecting predictive accuracy. Spatial variability of the SWE, especially in mountainous regions, may not be effectively captured.

This study holds relevance for ongoing scientific research focused on estimating SWEs 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 basin within the Himalayan region. Overall, the integration of remote sensing data into computational model enhances the accuracy of SWE simulations, particularly at the basin scale. This approach provides 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 basins 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 systems across the world.