Evaluating the dynamics of climate and anticipated snow cover changes on the hydrological response of rivers in the Hindu Kush Himalaya region

The Hindu Kush Himalayas (HKH) are renowned as Asia's water towers, serving as the origins for 10 major rivers. Snowmelt and lateral flow through groundwater are primary contributors to the lean flow of these rivers Mukherji et al. (2015). These riverine ecosystems are vital, sustaining the food and energy needs of 1.3 billion downstream people. However, ongoing climatic changes, socio-economic development, and urban migration are exerting increased pressure on natural resources, particularly water. A substantial portion of the annual stream flow in the Himalayan region is attributed to snow cover and glaciers Shrestha (2011). Snowfall and snowmelt are dominant hydrological processes in snow-fed basins Nazeer et al. (2022). During winter, extensive areas are blanketed by snow, which melts as temperatures rise toward the end of winter. In lean periods, the melting of snow and glaciers significantly contributes to water quantities. However, climate change is accelerating the melting of snow and glaciers and reducing snow accumulation due to shifts in precipitation patterns, resulting in increased variability in river discharge (Azizi & Asaoka 2020). A study by the Inventory of Carbon and Energy indicates an average decrease of 16% in the flow of water from snowmelt to rivers in high-mountain Asia between 1979–1999 and 1999–2019. Global mean surface temperature, encompassing both terrestrial and marine environments, has exhibited a consistent upward trend since 1861. Global warming is anticipated to elevate temperatures by 1.7–6.3 °C in the HKH region by the end of the 21st century (Lutz et al. 2016). The rise in minimum temperatures amplifies the snowmelt process and alters the snow cover area (SCA) by shifting the snow line to higher elevations. Higher elevations are expected to experience a greater increase in temperature. According to the IPCC AR5 assessment, there is a clear and significant decrease in both the range and extent of snow cover in the northern hemisphere. Satellite data indicate a high probability of approximately 10% reduction in snow-covered areas since the late 1960s (Brown & Robinson 2011). These consequences, coupled with future climate change impacts, are expected to exacerbate ecosystem degradation, disturbances in water quality, and the loss of snow cover, glaciers, and sea ice (Khadka et al. 2014). The snow cover and glaciers in the HKH are recognized as some of the fastest-thinning globally (Desinayak et al. 2021). Melt waters from snow and ice in the HKH region play a crucial role in providing water for downstream farming, hydropower generation, and domestic use (Scott et al. 2019). Consequently, the accelerated melting of snow cover and glaciers increases the incidence of natural hazards (Azizi & Asaoka 2020).

Numerous studies have been conducted to understand the natural hydrological response of river basins by establishing the relationship between climate and watershed characteristics in the hydrological process, utilizing rainfall runoff models (Talchabhadel & Karki 2019; Aryal et al. 2023; Ramahaimandimby et al. 2023; Devkota et al. 2024). Most research emphasizes that the application of hydrological models relies on available data and required outputs. For more accurate estimations, mathematical models are considered more realistic (Ghonchepour et al. 2021). In high mountainous areas with limited data availability, semi-distributed models or lump models are deemed more reliable (Bajracharya et al. 2018; Mishra et al. 2018; Pandey et al. 2019; Sinha et al. 2019). Researchers globally, as well as in the Nepalese Himalaya region, have successfully utilized the soil and water assessment tool (SWAT) for hydrologic modeling and water resources management within watersheds characterized by diverse climatic and topographic conditions. The SWAT model has demonstrated its efficacy in such applications. Upon calibrating the model parameters using a physically observed database, the SWAT model has shown significantly improved performance (Van Liew et al. 2007; Liu et al. 2014; Aryal et al. 2019; Marahatta et al. 2021; Talchabhadel et al. 2021; Malik et al. 2022; Devkota et al. 2024). These investigations highlight the robustness and reliability of the SWAT model in various geographical settings. In a study by Thapa et al. (2017) focusing on the Bagmati Basin, which is adjacent to the Marsyangdi River Basin (MRB), the SWAT model's outputs exhibited greater accuracy compared to alternative hydrological models such as the Hydrologiska Byrans Vattenbalansavdelning (Bergström 2006) and the block-wise use of the TOPMODEL model with Muskingum-Cunge routing (BTOPMC) (Takeuchi et al. 1999, 2007). This underscores the SWAT model's superiority in capturing the complexities of water balance components in diverse hydrological settings.

The MRB, a crucial water source in the HKH region, is currently undergoing a reduction in its SCA at a rate of 0.34% annually, as revealed by a snow cover map spanning the research region from 2000 to 2019. Given the expected alterations in snow dynamics in the future, it is imperative to develop a comprehensive understanding of the existing snow cover dynamics and their potential implications for basin hydrology. Despite the relative abundance of studies on glaciers in the Himalayas, research specifically focused on snow cover is scarce (Shrestha 2011). A conspicuous knowledge gap persists regarding the impact of changes in snow cover on basin hydrology, even though numerous studies have investigated the effects of land use and land cover (LULC) changes on watershed hydrology.

This paper endeavors to bridge this knowledge gap by projecting future snow cover in the MRB through logistic regression and establishing its impact on future river flow using a hydrological model. It is essential to establish a balanced understanding of the demand and supply dynamics at the whole river basin or sub-basin scale, considering factors such as socio-economic changes, urbanization, demographic shifts, climate variations, and alterations in agricultural practices. This is crucial for fostering sustainable water resources development. To attain this understanding, the first step is to identify and quantify the driving forces (climate and SCA change) that contribute to imbalances in resources (Mukherji et al. 2015). In this context, our study addresses these challenges by examining two key driving forces, namely climate and snow cover variation, utilizing the semi-distributed, physically based SWAT model (Arnold et al. 1998; Neitsch et al. 2011) over a snow-covered watershed. The study aims to analyze the potential effects of future changes on the hydrology of the basin in both contexts. The specific goals of this research are threefold: (i) to assess the impact of climate change on the hydrology of the basin in the near future (NF), mid-future (MF), and far future (FF); (ii) to quantify the variation in SCA in the future and its implications for river hydrology; and (iii) to evaluate the integrated impact of future SCA variation and climate change on the river hydrology of MRB. This research addresses the critical need to comprehend how shifts in snow cover could influence basin hydrology, thereby having far-reaching implications for water resource planning, development, and management strategies in the region. Establishing a robust knowledge base on the intricate relationship between snow cover fluctuations and hydrology is imperative for effectively adapting to the challenges presented by climate change in the Himalayas.

Topographic, soil, land use/cover maps, hydrological, and meteorological data were collected from various sources as outlined in Table 1. These datasets were then preprocessed for analysis. The digital elevation model (DEM) was used to delineate the watershed and generate the river stream network. Higher-resolution LULC maps included areas of mixed forest, water bodies, snow and glacier-covered regions, barren land, built-up areas, and cropland. Notably, snow cover accounted for 20% of the total area.

The simulation of hydrological and environmental processes within the river basin was conducted using the continuous-time, semi-distributed, process-based SWAT model, as detailed by Arnold et al. (1998). Unlike regression models, SWAT requires comprehensive information about land management practices, terrain, weather, soil properties, and vegetation within the watershed. The study area was divided into hydrologic response units (HRUs), elevation bands, and sub-basins to capture the spatial heterogeneity of the basin. Specifically, 56 sub-basins and 2,593 HRUs were generated using spatially distributed data for LULC, soil, and slope, with zero threshold values applied in a GIS environment. The slope classes were categorized into five ranges: 0–5%, 5–10%, 10–20%, 20–50%, and greater than 50%. To replicate snowmelt and orographic temperature and precipitation distribution, three elevation bands were created, each with a mean elevation of 500 m. Meteorological input for model development included time-series data as specified in Section 2.2. The auto-calibration process employed the SUFi-2 algorithm and SWAT-CUP to refine and reduce the effective range of each parameter. Global sensitivity analysis was used to assess parameter sensitivity, with the most sensitive parameters identified by the highest t-stat value and the lowest p-value (less than 0.05). Calibration involved parameterizing the model to specific conditions using the first 37 parameters. Given the sensitivity of snow parameters in the mountainous basin, seven snow parameters (SMFMX, SMFMN, TIMP, SNOCOVMX, SNO50COV, SFTMP, and SMTMP) were included. The model was calibrated using daily discharge data at the basin's outlet for the years 1988–2005 and validated for the years 2006–2015.

Performance metrics, including the coefficient of determination (R²), Nash–Sutcliffe simulation efficiency (NSE), percentage bias (PBIAS), R-factor, and P-factor, were employed to evaluate the model's performance (Table 2). The calibration process continued until acceptable performance metric values were achieved, after which the same parameters were used for validation. SWAT-CUP-derived parameter values were manually assigned in the SWAT model after calibration and validation. Additionally, the most sensitive parameters identified during auto-calibration were adjusted first during manual calibration for optimal statistical parameter values (NSE, PBIAS, and R²). The model's performance evaluation included daily and monthly simulations, assessing hydrographs, scatter plots, flow duration curves, statistical parameters, and water balance (actual vs. simulated) at daily, monthly, and annual scales. Various visual assessments were conducted, such as examining hydrograph peaks, time to peak, hydrograph shape, and base flow, to comprehensively evaluate the model's performance (Pandey et al. 2020; Devkota et al. 2024).

The analysis incorporates the latest scenarios known as shared socio-economic pathways (SSPs), as highlighted in the IPCC's Sixth Assessment Report (AR6) and the Coupled Model Intercomparison Project Phase 6 (CMIP6). Each dataset output from the GCMs represents distinct potential future scenarios. For this study, we selected the high-end emission scenario (SSP585), which reflects a fossil fuel-driven development ideology, and the middle-of-the-road scenario (SSP2), combined with the moderate emission scenario (SSP245). GCM models provide reliable insights into past, present, and future climate patterns. The baseline period from 1980 to 2009 serves as a reference for analyzing future periods, which are further divided into NF (2021–2045), MF (2046–2070), and FF (2071–2095).

Eleven GCMs were selected based on a literature review to anticipate future climate conditions in the HKH region and similar basins in Nepal. These GCM runs were evaluated using performance measures (PBIAS, RSR, and NSE) against observed data. The World Research Climate Programme provided GCM outputs in NetCDF format, which were downscaled for compatibility with the SWAT model using R-studio. A prior performance method assessed the appropriateness of each GCM, resulting in an aggregate rating and ranking based on performance scores, as shown in Table 3. Recognizing the diversity of potential scenarios depicted by different GCMs, the study employed an ensemble of the top three selected GCMs – ACCESS-ESM1-5, EC-Earth3, and CanESM5. This ensemble approach aimed to minimize uncertainties in anticipating future climatic conditions, aligning with practices in hydro-climatic studies that often utilize multiple GCMs to provide a more comprehensive understanding. Quantile mapping bias correction techniques were employed to rectify systematic distributional biases in climatic outputs derived from climate models. Among these techniques, the robust empirical quantile method proved effective for bias correction in climate data. This approach enabled the derivation of future time-series of climate variables from the raw GCMs under the SSP245 and SSP585 trajectories.

The three selected GCMs for precipitation, maximum temperature, and minimum temperature were evaluated based on their average values. Statistical performance indicators (Table 3) reveal that ACCESS-ESM1-5 scored 2.57, EC-Earth3 scored 2.56, and CanESM5 scored 2.53. These ratings reflect the performance and suitability of each GCM in accurately capturing and simulating the specified climatic parameters. The outcomes provide a concise overview of the prioritized GCMs, guiding further considerations in the study.

The Dyna-CLUE model, an advanced iteration of the CLUE-S model, was specifically designed to simulate land use changes, particularly regarding snow cover variation, in a spatially explicit manner based on an empirical assessment of location suitability. This dynamic model facilitates the simulation of interactions and competitions between the temporal and spatial dynamics of various land use types, featuring two distinct modules: one for demand and the other for spatial allocation. The historical changes in land use cover, especially SCA, and the evaluation of projected demand were assessed using LULC planning scenarios or historical patterns of LULC types (Verburg & Overmars 2009). Concurrently, the second module transformed the spatial LULC requirements of the study region.

In this study, the Dyna-CLUE model was utilized to project the snow cover map at a regional scale. The iterative determination of optimal results in the model relied on four inputs: geographic features, land-use requirements (demand), conversion settings specific to each land-use type, and spatial policies and constraints. The LULC map of 2010 for the watershed served as the base map to project the yearly snow cover map up to 2100. Logistic regression models (SPSS Model) were constructed to estimate location suitability requirements based on the spatial interaction of each LULC type with a set of driving forces for LULC change. These driving factors included elevation, slope, aspect, rainfall, T_max, T_min, population density, distance to river, distance to road, and soil map. To facilitate the Dyna-CLUE model, all layers in the GIS environment were required to overlap properly and be converted into raster to ASCII format. In the trend scenario, no spatial policies were implemented, and the restricted region (forest to other) encompassed the Annapurna Conservation Area, where changes were not permitted. The rules of conversion between LULC types were determined by conversion parameters, including conversion flexibility and the conversion matrix. A specific LULC type could convert to other types based on its conversion flexibility, with the conversion elasticity ranging from 0 (greatest conversion possibility) to 1 (no conversion possibility) (Trisurat et al. 2010). Through a trial-and-error approach, values of conversion parameters (0.6, 1, 0.8, 1, 0.9, 0.7, 0.8, and 0.4) were allocated to forest, water body, snow and glacier, riverbed, built-up area, cropland, barren land, and grassland, respectively (Lü et al. 2020). Water bodies and riverbeds were given high elasticity values due to their inherent resistance to transformation into another class. The conversion matrix employed in the MRB was established after a comprehensive examination of past patterns, specifically the changes that occurred between 2000 and 2010. A value of ‘1’ indicates a viable conversion, while a value of ‘0’ indicates an impossible conversion in the conversion matrix. The Dyna-CLUE model has not yet been implemented in the MRB. Therefore, conversion parameters were prepared based on an analysis of previous research on studies with comparable features and expert judgment. The following references were consulted: (Verburg & Overmars 2009; Trisurat et al. 2010; Tizora et al. 2018; Lü et al. 2020; Khoi et al. 2021). Future demands were estimated using a simple extrapolation method based on annual historical fluctuations in the area of each LULC type, particularly snow cover variation. The receiver operating characteristic (ROC) was employed as a tool to evaluate the goodness-of-fit of logistic regression models, with the area under the curve ranging from 0.5 (complete randomness) to 1.0 (perfect match), serving as a summary indicator to assess the model's overall effectiveness. The projected LULC map (2019) and the LULC map (2019) from ICIMOD were compared to validate the model using Kappa statistics for a reliability assessment of land-use change simulations, where a Kappa statistic of 0 denotes chance agreement, and a statistic of 1 denotes perfect agreement.

In this study, we aimed to examine the impact of climate change on stream flow and various water balance components within a watershed. We employed a comprehensive approach, utilizing a calibrated SWAT model to simulate anticipated future conditions based on quantified climate data. We maintained the baseline snow cover map, derived from LULC data, at constant levels throughout the study. The well-calibrated SWAT model was used to estimate stream flow and water balance for three distinct time periods: NF, MF, and FF. These estimates were rigorously compared against the baseline stream flow and water balance components. To specifically assess the impact of changing snow cover, we used the 2010 SCA map as a reference and incorporated future snow cover maps for 2030, 2050, and 2090, calibrated by the Dyna-CLUE model. By maintaining climatic data consistent with the baseline period (1980–2017), we systematically compared SWAT model outcomes using future snow cover maps against the 2010 baseline information. This comparative analysis enabled us to discern the potential effects of altered snow cover on hydrological dynamics and the overall water balance in the study area. The SWAT model was supplied with projected snow cover maps for 2030, 2050, and 2090, along with corresponding climate data for the NF (2021–2045), MF (2046–2070), and FF (2071–2095) periods. Following the model run, we conducted a meticulous comparison, contrasting the projected future stream flow and water balance components with the baseline hydrology of the study region. This comprehensive evaluation aimed to gauge the integrated influence of future changes in SCA and climate change on the hydrology of the basin, providing valuable insights into potential shifts in water resources under various future scenarios.

The model shows a more accurate simulation of low flows compared to high flows, a limitation acknowledged in previous studies (Mishra et al. 2018). Difficulty in replicating peak flows may be attributed to the complexity introduced by days with multiple storms, posing challenges for the existing curve number technique to predict runoff effectively. As depicted in Figure 3, a significant correlation between simulated and observed daily and monthly variables is evident during both calibration and validation stages, with NSE, PBIAS, and R² values within acceptable ranges for daily and monthly simulations, respectively (Table 4). Combined graphical and statistical analyses affirm the model's capability to achieve acceptable runoff simulation during calibration and validation stages. These results suggest that the SWAT model can be effectively utilized for hydrological modeling in the MRB, serving as a valuable tool for assessing water resources and studying the basin's water balance.

During the baseline period, over 75% of the average annual rainfall occurs during the monsoon season (June–September), with 16% falling in the pre-monsoon season (March–May), 4% in the post-monsoon season (October–November), and 5% in the winter season (December–February). This highlights the critical importance of monsoon rainfall in the MRB and HKH regions. The average annual precipitation for each future period shows an upward trend relative to the baseline for both scenarios. Table 5 illustrates the percentage increase in seasonal precipitation for NF, MF, and FF compared to the baseline value under both SSP245 and SSP585 scenarios, revealing a rising trend in daily precipitation from the baseline period (1988–2019) to NF, MF, and FF. Regarding temperature, the ensemble (arithmetic mean) of T_max and T_min from the top three selected GCM models was calculated and used for further analysis. As depicted in Table 5, the average T_max values increase by 1.49 °C in NF, 2.51 °C in MF, and 3.24 °C in FF for SSP245, and by 1.65 °C in NF, 3.50 °C in MF, and 6.05 °C in FF for SSP585, relative to the baseline data. Meanwhile, T_min shows an increase of 1.90 °C in NF, 2.73 °C in MF, and 3.39 °C in FF for SSP245, and 2.04 °C in NF, 3.59 °C in MF, and 5.61 °C in FF for SSP585, compared to the baseline values.

These projected decadal and seasonal temperature variations, particularly in maximum and minimum values, intensify the snowmelt process and shift the snow coverage line in the higher Himalaya region. The increased snowmelt rate enhances pre-monsoon flow, while the monsoon season experiences increased runoff volume from the catchment due to reduced snow coverage. These changes have multifaceted impacts, altering recharge rates and lateral flow at the catchment level, and highlight the complex interplay of phenomena resulting from the evolving climate context (Thapa et al. 2021).

However, in future climate scenarios, the increase in precipitation patterns leads to a corresponding rise in stream flow values during this season. Despite the initially low base flow, the numeric values show a high percentage change, though the impact is not significantly pronounced. The increased water availability during the lean period, however, has a more substantial impact on downstream water resource projects. In terms of the water balance component, the average annual water yields in NF, MF, and FF, which estimate the amount of freshwater entering streams and rivers from precipitation (including rain, snow, and snowmelt), rise by 7, 34.74, and 27.81% under SSP245, and by 19.95, 50.33, and 74.18% under SSP585. This increase is influenced by future precipitation increases of 24.80% in SSP245 and 62.60% in SSP585. The two climate scenarios predict increased precipitation and reduced evapotranspiration, leading to higher surface runoff, lateral flow, and groundwater contribution in MF. However, due to the projected temperature rise, evapotranspiration is expected to increase by 16.07% under SSP245 and by 29.02% under SSP585 in FF. Under SSP245, water yield increases the most in the pre-monsoon season for each subsequent period. This is likely due to the highest expected rise in temperature and precipitation during the pre-monsoon season compared to the baseline under both scenarios for each future period. Post-monsoon and winter surface runoff is lower due to increased snow accumulation, which acts as a temporary reservoir and melts during the spring thaw. Despite a gradual decrease in snowmelt contribution, our analysis suggests a limited discernible impact on overall stream flow. The watershed, being predominantly rainfall-driven, indicates that changes in snowmelt are not significantly influencing the broader stream flow dynamics, though this aspect is not overlooked. This nuanced understanding contributes to a more comprehensive grasp of the hydrological response to climate change scenarios in the MRB.

The impact of future changes in SCA on water balance components reveals a complex pattern. In the NF, there will be a slight increase in annual water yield, followed by a decrease in the long term. Conversely, evapotranspiration, which involves the loss of water through evaporation and plant uptake, is projected to decrease initially and then increase in the distant future. Annual surface runoff, representing water flowing over the land, is expected to rise by 11.4% in 2030 and 8.07% in 2050, but then decrease by 2.18% in 2090 (Table 7). This trend can be attributed to a reduction in snow cover, replaced by more bare ground and woodland.

Likewise, average annual water yield exceeded baseline levels for each future period, indicating a continuous rise in precipitation over time under both SSP245 and SSP585 scenarios, with rates of 8.41% (NF), 34.92% (MF), 28.14% (FF), and 20.5% (NF), 50.5% (MF), and 74.9% (FF), respectively. Surface runoff and groundwater flow were anticipated to increase more in the mid and NF compared to the far future, accompanied by a 15.48% increase in evapotranspiration under the SSP245 scenario. As temperatures rise and snow cover diminishes, increased evapotranspiration from forest and barren lands impacts both overland flow and groundwater recharge. Water yield exhibited an increasing trend in each season for every future period, albeit at varying rates, likely linked to precipitation variations. The combination of higher temperatures and increased precipitation may accelerate snowmelt, influencing runoff timing. Reductions in SCA and decreased snowfall rates contributed to an 11.82 and 38.36% decrease in snow water equivalent (SWE) in the far future under SSP245 and SSP585 scenarios, respectively. The decline in SWE value was particularly pronounced in post-monsoon and winter seasons compared to pre-monsoon and monsoon seasons, reflecting higher precipitation and temperature increases during the latter seasons. This shift in snowmelt dynamics underscores the changing role of snow in future river flow scenarios, driven by decreasing SCA and climate change-induced alterations in precipitation and temperature patterns.

Climate change continues to increasingly impact the HKH Region, raising concerns about its potential repercussions on river systems. These include altered flow regimes, shifts in the timing and intensity of river discharge, and implications for downstream water resource management. Understanding the dynamics of climate and projected changes in snow cover on river hydrology is crucial for informed water resource management and adaptation strategies. This paper explores historical climate patterns and anticipated changes in snow cover within the HKH Region, focusing on our study area. We developed a methodological framework using the SWAT model, coupling historical climate data, projected snow cover changes, and their integrated impact on basin water balance. The robust performance of the SWAT model, validated through calibration and validation metrics, confirms its effectiveness in simulating diverse runoff components. Our projections indicate a concerning decline in SCA by 2.03, 4.03, and 14.20% in the years 2030, 2050, and 2090, respectively.

Examining climate change alone reveals significant increases in average annual stream flow, projecting increments of 84 and 130% in FF under SSP245 and SSP585 scenarios, respectively, primarily due to expected rises in precipitation. This increase is most pronounced during the pre-monsoon season, with projected increases of 150 and 183% in FF under SSP245 and SSP585 scenarios, respectively, highlighting the season's sensitivity to future precipitation changes. Isolating the impact of projected SCA alone shows a modest increase in annual average stream flow, indicating rises of 4.1, 4, and 3.5% relative to baseline in the NF, MF, and FF periods, respectively. However, the anticipated decrease in stream flow during winter by 20.77% raises concerns for water availability during the lean period. Conversely, reduced snow cover prompts an increase in groundwater contribution and base flow, suggesting potential benefits for future groundwater-related projects. Analyzing the combined impact on annual average stream flow under SSP245 and SSP585 scenarios reveals increases of 59% (NF), 85% (MF), 85% (FF), and 70% (NF), 104% (MF), and 131% (FF), respectively, attributed to rising precipitation rates. Similarly, average annual water yield surpasses baseline levels for each future period, indicating a continuous increase in precipitation under both scenarios. The observed rise in surface runoff, groundwater, and evapotranspiration in NF and MF, accompanied by a drop in FF (SCA2090), underscores the complex interplay of rising temperatures, diminishing snow cover, and changing land use.

In conclusion, this study addresses a significant knowledge gap regarding the specific implications of future snow cover variations on MRB hydrology. By unraveling the intricate relationships between climate change, snow cover variability, and hydrological responses, our research offers valuable insights for water resource management and emphasizes the necessity for adaptive strategies in navigating evolving hydrological conditions. As we confront an uncertain climatic future, the findings from this study provide a foundation for evidence-based decision-making in sustainable water resource management, particularly in mountainous regions such as the HKH Region.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.