Assessing Climate Change Impacts on Streamflow and Baseflow in the Karnali River Basin, Nepal: A CMIP6 Multi-Model Ensemble Approach Using SWAT and Web-Based Hydrograph Analysis Tool

Abstract

Our study aims to understand how the hydrological cycle is affected by climate change in river basins. This study focused on the Karnali River Basin (KRB) to examine the impact of extreme weather events like floods and heat waves on water security and sustainable environmental management. Our research incorporates precipitation and temperature projections from ten Global Circulation Models (GCMs) under the Coupled Model Intercomparison Project Phase 6 (CMIP6). We applied thirteen statistical bias correction methods for precipitation and nine for temperatures to make future precipitation and temperature trend projections. The research study also utilized the Soil and Water Assessment Tool (SWAT) model at multi-sites to estimate future streamflow under the Shared Socioeconomic Pathway (SSP) scenarios of SSP245 and SSP585. Additionally, the Web-based Hydrograph Analysis Tool (WHAT) was used to distinguish between baseflow and streamflow. Our findings, based on the Multi-Model Ensemble (MME), indicate that precipitation will increase by 7.79–16.25% under SSP245 (9.43–27.47% under SSP585) and maximum temperatures will rise at rates of 0.018, 0.048, and 0.064 °C/yr under SSP245 (0.022, 0.066, and 0.119 °C/yr under SSP585). We also anticipate that minimum temperatures will increase at rates of 0.049, 0.08, and 0.97 °C/yr under SSP245 (0.057, 0.115, and 0.187 °C/yr under SSP585) for near, mid, and far future periods, respectively. Our research predicts an increase in river discharge in the KRB by 27.12% to 54.88% under SSP245 and 45.4% to 93.3% under SSP585 in different future periods. Our finding also showed that the expected minimum monthly baseflow in future periods will occur earlier than in the historical period. Our study emphasizes the need for sustainable and adaptive management strategies to address the effects of climate change on water security in the KRB. By providing detailed insights into future hydrological conditions, this research serves as a critical resource for policymakers and stakeholders, facilitating informed decision-making for the sustainable management of water resources in the face of climate change.

1. Introduction

Over the past century, the global average temperature has risen by 0.74 °C with increased extreme precipitation events and erratic weather patterns, as reported by the Intergovernmental Panel on Climate Change (IPCC) [1]. This warming trend is even more pronounced in the Himalayan region, which is especially susceptible to the impacts of climate change. Specifically, research indicates that the mountainous areas of Nepal are experiencing a significant warming rate of 0.057 °C per year [2]. Additionally, studies on the Hindu Kush Himalayan region and the Tibetan Plateau have recorded a remarkable warming rate of 0.195 °C per decade from 1951 to 2014, making these regions among the most severely affected by global warming [3]. Due to this climate change, there has been a rise in extreme phenomena such as floods and droughts, which could have severe consequences for various socioeconomic sectors [4,5]. The recent indicators of climate change have markedly impacted glaciers and water resources in the Himalayas, leading to substantial effects on river discharge and groundwater replenishment, as highlighted in studies [6,7]. Alterations in riverflow from the Himalayan rivers result in significant and enduring effects on the water resources vital for the sustenance of communities in and around these areas, as well as downstream, as indicated by [7,8,9,10,11]. Hence, measuring how climate change affects extreme precipitation and the resulting extreme flow events is crucial. This is vital for creating more effective adaptation strategies, improved mitigation measures, and enhanced resilience to water-related hazards.

The Karnali River Basin (KRB), situated in the western regions of Nepal, has been the focus of many studies previously investigating the effects of climate change on riverflows, with a particular emphasis on high-flow events [12,13]. The discharge in the outlet of the KRB increased by 10.9% in 2070–2099 under RCP 8.5 scenario [12] and by 1.6% in 2071–2095 at the Chisapani outlet of the KRB under the same scenario [13]. Furthermore, this hydroclimatic research has primarily utilized data from the third and fifth phases of the Coupled Model Intercomparison Project (CMIP). These studies found temperature to increase throughout the basin, precipitation to show high spatial variability [14], and river discharge to follow a similar pattern as precipitation with each other on a yearly and seasonal basis [14,15].

Typically, future climate projections are derived from the results of advanced tools known as GCMs, which depict atmospheric processes using specific mathematical formulations [14]. In hydrological studies, GCM outputs are valuable for providing essential data for hydrological models, which simulate the interplay between climate change and hydrological dynamics [16]. In the latest CMIP6 suite of GCMs, there has been a reduction in the uncertainty related to climate models, owing to improved knowledge of the fundamental physics governing atmospheric processes in the last few decades [15,17]. Compared to CMIP5, CMIP6 models have improved spatial resolution, improved physical parameterizations, and incorporated more physical processes and components [18]. CMIP6 uses Shared Socioeconomic Pathways (SSPs), which include various scenarios of future socioeconomic conditions with respect to population, economic, health, and education. Combining the SSPs of CMIP6 and the Representative Concentration Pathways (RCPs) of CMIP5 provides future climate analysis while considering both climatic and societal change [19]. The comparative studies in the Himalayan and nearby regions between CMIP5 and CMIP6 have found the latter to perform better. CMIP6 has been more sensitive to anthropogenic greenhouse gas emissions over the South Asian region [20], and CMIP6 represents the spatial pattern of precipitation in the Himalayan region better than CMIP5 [15]. Other studies on Himalayan and nearby regions have also reported CMIP6 models to capture the pattern and intensity of temperature and precipitation better than the previous versions of CMIP models [20,21]. SWAT is used to build a hydrological model to simulate future river discharge in various gauge locations. SWAT, a continuous-time, watershed-scale model [22], has been applied widely in the Himalayan region to simulate riverflow [12,13,23,24]. SWAT performed very well in simulating streamflow in the KRB [23]. Previous studies on multi-site calibrations found that the SWAT approach enhanced model performance more than a single-site calibration, which is a better-suited model for a large, spatially heterogeneous watershed such as the KRB [13,25]. Streamflow describes quickflow, interflow, and baseflow [26,27,28]; of these three flows, baseflow is not associated with a single rainfall event [29] and sustains the longest between rainfall events, making baseflow significantly crucial for irrigating agricultural lands during dry seasons. Despite Karnali province’s highest GDP contribution from agriculture and forest sectors, minimal studies have been conducted, and none are from the KRB region [30]. Further, mountainous river basins are highly susceptible to the impacts of climate change. Previous studies on the KRB have predicted an increase in river discharge in the future. However, these forecasts were made based on a hydrological model calibrated at a single location and integrated with climate variable projections from CMIP5 [12]. A multi-site calibration would provide a better understanding of streamflow dynamics, particularly in geographically diverse basins such as the KRB.

Our study aims to understand how the hydrological cycle is affected by climate change in river basins. This study combines the multi-site calibrated SWAT model with the latest CMIP6 GCMs to quantitatively evaluate climate change’s impact on hydrology and the segregation of baseflow using the Web-based Hydrograph Analysis Tool (WHAT), which has been found to be suitable in Bagmati and Koshi Basins (Central Nepal) [31]. Our methodology is structured into five main parts: (i) selecting four suitable GCMs each for precipitation, maximum temperature, and minimum temperature; (ii) choosing the right bias correction method for each of these climate variables; (iii) conducting multi-site calibration and validation of the SWAT hydrological model, (iv) projecting future streamflow at multiple sites, and (v) separating baseflow with the help of the WHAT.

This paper makes significant progress in assessing the impact of climate change on the hydrology of high-altitude river basins like the KRB using a CMIP6 multi-model ensemble approach. It offers new insights into how temperature changes, especially the pronounced increase in minimum temperatures, impact these sensitive areas. Our findings reveal that rising temperatures lead to increased snowmelt and higher river discharge levels, altering water availability and the timing of peak discharges. Such changes have important implications for water resource management, agriculture, and disaster risk reduction. This research enhances our ability to address the challenges climate change poses on water resources by providing a robust framework for predicting hydrological changes. It presents comprehensive assessments for future scenarios, aiding policymakers, stakeholders, and researchers in sustainable water management strategies. Ultimately, this study equips local authorities with the knowledge to better prepare for and adapt to climate-induced changes, offering valuable insights for mitigating flood and drought risks in critical ecosystems.

2. Materials and Methods

2.1. Study Area

The Karnali River basin of western Nepal comprises a watershed area of 42,331 km² (Source: Watershed area calculation in GIS). The KRB extends from Dhaulagiri Mountain in the east to Nanda Devi Mountain in the west.

Though a small portion of the KRB lay in the Tibetan region, this research focused only on the basin within Nepal due to the inaccessibility of transboundary data. The KRB has six major watersheds: West Seti, Kawari, Humla Karnali, Mugu Karnali, Tila, and Bheri. They are all in Nepal except for Humla Karnali, which originates from China. The map of the KRB and meteorological stations used in this study are shown in Figure 1. The elevation of the KRB ranges from 200 to 7720 m above mean sea level. KRB lies within 28.33–30.45° N and 80.55–83.68° E. In the basin, 16% and 33% of the total area are occupied by agricultural land and forest, respectively. The Karnali River Basin has 742 glacial lakes and 1459 glaciers [32]. Based on the modified Köppen–Geiger classification, the northern mountainous region has a polar tundra climate. In contrast, the more level topographic region in the south has a temperate climate with a dry winter and hot summer [33]. Approximately 2% of the basin’s area is covered by permanent glaciers, and about 12% has snow and ice cover. The annual average precipitation in the basin is 1479 mm, and 77% occurs from June to September. The average maximum and minimum temperatures are 25 °C and 13 °C [34].

2.2. Data Collection

Meteorological and Hydrological Data

The Department of Hydrology and Meteorology provided all the weather information needed for the study, including daily precipitation, maximum and minimum temperature, solar radiation, wind, and humidity between 1980 and 2014. This study used data from 14 precipitation stations and 10 temperature stations to forecast future precipitation and temperature patterns. Additionally, we gathered data from 17 humidity stations, 17 solar radiation stations, and 15 wind speed stations, which served as inputs for the SWAT model, alongside precipitation and temperature data. The collected data span from 1980 to 2014. For calibration and validation of the rainfall–runoff model, we employed data from five stream gauge stations within the KRB, specifically: Karnali River at Asaraghat (Q240), Karnali River at Benighat (Q250), Thulo Bheri River at Rimna (Q265), Bheri River at Samaijighat (Q269.5), and Karnali River at Chisapani (Q280). All data were obtained from the Department of Hydrology and Meteorology (DHM) in Nepal. Details of the selected meteorological station are presented in

Table 1. Meteorological station information.

SN	Station Name	Index No	District	Longitude (°)	Latitude (°)	Elevation (m)	Annual Precipitation (mm)	Annual Max. Temperature (°C)	Annual Min. Temperature (°C)
1	Dadeldhura	104	Dadeldhura	80.58	29.30	1879	-	33.2	−5
2	Pipalkot	201	Bajhang	80.84	29.61	1455	2133.7	-	-
3	Silgsdhi	203	Doti	80.98	29.26	1309	-	39	−0.5
4	Martadi	204	Bajura	81.48	29.45	1598	2220.6	-	-
5	Kola Gauna	214	Doti	80.70	29.12	1364	2316.6	-	-
6	Mangalsen	217	Achham	81.25	29.13	1310	1598.8	-	-
7	Rara	307	Mugu	82.08	29.54	2989	834.2	-	-
8	Dipal Gaun	310	Jumla	82.22	29.26	2422	919.2	34.9	−14
9	Simikot	311	Humla	81.81	29.97	2993	828.6	29.5	−17.5
10	Dunai	312	Dolpa	82.89	28.95	2098	375.9	36.3	−7
11	Pusma Camp	401	Surkhet	81.23	28.87	953	1738.1	39	0
12	Dailekh	402	Dailekh	81.70	28.83	1394	1958.6	39.6	0
13	Birendra Nagar	406	Surkhet	81.63	28.58	720	1737.6	41.8	0.5
14	Maina Gaun	418	Mugu	82.26	28.96	1913	1908.4	-	-
15	Rukumkot	501	Rukum	82.62	28.61	1568	1754.3	-	-
16	Chaurjhari Tar	513	Rukum	82.21	28.65	863	1259.3	42	0.5
17	Musikot	514	Rukum	82.46	28.61	1412	-	41.9	−0.5

.

2.3. CMIP6-Based GCMs Data

This research analyzed 14 GCMs under two SSP scenarios, SSP245 and SSP585, from the CMIP6 database. We acquired historical and future daily atmospheric data from these selected GCMs via the Earth System Grid Federation (ESGF) website (https://esgf-node.ipsl.upmc.fr/search/cmip6-ipsl/, accessed on 20 April 2022). This platform provided raw data from GCMs for historical and future scenarios, covering a range of atmospheric climate variables, for example, temperature, precipitation, relative humidity, and geopotential height at various time scales. Details about the model’s name, the responsible institution, and the resolutions of each GCM used in this study are presented in

Table 2. Information on CMIP6 models used in this study, including model name, institutions, and atmospheric resolution.

SN	Model	Institution	Resolution
1	ACCESS-CM2	Commonwealth Scientific and Industrial Research Organisation (CSIRO) and ACCESS (Australian Research Council Centre of Excellence for Climate System Science.	1.25° × 1.875°
2	ACCESS-ESMI	Commonwealth Scientific and Industrial Research Organisation (CSIRO) and ACCESS (Australian Research Council Centre of Excellence for Climate System Science	1.25° × 1.875°
3	BCC-CSM2-MR	Beijing Climate Center, Beijing	1.125° × 1.125°
4	EC-Earth3	EC-Earth Consortium	0.35° × 0.35°
5	FGOALS-f3-L	Chinese Academy of Sciences Flexible Global Ocean–Atmosphere–Land System model	1° × 1°
6	INM-CM4-8	Institute for Numerical Mathematics, Russia	2° × 1.5°
7	IPSL-CM6A-LR	Institute Pierre Simon Laplace (IPSL), Paris	2.5° × 1.27°
8	INM.INM-CM4-8	Institute for Numerical Mathematics, Russia	2° × 1.5°
9	INM-CM5-0	Institute for Numerical Mathematics, Russia	2° × 1.5°
10	MPI-ESM1-2-HR	Max Planck Institute for Meteorology (MPI-M), Germany	0.94° × 0.94°
11	MRI-ESM2-0	Meteorological Research Institute, Ibaraki, Japan	1.125° × 1.125°
12	MPI-ESM1-2-LR	Max Planck Institute for Meteorology (MPI-M), Germany	1.875° × 1.86°
13	MIROC6	Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Kanagawa	1.4° × 1.4°
14	NorESM2-MM	Norwegian Climate Center, Norway	2.5° × 1.89°

. Our analysis focused on precipitation, maximum temperature, and minimum temperature as key predictors.

2.4. Digital Elevation Model (DEM), Land Use, and Soil Data

The SWAT model necessitates the integration of various data types, including meteorological, hydrological, and land characterization datasets. A DEM, along with land use and soil data, are crucial to this integration. For this study, the DEM, crucial in delineating the watershed’s topography, was sourced from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Digital Elevation Model Version 2. This model, with a spatial resolution of 30 m × 30 m, was accessed from ASTER’s website (https://asterweb.jpl.nasa.gov/GDEM.ASP, accessed on 15 June 2022). The DEM’s utility extends to extracting key watershed features, including the catchment area, drainage layout, and the terrain’s slope, aspect, and length. It also aids in determining drainage properties like channel dimensions and slope.

Land use and land cover data reflect the Earth’s systemic characteristics and are imperative for planning, management, and hydrologic modeling. The land use map specific to this study was developed by the International Centre for Integrated Mountain Development (ICIMOD) in 2010 and was retrieved from ICIMOD’s resource database (https://rds.icimod.org/Home/DataDetail?metadataId=9224, accessed on 15 June 2022). The land use types in the KRB are represented in Figure 2. Furthermore, the SWAT model’s functionality hinges on detailed soil characteristics, encompassing texture, bulk density, organic carbon content, hydraulic conductivity, and water retention capabilities across various soil layers. These soil attributes pertinent to Nepal were obtained from the International Soil Reference and Information Centre (ISRIC), which is the World Soil Information’s Soil and Terrain (SOTER) database for Nepal, available at ISRIC’s website (http://www.isric.org, accessed on 15 June 2022) [35]. The primary soil types identified within the catchment area include Gelic LEPTOSOL, Eutric REGOSOLS, Humic CAMBISOLS, Eutric CAMBISOLS, and Chromic CAMBISOLS, which are represented in Figure 3.

The slope is one of the most important factors in the simulation of hydrological processes. ArcSWAT requires slope class/classes to discretize the sub-watersheds into different hydrologic response units (HRUs) by integrating land slope classes with different land use and soil types in the sub-watersheds. In this study, a slope map of the study area was prepared from the ASTER DEM and classified into three different classes. High relief and steep slopes characterize the Karnali basin. The elevation varies from a minimum of 136 m above mean sea level (msl) to a maximum of 7707 m msl within the span of about 300 km only. Most of the area has is within a 25–75% slope range. However, about 13.5% of the area has a slope in the 75–100% range, and about 4.7% of the area has a slope even greater than 100%.

2.5. Methodology

The methodological flow chart is shown in Figure 4.

2.5.1. Robust Selection Criteria for GCMs

The selection of GCMs was strategically informed by their relevance to climate change research in the South Asian region. This study selected thirteen GCMs based on the availability of daily data for precipitation and maximum and minimum temperatures under historical conditions and four future scenarios (SSP126, SSP245, SSP370, SSP585), applying bias correction to their outputs. The study found that South Asia will experience a warmer (3–5 °C) and wetter (13–30% increase) climate in the 21st century [36].

In our study, nine GCMs were chosen for daily precipitation data and ten were selected for maximum and minimum temperature, covering both historical periods and future scenarios under SSP245 and SSP585. The historical period, spanning from 1980 to 2014, provided a baseline for assessing each GCM’s performance, utilizing daily mean values of the respective variables. To evaluate the performance of these models, four established criteria were employed, namely the Root Mean Square Error—observations Standard Deviation Ratio (RSR), the coefficient of determination (R²), the Percent Bias (PBIAS), and the Nash–Sutcliffe Efficiency (NSE). These criteria formed the basis of a performance rating system for each GCM [37]. The overall ranking of the GCMs was determined by averaging their ratings across all stations. In the subsequent phase, the top four GCMs for precipitation were selected from the initial pool of nine. The leading four GCMs for maximum and minimum temperatures were chosen from the ten available models. These models were selected based on the cumulative average of their performance ratings, using the four criteria mentioned earlier.

The detailed scoring of each GCM based on these performance metrics is presented in

Table 3. Rating system of performance metrics [37].

Rating	NSR	RSR	PBIAS	R2	Ratings
Very Good	0.75–1.00	0.00–0.50	<10	75–100	5
Good	0.55–0.75	0.50–0.6	15-Oct	65–75	4
Satisfactory	0.40–0.55	0.60–0.70	15–25	50–65	3
Unsatisfactory	0.25–0.40	0.70–0.80	25–35	40–50	2
Poor	≤0.25	>0.80	≥35	<40	1

. The final step involved selecting the four highest-rated GCMs for each variable—precipitation, maximum temperature, and minimum temperature. These models were then subjected to further bias correction and integrated into a multi-model ensemble, enhancing the robustness and reliability of the climate projections.

2.5.2. Bias Correction

Bias correction is pivotal in reducing systematic errors inherent in climate models, thereby improving their congruence with observational data. This technique modifies the daily outputs of GCMs by rectifying differences in mean and variability compared to observed data from the reference period (1980–2014) [38,39]. Various techniques have been employed for the bias correction of precipitation and temperature data. A study conducted in the Himalayan region of India demonstrated that quantile mapping effectively minimized biases in both precipitation and temperature [40]. Similarly, research in the Karkheh River basin, an area with diverse topography, showed that several quantile mapping methods successfully reduced biases in these climatic variables [41]. Our study area, also situated at a high altitude, employs quantile mapping for bias correction of raw GCMs.

Quantile Mapping (QM) can be classified into Parametric Transformation (PT), Non-Parametric Transformation (NPT), and Distribution Derivation Transformation (DDT). In this present research, 13 different Quantile Mapping (QM) bias correction methods were employed for adjusting precipitation in GCMs and 9 for maximum and minimum temperatures. Of the thirteen QM techniques used, four fall under DDT (Bernoulli Exponential, Bernoulli Gamma, Bernoulli Weibull, Bernoulli Log-normal); four under PT (Parameter Transformation function—exponential asymptote, Parameter Transformation function—linear, Parameter Transformation function—power, Parameter Transformation function—scale); and five under NPT (Non-parametric quantile mapping using empirical quantiles—linear, non-parametric quantile mapping using empirical quantiles—tricub, Non-parametric quantile mapping using robust empirical quantiles—linear, Non-parametric quantile mapping using robust empirical quantiles—tricub, Quantile mapping using a smoothing spline).

In this study, X₀ and X_m signify the observed and modeled climatic variables, such as precipitation, maximum temperature, and minimum temperature. To articulate the relationship between these variables, various formulations are used:

The general expression representing the transformation of modeled data into observed data:

$$X_0=h\left(X_m\right)$$

(1)

The inverse cumulative distribution function, mapping modeled data to its equivalent in the observed dataset:

$$X_0=F_0^{-1}\left(F_m\right(X_m\left)\right)$$

(2)

The power transformation is used to modify the modeled data by raising it to a power, scaled by a factor:

$$X_0=b{X}_m^c$$

(3)

A linear approach is depicted by:

$$X_0=a+b{X}_m$$

(4)

where the observed data is a linear function of the modeled data.

Exponential tendency of an asymptote is:

$$X_0=(a+b{X}_m)(1-e^{{-X}_m/r})$$

(5)

Simple scaling can be defined as:

$$X_0=b{X}_m$$

(6)

Power transformation with parameter x

$$X_0=b{(X_m-x0)}^c$$

(7)

Similarly, an exponential tendency to an asymptote with a parameter x is given by:

$$X_0=(a+b{X}_m)(1-e^{{-(X}_m-x0)/r})$$

(8)

These formulas, involving parameters a, b, c, x, and r, are subject to calibration to fit the specific characteristics of the climatic data under study.

To assess the performance of these bias correction methods, each of these bias correction methods was applied to four selected raw GCMs, and the performance of each method was evaluated based on four established criteria: RSR, R², PBIAS, and NSE [37]. These criteria constitute a performance rating system for each bias correction method. Table 3 presents the performance rating criterion for the performance indicators.

2.5.3. Formation of Multi-Model Ensemble (MME)

The reduction in uncertainties in forecasting a specific climate variable could be significantly enhanced by employing the mean time series derived from an MME of high-performing GCMs [40]. The literature presents various methodologies for deriving mean time series from an ensemble of select, better-performing GCMs. These methodologies range from straightforward arithmetic means to more complex machine learning algorithms [42]. Our current study used the simple mean approach to calculate the mean time series, which was specifically applied to variables such as precipitation, maximum temperature, and minimum temperature. These calculations were based on an ensemble comprising the four top-ranked GCMs. This choice of a simple mean method was based on its simplicity rather than other available techniques, and was grounded in its effectiveness and clarity in aggregating outputs from multiple models, thereby providing a consolidated and comprehensible representation of these key climate variables. Integrating these top-performing models aims to yield a more reliable and accurate mean time series, enhancing the robustness of our climate projections.

$$\mathrm{Simple} \mathrm{Mean}=\frac{1}{n}{\sum }_{j=1}^n{GCM}_j$$

(9)

where n refers to the number of GCMs considered for the development of MMEs, which is four in the present study, and GCM_j refers to the simulations of the climate variable of interest (i.e., precipitation, maximum temperature, and minimum temperature) produced by the j^th GCM.

2.5.4. Soil and Water Assessment Tool (SWAT)

In watershed hydrology, the SWAT model is pivotal for dissecting a watershed into smaller sub-watersheds, further delineated into hydrological response units (HRUs). These HRUs are integral to accurately representing the basin’s physical diversity, as distinct combinations of land use patterns, soil types, and slope gradients characterize them. This setup facilitated the calculation of soil water balance within each HRU, providing an in-depth view of the hydrological processes at play. The flow and water quality parameters were then methodically channeled from these HRUs to the sub-watersheds and eventually to the watershed’s main outlet. This meticulous representation of soil water balance in each HRU has been well documented in previous studies [22,43].

Furthermore, the SWAT model encompasses a comprehensive range of hydrological processes, including infiltration, percolation, evaporation, plant uptake, lateral flows, groundwater flows, snowfall, and snowmelt dynamics [43]. The model also simulated lateral flow using the kinematic storage model and calculated return flow by considering a shallow aquifer [22]. Additionally, channel flood routing was addressed using the Muskingum method, with necessary adjustments for transmission losses, evaporation, and return flow in the outflow estimation from a channel [44]. The water balance equation is central in the SWAT model [42,43], which integrates these varied hydrological components for a cohesive understanding of watershed hydrology:

$${\mathrm{S}\mathrm{W}}_{\mathrm{t}}={\mathrm{S}\mathrm{W}}_0+{\sum }_{\mathrm{j}=1}^{\mathrm{t}}({\mathrm{R}}_{\mathrm{d}\mathrm{a}\mathrm{y}}-{\mathrm{Q}}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}-{\mathrm{E}}_{\mathrm{a}}-{\mathrm{W}}_{\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{p}}-{\mathrm{Q}}_{\mathrm{g}\mathrm{w}})$$

(10)

where ${\mathrm{S}\mathrm{W}}_{\mathrm{t}}$ is the final soil water content (mm); ${\mathrm{S}\mathrm{W}}_0$ is the initial soil water content on day j (mm); ${\mathrm{R}}_{\mathrm{d}\mathrm{a}\mathrm{y}}$ is the amount of precipitation on day j (mm); ${\mathrm{Q}}_{\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}}$ is the amount of surface runoff on day j (mm); ${\mathrm{E}}_{\mathrm{a}}$ is the amount of evapotranspiration (ET) on day j (mm); ${\mathrm{W}}_{\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{p}}$ is the amount of water entering the vadose zone from the soil profile on day j (mm); and ${\mathrm{Q}}_{\mathrm{g}\mathrm{w}}$ is the amount of return flow on day j (mm).

2.5.5. Application of SWAT

The model has been designed to simulate daily streamflow. The ArcGIS 10.1 platform’s SWAT Model was created using Arc SWAT 2012. The basin area was divided into 45 sub-basins, which were further sub-divided into 245 HRUs. A threshold of 10% was set for slope, soil type, and land use. From 1995 to 2003 and 2004 to 2008, respectively, were used for calibration and validation for Chisapani, and from 1993 to 1994, there was a warming-up period. From 2001 to 2006 and 2007 to 2010 were used for calibration and validation for all the other four-gauge stations (Asaraghat, Benighat, Samaijighat, and Rimna), and from 1999 to 2000, there was a warming-up period. Model calibration is necessary for any hydrologic simulation to be used correctly and for the first testing and fine-tuning of observed data. The SWAT model for runoff was calibrated using daily observed runoff data collected at the watershed outlet from 1995 to 2003 for Chisapani and 2001 to 2006 for Asaraghat, Benighat, Samaijighat, and Rimna.

2.5.6. Multi-Site Calibration and Validation Procedure

Model calibration and uncertainty analysis were conducted using the semi-automated Sequential Uncertainty Fitting (SUFI-2) algorithm [45]. This algorithm employs a global search method via Latin hypercube sampling to determine optimal parameter ranges and conduct global sensitivity analysis. A multi-site and multi-variable calibration approach was applied in this study to better represent spatial heterogeneity in the KRB. The calibration and validation were first performed at upstream stations and then gradually moved toward downstream stations, as shown in

Table 4. Sub-basin number, HRU number, calibration, and validation period for different gauge stations.

Watershed	Sub-Watershed	Area (km2)	No. of Sub-Basin	No. of HRUs	Warm-Up Period (yrs)	Calibration Period	Validation Period
KRB	Rimna	2712.70	1	7	2	2001–2006	2007–2010
	Samaijighat	12,615.75	11	66	2	2001–2006	2007–2010
	Benighat	19,467.94	21	120	2	2001–2006	2007–2010
	Asaraghat	17,668.20	18	108	2	2001–2006	2007–2010
	Chisapani	42,086.67	45	248	2	1995–2003	2004–2008

.

The SWAT model parameters were varied at different spatial levels: HRUs, sub-basins, and basins. Since sub-catchments may have different basin characteristics, assigning the same basin parameter values to the whole catchment may limit the calibration process. To allow each sub-catchment to have its specific basin parameter values and avoid the limitation of calibration process, instead of fixing the parameter set as discussed above, we used the simulated flow obtained from that parameter set as input for the downstream station. For example, the output resulting from parameter set 1 at Rimna was used as input for the downstream Chisapani station, as shown in Figure 5. In this way, the basin parameter value obtained at the upstream station was not changed while calibrating for the downstream station. The same number of initial parameters and their ranges were used to start the calibration at each station.

2.5.7. Model Performance Evaluation

In this study, the performance of the SUFI2 model was evaluated using several indicators. The primary measure used for calibration and validation was the Nash–Sutcliffe (NS) coefficient [46]. Additionally, the R², PBIAS [47], and the ratio of the RSR served as supplementary criteria for assessment. Equations (11)–(14) were employed to compute these performance indices. The R² value indicates the proportion of variance in the observed data explained by the model simulations. The NSE measures the degree to which the plots of observed versus simulated data were aligned with a 1:1 line. The RMSE was derived from the variance in the modeled values plus the square of the bias; a lower RMSE signifies better model performance. An RMSE value of 0.0 indicates perfect replication of observed streamflow. PBIAS evaluates the general tendency of the simulated data to be larger or smaller than the observed data:

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{(P_i-O_i)}^2}{n}}$$

(11)

$$NSE=\frac{{\sum }_{i=1}^n{(O_i-\overline{O})}^2-{\sum }_{i=1}^n{(P_i-O_i)}^2}{{\sum }_{i=1}^n{(O_i-\overline{O})}^2}$$

(12)

$$R^2=\left[\frac{{\sum }_{i=1}^n(O_i-\overline{O})(P_i-\overline{P})}{\sqrt{{\sum }_{i=1}^n{(O_i-\overline{O})}^2}\sqrt{{\sum }_{i=1}^n{(P_i-\overline{P})}^2}}\right]$$

(13)

$$PBIAS=\left[\frac{{\sum }_{i=1}^n\left(O_i-P_i\right)\ast 100}{{\sum }_{i=1}^n\left(O_i\right)}\right]$$

(14)

where $O_i$ is the i^th observed value, $\overline{O}$ is the mean observed value, $P_i$ is the i^th model predicted value, $\overline{P}$ is the mean model predicted value, and n is the number of observations in the period under consideration.

2.6. Future Precipitation, Maximum Temperature, Minimum Temperature, and Streamflow Projection

The selected bias correction method was applied to the top four performing GCMs of each variable for future periods (2022–2100) to project precipitation, maximum temperature, and minimum temperature. We used a multi-site calibrated and validated SWAT model to evaluate future streamflow patterns. To study the temporal trends of future streamflow, we segmented the future timeframe into three distinct periods: Near Future (NF, 2022–2047), Mid Future (MF, 2048–2073), and Far Future (FF, 2074–2100). This integration was performed to forecast future streamflow at multiple sites (i.e., five gauge stations within the KRB) under the SSP245 and SSP585 scenarios to analyze spatial variations in streamflow. The projected streamflow from the SWAT model was subsequently compared with the baseline observed streamflow to assess streamflow change.

2.7. Baseflow Separation

Baseflow is the portion of streamflow that contains delayed sub-surface flow, generally maintained by groundwater discharge. Most streamflow during the dry season is accounted for by baseflow, which is vital to aquatic ecosystems and water quality. Without field data, many non-tracer-based techniques have been developed to estimate baseflow from streamflow, such as graphical analysis methods, numerical simulation methods, and digital filter methods [48]:

$$q_t=\alpha \times q_{t-1}+\frac{(1+\alpha )}{2}\times (Q_t-Q_{t-1})$$

(15)

The Web-based Hydrograph Analysis Tool (WHAT) system was developed by incorporating digital filter methods (Equations (15) and (16)) [27,49,50] for baseflow separation with the iSep system [51]. We employed a simple filter through a WHAT to distinguish baseflow from streamflow. This tool was chosen as the optimal baseflow separation technology for river basins within similar geographical locations [31]:

$$b_t=\frac{(1-{BFI}_{max}\times \alpha +b_t-1\times (1-\alpha )\times {BFI}_{max}\times Q_t}{1-\alpha \times {BFI}_{max}}$$

(16)

where α is the filter parameter, Q_t is the total streamflow at the t time step, b_t is the filtered baseflow at the t time step, b_t−1 is the filtered baseflow at the t − 1 time step, and the maximum value of the long-term ratio of baseflow to total streamflow is called the BFImax.

The three baseflow separation modules of the WHAT are the local minimum method (LLM), BFLOW filter, and Eckhardt filter. Among these three, we made a use of Eckhardt filter. Users must supply BFI_max (the maximum value of the long-term ratio of baseflow to total streamflow) and filter parameters to use the Eckhardt filter, which considers aquifer characteristics. Baseflow from streamflow was separated by the filter using Equation (16). Because a poor choice in BFImax can produce inaccurate baseflow separation findings, this approach predefines sample BFImax values for a range of hydrological and hydrogeological scenarios. The recommended values for this approach were 0.80 for perennial and ephemeral streams with porous aquifers and 0.25 for perennial streams with hard rock aquifers. The BFImax, or the value for perennial streams with permeable aquifers, was considered to be 0.80 [50].

3. Results and Discussion

3.1. GCM Selection

The CMIP6 GCMs were compared with the observed data to evaluate the ability of the GCMs to replicate the observed data and rate them based on performance metrics. Out of ten GCMs, six were discarded due to their high biases while comparing with the historical observed data (1980–2014). The four highest-performing GCMs for precipitation identified from this study were INM-CM5-0, INM-CM4-8, MPI-ESM1-2-LR, and ACCESS-ESM1-5. For maximum temperature, the top models found were INM-CM5-0, ACCESS-CM2, NorESM2-MM, and INM-CM4-8, while for minimum temperature, the leading models were MPI-ESM1-2-LR, NorESM2-MM, ACCESS-CM2, and MRI-ESM2-0. Notably, MPI-ESM1-2-HR secured the top rank for both maximum and minimum temperatures, and FGOALS-g3 was ranked fourth for maximum temperature. However, these models were not selected for the final analysis because their datasets were not extended to 2100. Consequently, they were substituted by the next highest-rated GCMs: INM-CM4-8 and INM-CM5-0 for maximum temperature and MPI-ESM1-2-LR for minimum temperature (Table 4), which also agreed with a study that reported that MPI-ESMI-2-LR and INM-CM5-0 would be suitable in simulating both precipitation and temperature over South Asia [51], and MPI-ESM1-2-HR would be suitable for simulating temperature in the Upper Indus Basin (Western Himalayas) [52]. The performance ratings for each GCMs, regarding precipitation, maximum temperature, and minimum temperature, are systematically tabulated in

Table 5. Performance rating of CMIP6 GCMs.

Precipitation	Rating	Max Temperature	Rating	Min Temperature	Rating
INM-CM5-0	2.464	MPI-ESM1-2-HR	3.050	MPI-ESM1-2-HR	3.575
INM-CM4-8	2.446	ACCESS-CM2	2.700	NorESM2-MM	3.475
MPI-ESMI-2-LR	2.339	NorESM2-MM	2.650	ACCESS-CM2	3.300
ACCESS-ESM1-5	2.321	FGOALS-g3	2.550	MRI-ESM2-0	3.300
BCC-CSM2-MR	2.196	INM-CM4-8	2.475	FGOALS-f3-L	3.025
MIROC6	2.125	INM-CM5-0	2.350	MPI-ESM1-2-LR	2.750
NorESM2-MM	1.875	BCC-CSM2-MR	2.350	INM-CM4-8	2.550
ACCESS-CM2	1.214	IPSL-CM6A-LR	2.125	INM-CM5-0	2.500
MRI-ESM2-0	1.214	MRI-ESM2-0	2.125	EC-Earth3	2.325
		EC-Earth3	1.000	IPSL-CM6A-LR	2.250

.

The results of our study reveal that the ratings for the top four highest-rated GCMs in simulating precipitation and maximum temperature were below 3, except for MPI-ESM1-2-HR for maximum temperature. In contrast, GCMs for minimum temperature received ratings between 3 and 4. A plausible reason for these lower ratings would have been attributed to the limited capacity of GCMs to accurately represent altitude-dependent climatic processes and other hydrological phenomena at a small basin scale, primarily due to their coarse resolution [53,54]. Consistent with these findings, previous studies have observed that GCMs are notably less effective in simulating precipitation and demonstrate marginally better performance in simulating minimum temperature than maximum temperature [21,55].

The top-performing models were selected based on their exceptional performance metrics for precipitation, maximum temperature, and minimum temperature across the KRB. Their accuracy and reliability in capturing the nuances of climate variability significantly contribute to the trustworthiness of climate projections in the region. Including these top performing GCMs in the KRB climate studies ensures more reliable and precise predictions of future climatic conditions. The reliability of these selected models strengthens the foundation for policymaking and planning processes by providing stakeholders with more confident forecasts, thereby enhancing resilience and sustainability in the face of climatic uncertainties.

3.2. Bias Correction Method Selection

Thirteen bias correction methods for precipitation and nine for maximum and minimum temperatures were implemented to correct bias for the raw simulated data from the four selected GCMs models. The study compared the bias-corrected GCM data with the observed data to evaluate how well it reproduces. The results showed that the bias-corrected GCM data are more consistent with the observed data than the raw GCM data. Figure 6, Figure 7 and Figure 8 illustrate how the selected bias correction methods bring the climatic variables closer to the observed values. The results from this present study showed that the best bias correction methods are (1) Bernoulli Weibull for precipitation, (2) Parametric transformation function—linear for maximum temperature, and (3) Smoothing spine for minimum temperature (

Table 6. Performance rating of the bias correction methods.

Bias Correction Methods	Rating
	Precipitation	Maximum Temperature	Minimum Temperature
Bernoulli Exponential	1.90
Bernoulli Gamma	2.43
Bernoulli Weibull	2.56
Bernoulli Log-normal	2.01
Non-parametric quantile mapping using empirical quantiles—linear	2.02	3.68	4.06
Non-parametric quantile mapping using empirical quantiles—tricub	2.08	3.71	4.06
Parameter Transformation function—exponential asymptote	1.72	4.33	4.01
Parameter Transformation function—linear	2.24	4.37	4.07
Parameter Transformation function—power	2.27	3.10	3.14
Parameter Transformation function—scale	1.56	3.96	3.43
Non-parametric quantile mapping using robust empirical quantiles—linear	2.05	3.68	4.06
Non-parametric quantile mapping using robust empirical quantiles—tricub	2.10	3.75	4.09
Quantile mapping using a smoothing spline	2.24	4.168	4.20

). This is consistent with a study conducted on the Bagmati irrigation system [56], which is nearest to the KRB, that reports the Bernoulli Weibull method as the best bias correction method for CMIP6-based GCMs precipitation. Our finding is also consistent with the studies conducted on the Kankai and Tamor River Basins in Nepal [24].

Through the application of these methods, we observed a significant reduction in the discrepancies between the raw data obtained from GCMs and the actual observed climatic variables. The reduction in bias makes the climate models more reliable and robust, providing a more accurate analysis of streamflow for future use. The effectiveness of these bias correction methods lies in their ability to mitigate uncertainties associated with climatic variables when they are applied in the context of basin hydrology. This reliability enhancement and uncertainty reduction are crucial for water resource management and environmental planning. They provide a more solid foundation for making informed decisions, particularly in designing and implementing adaptive strategies to mitigate the impacts of climate variability and change. The implications of these corrections are profound, enabling more accurate forecasting and planning, which is essential for sustainable water resource management and environmental conservation efforts.

Three representative stations were chosen to reflect the basin’s diverse topography: a hilly region station (Station-310), a Himalayan region station (Station-311), and a plains region station (Station-406). In each case, after implementing bias correction methods, the model data are more consistent with the observed data than the raw data shown in Figure 6, Figure 7 and Figure 8.

The accuracy of the selected bias correction methods for all stations across the KRB are detailed in Figure S1.

3.3. Formation of MME and Projected Trend of Climatic Variables

This study meticulously selected the four top-performing GCMs for precipitation, maximum temperature, and minimum temperature. It integrated them into an MME to assess changes in precipitation and temperature against the baseline of historical observations. There are several things to consider while using a multi-model ensemble approach; for example, individual models exhibit higher uncertainties; therefore, the results derived from such models may not be reliable. Implementing an ensemble approach can reduce the uncertainties associated with these models [57]. Different GCMs represent different but not all the existing climatic processes. The resulting errors are assumed to be partly canceled by combining these models to form a multi-model ensemble (MME) [58]. Studies have found that MMEs represent the trend in precipitation and temperature better than a single model [51,59,60]. However, the selection of an appropriate number of models is essential. The ensembles that contain many poorly performing models give significantly different outputs than the ensembles of only well-performing models [51]. There is no clear guideline for determining the ideal number of general GCMs to create a multi-model ensemble that minimizes uncertainties in climate forecasts. However, it is acceptable to select between three and ten of the highest-ranked GCMs [61,62]. Our study focused on the top four ranked GCMs to develop multi-model ensembles for various climate variables. We generated daily time series for these ensembles under the CMIP6 historical experiment and the future emission scenarios SSP245 and SSP585. This was accomplished using a straightforward mean technique, calculating the average daily bias-corrected outputs from these top four GCMs [61].

In this study, we opted three representative stations strategically to reflect the basin’s diverse topography: a hilly region station (Station-310), a Himalayan region station (Station-311), and a plains region station (Station-406). The findings indicate an upward trend in precipitation and maximum and minimum temperatures across all three regions under both SSPs, with increasing rates from the plains to the Himalayan region. Minimum temperature is anticipated to exhibit a more pronounced rise than maximum temperature. Our study utilized a multi-model ensemble of the four selected GCMs to project precipitation and maximum and minimum temperatures in the KRB for 2022–2100. The analysis indicates an increasing trend in precipitation and maximum and minimum temperatures over time, aligning with previous research in the KRB region [12,62]. Additionally, our results reveal a more rapid increase in minimum temperature than maximum temperature across the Himalayan, hilly, and plains regions. This trend is consistent with findings from other river basins in Nepal [63,64,65]. However, this contradicts the result of a study that reported a faster increase in maximum temperature in the Himalayan region [62]. Furthermore, our analysis suggests a higher trend in precipitation and maximum and minimum temperatures in highland areas compared to lowland regions. This observation corroborates the findings from [66,67]. This divergence highlights the complex and varied climatic patterns across different topographical zones within the region.

A detailed summary of the changes in precipitation and maximum and minimum temperatures for the three future periods under both SSP scenarios is compiled in

Table 7. Net change in projected precipitation, maximum temperature, and minimum temperature at station-310, -311, and -406 and an average for all stations (All).

SSP245	Precipitation (%)			Maximum Temperature (°C/yr)			Minimum Temperature (°C/yr)
Stations	NF	MF	FF	NF	MF	FF	NF	MF	FF
310	+4.73	+11.21	+13.72	+0.0003	+0.032	+0.054	+0.086	+0.117	+0.134
311	+1.67	+10.61	+14.47	+0.042	+0.073	+0.091	+0.107	+0.128	+0.14
406	−2.19	+5.05	+7.58	+0.0001	+0.027	+0.041	+0.039	+0.073	+0.092
All	+7.79	+11.65	+16.25	+0.018	+0.048	+0.064	+0.049	+0.08	+0.097
SSP585	Precipitation (%)			Maximum Temperature (°C/yr)			Minimum Temperature (°C/yr)
Stations	NF	MF	FF	NF	MF	FF	NF	MF	FF
310	+6.48	+13.49	+24.17	+0.003	+0.041	+0.097	+0.09	+0.153	+0.25
311	+6.12	+15.94	+28.78	+0.044	+0.088	+0.144	+0.112	+0.152	+0.201
406	+0.69	+8.41	+20.75	+0.024	+0.062	+0.109	+0.044	+0.103	+0.178
All	+9.43	+16.51	+27.47	+0.022	+0.066	+0.119	+0.057	+0.115	+0.187

.

Additionally, Figure 9, Figure 10 and Figure 11 depict the trend graphs for the three representative stations, while Figure S2 provides extended graphical data for the remaining stations.

As observed in our study, the differential increase in minimum temperatures has significant implications for water resources, plant life, and agriculture within a basin that spans from mountainous regions to plain areas. In mountainous regions, high minimum temperatures can cause a decrease in snowpack and earlier snowmelt, which can alter the timing and magnitude of runoff into rivers and streams. This shift can impact water availability for downstream ecosystems and human use, increasing the risk of droughts and floods in both highland and lowland areas. For plant life, warmer nights can affect physiological processes such as respiration [68], which could potentially alter growth patterns and the distribution in vegetation zones from the mountains to the plains. In agriculture, these temperature changes can disrupt the growing seasons, pest and disease patterns, and crop viability [69], necessitating adjustments in cropping systems, planting schedules, and water management practices. Consequently, it is crucial to adapt to these changes to sustain water resources and ensure food security throughout the regions. This requires integrated water resource management and adopting climate-resilient agricultural practices tailored to the diverse conditions found in mountainous and plains regions.

3.4. Seasonal Changes in Precipitation, Maximum Temperature, and Minimum Temperature

The vast geography of the KRB suggests a diverse spatial distribution in projected precipitation patterns. This study delves into this heterogeneity by analyzing future precipitation and maximum and minimum temperatures at three strategically chosen stations, each representing a distinct physiographic region within the basin. These include a mountain region station (Station-311, at 2993 m above sea level), a hilly region station (Station-310, at 2422 m above sea level), and a plains region station (Station-406, at 720 m above sea level). The study focuses on changes in precipitation and maximum and minimum temperatures, benchmarked against the baseline period average from 1980–2014. It examines variations in average annual precipitation, maximum temperature, and minimum temperature over the near, mid, and far future. Additionally, the study provides a comprehensive analysis of the rate of precipitation, maximum temperature, and minimum temperature changes throughout the year, dissected into the four distinct seasons: pre-monsoon (MAM), monsoon (JJAS), post-monsoon (ON), and winter (DJS).

3.4.1. Seasonal Changes in Precipitation across the KRB

Figure 12 illustrates the seasonal variation in precipitation across the mountain (Station-311), hilly (Station-310), and plains (Station-406) regions of the KRB, as forecast by an MME under both the SSP245 and SSP585 scenarios. In the mountain region (station-311), according to these scenarios, there will be a progressive increase in precipitation over time for all seasons except the post-monsoon season. This increase is expected to be the highest in the far future, followed by the mid future, and then the near future, compared to the baseline period. These models consistently predict a general increase in precipitation during pre-monsoon and monsoon and decrease in other seasons. Under both SSP245 and SSP585 scenarios, the pre-monsoon season is expected to experience a progressive increase in precipitation, ranging from 24% to 68%, extending from the near to the far future. Similarly, there is a projected increase in precipitation of 22% to 63% for the monsoon season, also following a progressive pattern. In contrast, precipitation is anticipated to decrease progressively by 58% to 63% during the post-monsoon season. Furthermore, the winter season will also gradually decrease in precipitation, ranging from 41 to 46%.

Compared to the mountainous areas, the hilly regions (Station-310) are also projected to show parallel trends in future precipitation changes over various periods, as per the multi-model ensemble under both SSP245 and SSP585 scenarios. Notably, the extent of precipitation change in these hilly terrains is anticipated to be less pronounced than in the mountainous areas. Throughout most seasons, except the post-monsoon and the near future of the pre-monsoon, a gradual increase in precipitation is forecast. The most significant rise is predicted for the far future, with moderately less in the mid future, and the least in the near future, relative to the baseline period. During the pre-monsoon season, precipitation changes are expected to vary between a decrease of 8% and an increase of 24%, progressively unfolding from the near to the far future. The monsoon season, in contrast, is likely to see an increase in precipitation ranging from 27 to 54%. Conversely, the post-monsoon season is set to experience a decline in precipitation by 11 to 27%. The winter season is forecast to undergo a modest alteration in precipitation, fluctuating between a 2% decrease and an 11% increase.

For the plains region (Station-406), as outlined in both SSP245 and SSP585 scenarios, a steady escalation in precipitation is predicted across all seasons, apart from the post-monsoon period. This trend is projected to intensify most substantially in the far future, followed by moderate increases in the mid future, and the smallest increments shortly when contrasted with the baseline period. In the pre-monsoon season, an anticipated progressive shift in precipitation is forecast, ranging from an increase of 60% to 117%, extending from the near future toward the far future. For the monsoon season, a varied change in precipitation is expected, with projections ranging from a decrease of 3% to an increase of 19%. Additionally, during the post-monsoon season, a consistent decrease in precipitation is anticipated, falling between 10 and 31%. Furthermore, the winter season is predicted to witness a steady rise in precipitation, estimated to be between 39 and 50%. Our findings align with prior studies on the Karnali–Mohana Basin [13] and the Gandaki River Basin [64], indicating an expected increase in precipitation during the pre-monsoon season.

3.4.2. Seasonal Changes in Maximum Temperature across the KRB

Figure 13 displays the seasonal variation in maximum temperature across the mountainous (Station-311), hilly (Station-310), and plains (Station-406) regions of the KRB, as predicted by a multi-model ensemble under both the SSP245 and SSP585 scenarios. In the mountain region (Station-311), a consistent increase in maximum temperature is projected for all seasons. This increase is expected to be most pronounced in the far future, followed by a moderate rise in the mid future and the least significant increase in the near future, compared to the baseline period. Under both SSP245 and SSP585 scenarios, the pre-monsoon season is anticipated to experience a steady rise in maximum temperature, ranging from 1.9 to 5.4 °C, progressively from the near to the far future. The monsoon season is also expected to follow a similar trend, with maximum temperatures increasing from 0.5 to 2.4 °C. During the post-monsoon season, a consistent increase in maximum temperature is forecast, varying between 0.2 to 2.8 °C. Furthermore, the winter season is predicted to undergo gradual changes in maximum temperature, ranging from −0.1 to 4.7 °C.

In the hilly region (Station-310), a trend similar to that of the mountainous regions is observed, albeit with a slightly less significant change in magnitude. The scenarios project a consistent increase in maximum temperature for all seasons in this region. This rise is forecast to be the greatest in the far future, followed by a moderate increase in the mid future and a smaller escalation in the near future compared to the baseline period. According to the SSP245 and SSP585 scenarios, the pre-monsoon season is expected to see gradual changes in maximum temperature, ranging from −0.3 to 3 °C, progressing steadily from the near future to the far future. During the monsoon season, an upward trend in maximum temperatures is anticipated, increasing from 0.3 to 1.8 °C. The post-monsoon season is projected to experience a continuous rise in maximum temperature, with variations between 0.2 to 2 °C. Furthermore, the winter season in the hilly region is predicted to undergo a steady alteration in maximum temperature, expected to vary from −2.6 to 3.3 °C.

The plains region (Station-406) is forecast to experience a continual increase in maximum temperature throughout all seasons. The most significant rise is anticipated in the far future, with a smaller increase expected in the mid future and the smallest increment occurring soon, compared to the baseline period. According to both SSP245 and SSP585 scenarios, the pre-monsoon season in the plains region is projected to undergo a steady change in maximum temperature, from −1.7 to 1.2 °C, progressively from the near future to the far future. During the monsoon season, an increase in maximum temperature is likely, varying between 1.2 to 2.8 °C. The post-monsoon season is expected to consistently rise in maximum temperature within the range of 0.5 to 2.7 °C. Additionally, the winter season is predicted to gradually increase in maximum temperature, ranging from 0.8 to 4.7 °C. The CMIP5 and CMIP6 GCMs predict an increase in maximum temperature in the Karnali–Mohana Basin and Gandaki River Basin, respectively, across all seasons [13,64].

3.4.3. Seasonal Changes in Minimum Temperature across the KRB

Figure 14 displays the seasonal variation in minimum temperature across the mountainous (Station-311), hilly (Station-310), and plains (Station-406) regions of the KRB, as predicted by a multi-model ensemble under both the SSP245 and SSP585 scenarios. In the mountain region (Station-311), these scenarios suggest a continuous increase in minimum temperature throughout all seasons. The most substantial rise is anticipated in the far future, followed by a less pronounced increase in the mid future and the smallest elevation compared to the baseline period. According to both SSP245 and SSP585 scenarios, the pre-monsoon season will steadily increase in minimum temperature, varying from 1.2 to 4.2 °C, gradually progressing from the near future to the far future. A similar upward trend is expected in the monsoon season, with minimum temperatures increasing from 1.6 to 4.2 °C. The post-monsoon season is projected to experience a consistent rise in minimum temperature, ranging between 1.8 to 4.2 °C. Additionally, the winter season is forecast to increase in minimum temperature gradually from 6.5 to 8.3 °C.

A similar trend is observed in the hilly region (Station-310) but with a slightly lower magnitude of change than in the mountainous regions. These scenarios indicate a steady rise in minimum temperature across all seasons in the hilly region. The increase is expected to be most significant in the far future, with a moderate rise in the mid future and a smaller increase in the near future concerning the baseline period. For the pre-monsoon season, the SSP245 and SSP585 scenarios project a gradual increase in minimum temperature, ranging from 1.6 to 5.3 °C, progressively from the near future to the far future. During the monsoon season, minimum temperatures are expected to follow an upward trajectory, increasing from 1.4 to 8.4 °C. In the post-monsoon season, a continuous increase in minimum temperature is anticipated, with a range between 2.6 to 5.9 °C. The winter season in the hilly region is also set to experience a steady increase in minimum temperature, estimated to range from 3.6 to 5.4 °C.

Reflecting trends close to the hilly region but with a distinct pattern, the plains region (Station-406) is also projected to increase minimum temperature across all seasons continuously. The most significant increase is forecast for the far future, diminishing slightly in the mid future and showing the least rise in the near future when compared with the baseline period. Under both SSP245 and SSP585 scenarios, the pre-monsoon season in the plains region is predicted to see a steady rise in minimum temperature, varying from 0.6 to 4.8 °C, extending progressively from the near future to the far future. The monsoon season is likely to witness an increase in minimum temperatures, ranging from 0.8 to 3.7 °C. The post-monsoon season is projected to consistently increase in minimum temperature, within the 1.4 to 5.1 °C range. Furthermore, the winter season is forecast to experience a gradual rise in minimum temperature from 1.3 to 5.3 °C.

3.4.4. Overall Impacts from Seasonal Changes in Temperature and Precipitation Patterns

There is a general increase in precipitation throughout the year, except for the post-monsoon period, which shows a decrease. This trend is accompanied by a consistent rise in both maximum and minimum temperatures across all seasons. Notably, the rate of increase in minimum temperatures is higher than that for maximum temperatures. Although precipitation is decreasing during the post-monsoon season, streamflow is expected to increase throughout all seasons. This could be attributed to the faster melting of snow in the Himalayan regions. The relationship between temperature increase and snowmelt contributes to increased streamflows. This highlights the complex dynamics of how climate change impacts hydrological patterns in the regions. Agricultural systems in the basin are likely to experience more stress due to the changing weather patterns. Higher temperatures can affect crop growth cycles, reduce yields of temperature-sensitive crops, and increase the spread of pests and diseases, which can ultimately lead to food insecurity. Furthermore, the lack of a clear trend in precipitation, especially the uncertain changes in the post-monsoon season, makes water resource management and agriculture planning more challenging. This requires more adaptive and flexible approaches to cope with the increased variability and uncertainty. The two scenarios, namely SSP245 and SSP585, indicate a noticeable warming trend in the KRB. The maximum and minimum temperatures have been increasing annually, and it is happening in all seasons. The warming trend is more evident in the SSP585 scenario, which reflects a higher greenhouse gas emission pathway compared to SSP245. However, the precipitation patterns do not show any clear and consistent trend. They indicate a slower rate of increase in most seasons, except the post-monsoon season, which does not follow a distinct pattern of change. The temperature increases that have been observed in the KRB have significant implications for the vulnerability and resilience of the basin to climate change. The rise in temperature, particularly the faster increase in minimum temperatures, can have ecological and hydrological impacts. For example, higher minimum temperatures can accelerate the melting of glaciers in the mountainous regions of the basin, which can affect the long-term availability of water for downstream communities, agriculture, and hydropower production. The timing and magnitude of runoff may be altered, which could increase the risk of floods during the wet season and worsen water scarcity during the dry season.

To enhance the resilience of the KRB to the impacts of climate change, a multi-faceted approach is required. This approach includes implementing integrated water resource management practices to optimize the available water resources, adopting climate-resilient agricultural practices such as drought-resistant crop varieties and improved irrigation techniques. Furthermore, investing in infrastructure and community-based initiatives to reduce the risk of climate-induced disasters is crucial. Strengthening the capacity of local communities to adapt to climate change through education, awareness, and access to technology and resources will also be essential in building a more resilient Karnali River Basin. In addition, the study highlights the significance of upgrading infrastructure to withstand severe weather conditions and promoting community involvement and capacity development to equip local communities with the knowledge and resources to build resilience to climate change. This holistic approach is essential for protecting the environmental, economic, and social well-being of the KRB concerning climate change. Policymakers can encourage a proactive approach to climate resilience by informing local communities about the expected changes. This can help ensure that adaptation measures are customized to the unique cultural and geographical characteristics of the KRB. In doing so, we can protect the environment and livelihoods from the impending impacts of climate change.

Comprehensive details on the monthly and seasonal changes in precipitation, maximum temperature, and minimum temperature under the SSP245 and SSP585 scenarios are available in Figure S3.

3.5. Hydrological Model Evaluation

During the auto-calibration process using SWAT-CUP for the KRB, NSE was selected as the primary objective function. The choice is rooted in the wide acceptance of NSE as a pivotal metric for gauging the effectiveness of hydrological models. The selected model’s efficacy lies in representing the hydrograph’s overall congruence accurately [37,70]. A strategy involving calibration at multiple sites was implemented to represent the diversity more accurately within the basin. This involved calibrating the model at five different hydrological stations, using observed daily discharge as the benchmark. The calibration process was divided into three phases: (1) a sensitivity analysis, (2) an automatic calibration, and (3) a manual calibration. The sensitivity analysis was followed using SWAT-CUP for automated calibration, involving an initial run of 500 iterations to refine the values for the most sensitive parameters. After the automatic calibration, a manual calibration was carried out, considering the specific characteristics and knowledge of the basin. In total, calibration was conducted at five hydrological stations. The process was carried out sequentially from the uppermost part of the basin to the lowermost, with calibration settings being fixed at upstream points before proceeding to downstream stations. Previous studies have shown that the SWAT model effectively assesses the impact of climate change on streamflow [12,23]. However, these studies have relied on calibrations at single sites. Our research employed a multi-site calibration approach to improve model performance and reduce uncertainties. Additionally, previous studies have not explored baseflow separation in the KRB. We used a WHAT to filter the baseflow to address this gap. This method was appropriate for isolating the baseflow in areas near the KRB that have similar geographical characteristics [31].

3.5.1. Sensitivity Analysis

The study conducted a comprehensive sensitivity analysis to identify key parameters influencing hydrological processes within five distinct basins in the Karnali Basin, Nepal. Utilizing the SUFI-2 algorithm, this analysis examined the sensitivity and uncertainty inherent in the SWAT model. A set of 30 parameters, as recommended in the SWAT-CUP manual, was selected for the calibration process after careful examination of the SWAT literature [39,71,72]. Parameters displaying larger absolute t-stat values were identified as more sensitive, indicating their significant impact on the model. Additionally, parameters with p-values approaching zero were deemed statistically significant. Those with a p-value below 0.05 were particularly noted as sensitive, signifying their substantial influence on the model’s performance, a concept underscored by Abbaspour’s research in 2014.

Bheri River at Samaijighat (Q269.5)—GW_DELAY emerged as the most sensitive parameter, followed by CH_N1 and SOL_K. The high sensitivity of these parameters suggests they play a crucial role in controlling groundwater delay, channel roughness, and soil permeability. Conversely, SMFMX, CH_N2, and SLSUBBSN were the least sensitive parameters, indicating a relatively minor impact on snowmelt factors, channel roughness in sub-basins, and slope length.

Thulo Bheri River at Rimna (Q265)—TLAPS was identified as the most sensitive, succeeded by PLAPS and RCHRG_DP, indicating their significant roles in temperature lapse rate, precipitation lapse rate, and deep aquifer recharge. Being the least sensitive, CH_N1, CN2, and REVAPMN imply lesser influence on channel roughness, runoff curve number, and shallow aquifer revap coefficient.

Analysis of the Karnali River at Benighat (Q250) indicated that TLAPS, followed by PLAPS and REVAPMN, were the most sensitive parameters, suggesting a strong influence on temperature and precipitation patterns, and shallow aquifer dynamics. GWQMN, SMFMX, and SHALLST were the least sensitive, indicating a minor effect on groundwater minimum flow, snowmelt, and shallow aquifer storage. Karnali River at Asaraghat (Q240), TLAPS, PLAPS, and SFTMP were identified as the most sensitive parameters, emphasizing their critical role in temperature, precipitation variability, and soil temperature. Conversely, SLSUBBSN, ESCO, and SHALLST, being the least sensitive, imply a reduced impact on slope length, soil evaporation compensation, and shallow aquifer storage. Finally, for the Karnali River at Chisapani (Q280), GW_DELAY, TLAPS, and PLAPS were the most sensitive, highlighting their importance in groundwater processes and temperature and precipitation variations. In contrast, SOL_AWC, ALPHA_BF, and LAT_TTIME, as the least sensitive parameters, suggest a lesser role in available water capacity, baseflow, and lateral flow travel time. The study revealed a notable variation in the ranking of sensitive parameters across the different sub-watersheds within the KBR. This variation underscores the unique hydrological characteristics of each basin. The optimal parameter ranges derived from the watershed calibration, using identical initial parameter ranges, are detailed in

Table 8. Calibrated values of the SWAT parameters for the hydrological stations.

Parameter	Change Type	Suggested Ranges			Gauge Stations
			Asaraghat	Benighat	Rimna	Samaijighat	Chisapani
CN2	r	35–98	60.2	58.9	50.7	53.9	53.3
ALPHA_BF	v	0–1	0.12	0.31	0.28	0.19	0.23
GW_DELAY	v	0–500	189.11	152.68	42.87	57.56	51.39
GWQMN	v	0–5000	1.41	1.39	1.28	1.39	1.3
GW_REVAP	v	0.02–0.2	0.03	−0.02	0.04	0.13	0.04
ESCO	v	0–1	0.87	0.84	0.86	0.88	0.84
CH_N2	v	−0.01–0.3	0.29	0.29	0.32	0.32	0.32
CH_K2	v	−0.01–500	86.52	94.98	75.19	75.51	65.71
ALPHA_BNK	v	0–1	0.56	0.63	0.5	0.52	0.57
SOL_AWC	r	0–1	0.5	0.51	0.52	0.58	0.51
SOL_K	r	0–2000	0.19	0.28	0.09	0.11	0.29
SOL_BD	r	0.9–2.5	−1.56	−1.17	−1.48	−1.12	−1.01
HRU_SLP	r	0–1	0.25	0.23	0.24	0.31	0.31
OV_N	r	0.01–1	1.02	1.06	1	0.99	1.07
SLSUBBSN	r	10–150	72.79	65.32	62.6	68.57	70.9
REVAPMN	v	0–500	661.4	693.56	737.55	678.21	692.21
RCHRG_DP	a	0–1	−0.08	−0.06	−0.15	−0.09	−0.06
SHALLST	r	0–50,000	30,916.34	28,810.26	25,901.14	18,576.11	33,973.12
CANMX	r	0–100	86.71	90.27	80.96	81.07	94.55
EPCO	r	0–1	0.76	0.76	0.7	0.71	0.71
LAT_TTIME	r	0–180	−1.13	−0.5	25.81	−78.45	−21.28
CH_N1	r	0.01–30	4.28	10.62	5.19	1.84	−2.92
SFTMP	v	−20–20	0.57	2.42	-	0.42	1.06
SMTMP	v	−20–20	2.93	1.83	-	3.39	2.06
SMFMX	v	0–20	4.85	6.87	-	8.52	6.07
SMFMN	v	0–20	3.5	0.75	-	0.52	0.86
TIMP	v	0–1	0.08	0.02	-	−0.32	−0.3
PLAPS	v	−1000–1000	0.02	0.02	0.05	0.02	0.02
TLAPS	v	−10–10	−6.77	−6.42	−7.14	−7.86	−8.09
SURLAG	v	0.05–24	-	-	-	2.4	2.41

. This analysis provides a foundational understanding of effective watershed management and modeling within the KRB.

3.5.2. Multi-Site Model Calibration and Validation

For the Bheri River at Samaijighat (Q269.5), the period from 2001 to 2006 was utilized for calibration, while 2007 to 2010 was used for validation. The calibration results yielded an NSE of 0.75, a R² of 0.82, and a Percent Bias (PBIAS) of 28.1. These values fall into the categories of “Good” for NSE, “Very Good” for R², and “Unsatisfactory” for PBIAS [37,73]. During the validation phase, the model showed NSE, R², and PBIAS values of 0.69, 0.71, and −9.6, respectively, rated as “Good” for both NSE and R² and “Very Good” for PBIAS. The model exhibited mixed accuracy in simulating baseflow during the calibration phase, with instances of both overestimations and underestimations, while the validation phase indicated a general trend of overestimation. In terms of peak flows, these were slightly underestimated during calibration but showed improved accuracy in the validation phase.

Thulo Bheri River at the Rimna (Q265)—the period from 2001 to 2006 was designated for calibration, and 2007 to 2010 was set aside for validation. During calibration, the NSE, R², and PBIAS values were 0.77, 0.77, and −0.6, respectively and were classified as either “Good” or “Very Good” [37,73]. In the validation phase, the model exhibited NSE, R², and PBIAS values of 0.51, 0.54, and 16.5, respectively, all rated as “Satisfactory”. Baseflow estimation proved mostly accurate during the calibration phase. However, there was a noticeable underestimation of peak flows, especially in 2005 and 2006. The validation period demonstrated a decline in model performance, characterized by overestimations in baseflow and significant underestimations in peak flow, as illustrated in Figure 15. This decrease in performance is attributed to the less-than-optimal data quality from the Rimna gauge station.

Karnali River at Benighat (Q250)—the period from 2001 to 2006 was allocated for calibration, while 2007 to 2010 was used for validation. The calibration phase yielded NSE, R², and PBIAS values of 0.84, 0.84, and −0.1, respectively, each of which was rated as “Very Good” [37,73]. During the validation phase, the model demonstrated NSE, R², and PBIAS values of 0.75, 0.78, and −12.1, respectively, falling into the “Good” category. In this case, the model effectively simulates both baseflow and peak flows during the calibration period. However, a slight overestimation in baseflow was observed during validation, as indicated by the shift in the PBIAS score.

Karnali River at Asaraghat (Q240)—the period from 2001 to 2006, was designated for calibration, while 2007 to 2010 was used for model validation. During the calibration phase, the NSE, R², and PBIAS values were 0.81, 0.81, and −0.7, respectively, all categorized as “Very Good” [37,73]. In the validation phase, these values were 0.66, 0.71, and −12.5, respectively, and were rated as “Good”. The model effectively simulated baseflow during calibration, except for the years 2004, 2005, and 2006. However, during the validation period, baseflow was consistently overestimated. The simulation of peak flows was largely accurate throughout the calibration and validation periods. Finally, the Karnali River at Chisapani (Q280) was designated for calibration from 1995 to 2003, while 2004 to 2008 was set aside for validation. The calibration phase yielded NSE, R², and PBIAS values of 0.89, 0.89, and 1.9, respectively, each rated as “Very Good” [37,73]. In the validation phase, the model achieved NSE, R², and PBIAS values of 0.73, 0.74, and −5.7, respectively, classified as “Good” or “Very Good”. The model effectively simulates baseflow during both the calibration and validation periods. However, it tended to underestimate peak flows.

The model performance for the five gauge stations was examined by comparing the observation and simulated discharge, and is summarized in

Table 9. Performance of calibration and validation.

Station	Variable Period	NSE	R2	PBIAS
Asaraghat	Calibration period (2001–2006)	0.81	0.81	−0.7
	Validation period (2007–2010)	0.66	0.71	−12.5
Benighat	Calibration period (2001–2006)	0.84	0.84	−0.1
	Validation period (2007–2010)	0.75	0.78	−12.1
Rimna	Calibration period (2001–2006)	0.77	0.77	−0.6
	Validation period (2007–2010)	0.51	0.54	16.1
Samaijighat	Calibration period (2001–2006)	0.75	0.82	28.1
	Validation period (2007–2010)	0.69	0.71	−9.6
Chisapani	Calibration period (1995–2003)	0.89	0.89	1.9
	Validation period (2004–2008)	0.73	0.74	−5.7

. The scatter plot represented in Figure 16 shows a good correlation between the observed and simulated streamflow at all gauge stations within the KRB.

3.6. Impact of Climate Change on Streamflow

In Section 3.3, the MME provides projected data on precipitation, maximum temperature, and minimum temperature. These projections are integral inputs for the calibrated and validated SWAT model, enabling the projection of future streamflow at the five stream gauge sites. The monthly projections for discharge, alongside precipitation, maximum temperature, and minimum temperature for three future periods—Near Future (NF), Mid Future (MF), and Far Future (FF)—are presented under two scenarios, SSP245 and SSP585, as illustrated in Figure 17 and Figure 18, respectively. Significant increases in the discharge of various rivers are noted. For instance, under the SSP245 scenario, the discharge of the Bheri River at Samaijighat (Q269.5) is expected to rise by 86%, 99%, and 107% during the NF, MF, and FF periods, respectively. Under the SSP585 scenario, these increases are projected at 87%, 115%, and 144% for the same periods. Similarly, the Thulo Bheri River at Rimna (Q265) shows projected increases of 35%, 38%, and 43% during the NF, MF, and FF under SSP245, and 35%, 40%, and 49% under SSP585. The Karnali River exhibits notable increments as well. At Benighat (Q250), 79%, 97%, and 89% increases are projected under SSP245 for NF, MF, and FF, respectively, with even higher increases of 103%, 119%, and 135% under SSP585. At Asaraghat (Q240), projections indicate increases of 89%, 109%, and 132% under SSP245, and 117%, 135%, and 150% under SSP585. Lastly, the discharge at Chisapani (Q280) is expected to increase by 27%, 40%, and 54% during the NF, MF, and FF under SSP245 and by 45%, 70%, and 93% under SSP585. These projections suggest a significant overall increase in streamflow throughout the year. This trend is anticipated despite a decrease in post-monsoon precipitation, primarily due to heightened snowmelt because of rising temperatures. A previous study investigated how climate change affects water availability in the Karnali River Basin, utilizing CMIP5 model outputs and single-site calibration. The study found that the atmospheric temperature across the basin has increased, but there is significant variability in precipitation patterns depending on the location. Their finding anticipates that there will be an increase in annual precipitation and river discharge compared to baseline levels [21]. Our research, which leverages CMIP6 model outputs and multi-site calibration, reveals that minimum temperatures are rising more rapidly than maximum temperatures, with a uniform temperature increase across all monitoring stations. Precipitation patterns vary depending on the season and location. Historically, precipitation peaks at gauge stations have trailed river discharge peaks by approximately one month. This delay is attributed to snowmelt in the Himalayas, which sustains riverflows even after decreased precipitation. However, with projected increases in maximum and minimum temperatures, snowmelt may accelerate, leading to future peaks of precipitation and river discharge occurring within the same month [12].

3.7. Baseflow Separation

Recognizing the significance of baseflow is fundamental for effectively managing the river systems. Employing hydrograph analysis, this study differentiates baseflow from the varying streamflow patterns. Such differentiation is crucial for managing riverflows, maintaining water quantity and quality, ensuring reliable water supply, and preserving aquatic habitats. In this study, detailed baseflow separation was conducted at multiple locations, including the Karnali River at Asaraghat (Q240), Karnali River at Benighat (Q250), Thulo Bheri River at Rimna (Q265), Bheri River at Samaijighat (Q269.5), and Karnali River at Chisapani (Q280) using the WHAT. The seasonal BFI is systematically presented in Table 10, while Figure 19 and Figure 20 graphically depict the monthly variations in BFI, baseflow, and streamflow across each sub-basin, offering a comprehensive view of the hydrological patterns.

3.7.1. Observed

During the historical period, the lowest average monthly BFI is consistently observed in July across all watersheds. This index experiences a gradual increase from July, reaching its highest point in October. However, the Rimna watershed deviates from this pattern, achieving its highest BFI in December at 0.94. In the KRB, the average monthly BFI shows its most significant changes during the monsoon season, with values ranging from 0.62 to 0.91. This is followed by more stable periods: the pre-monsoon season with BFI ranging from 0.69 to 0.81, the winter season from 0.82 to 0.94, and the post-monsoon season showing the slightest variation, maintaining a consistent range of 0.82 to 0.94. Table 10 illustrates these BFI values across different seasons under the SSP245 and SSP585 scenarios.

Table 10. Seasonal trend in the BFI range for the KRB.

		Pre-Monsoon	Monsoon	Post-Monsoon	Winter
Historical		0.69–0.81	0.62–0.91	0.93–0.97	0.82–0.94
SSP245	Near	0.54–0.84	0.63–0.95	0.87–0.99	0.81–0.97
	Mid	0.60–0.82	0.65–0.94	0.87–0.99	0.82–0.96
	Far	0.62–0.81	0.66–0.97	0.87–0.99	0.80–0.97
SSP585	Near	0.57–0.83	0.63–0.95	0.88–0.99	0.81–0.97
	Mid	0.61–0.83	0.62–0.95	0.88–0.99	0.82–0.96
	Far	0.62–0.80	0.64–0.96	0.88–0.99	0.82–0.97

Table 10. Seasonal trend in the BFI range for the KRB.

		Pre-Monsoon	Monsoon	Post-Monsoon	Winter
Historical		0.69–0.81	0.62–0.91	0.93–0.97	0.82–0.94
SSP245	Near	0.54–0.84	0.63–0.95	0.87–0.99	0.81–0.97
	Mid	0.60–0.82	0.65–0.94	0.87–0.99	0.82–0.96
	Far	0.62–0.81	0.66–0.97	0.87–0.99	0.80–0.97
SSP585	Near	0.57–0.83	0.63–0.95	0.88–0.99	0.81–0.97
	Mid	0.61–0.83	0.62–0.95	0.88–0.99	0.82–0.96
	Far	0.62–0.80	0.64–0.96	0.88–0.99	0.82–0.97

3.7.2. SSP245

Figure 19 presents the monthly variation in streamflow, baseflow, and the BFI at five gauge stations under the SSP245 scenario. This scenario reveals a notable shift in the minimum average monthly BFI timing compared to the historical period, occurring one to three months earlier in future periods. Specifically, for Q240, the minimum average monthly BFI is observed in June at the NF station and May at the MF station, contrasting with the historical minimum in July. At Q250, this minimum occurs in June for both NF and MF and in May for FF. A similar trend is evident in Q265, with the minimum occurring in April across NF, MF, and FF. In the case of Q269.5, the minimum is observed in May at NF and April at MF and FF. Lastly, at Q280, the minimum is recorded for NF, MF, and FF in May. Regarding the peak average monthly BFI, it either aligns with the historical period or appears a month earlier. In Q250, Q269.5, and Q280, the peak is in October for NF, MF, and FF. For Q265, this peak is in November across NF, MF, and FF, while in Q240, it is observed in October at NF and in September at MF and FF. These data illustrate significant temporal shifts in BFI under the SSP245 scenario, indicating altered hydrological patterns in future periods. Within the KRB, the average monthly BFI exhibits its most substantial variations during the monsoon season. The BFI ranges from 0.63 to 0.95 in the NF region, 0.65 to 0.94 in MF, and 0.66 to 0.97 in FF. Following the monsoon season, the pre-monsoon period shows narrower ranges: 0.54 to 0.84 in NF, 0.60 to 0.82 in MF, and 0.62 to 0.81 in FF. In the winter season, the BFI stabilizes further, with values spanning from 0.81 to 0.97 in NF, 0.82 to 0.96 in MF, and 0.80 to 0.97 in FF. The least variation is observed during the post-monsoon season, where the BFI maintains a remarkably consistent range of 0.87 to 0.99 across NF, MF, and FF. This pattern of variability across seasons highlights the dynamic nature of hydrological responses in the KRB to seasonal changes.

Figure 19. Comparison of average monthly baseflow index, baseflow, and total streamflow at five gauge stations across the KRB under the baseline period and three different future periods (NN, MF, and FF) under the SSP245 scenario.

Figure 19. Comparison of average monthly baseflow index, baseflow, and total streamflow at five gauge stations across the KRB under the baseline period and three different future periods (NN, MF, and FF) under the SSP245 scenario.

3.7.3. SSP585

Figure 20 illustrates the monthly variations in streamflow, baseflow, and the BFI at five gauge stations under the SSP585 scenario. This scenario indicates a trend parallel to the SSP245 scenario, with the minimum average monthly BFI consistently occurring earlier than in the historical period. Specifically, for Q240, the minimum BFI is noted in June for NF and in May for MF and FF. In the case of Q250, this minimum is observed in June across all future scenarios. For the Rimna region, the minimum BFI shifts further ahead, recorded in April for NF and MF and even earlier in March for FF. For Q269.5, the minimum BFI is seen in May for NF and FF and in June for MF, while in Q280, it is in June for NF and MF and May for FF. The peak average monthly BFI either mirrors the historical trend in July or appears a month earlier. For Q250, Q269.5, Q280, and Q240, the peak is consistently observed in October across all future periods. At Q265, the peak manifests in November for all future scenarios. Within the KRB region, the average monthly BFI exhibits its most significant fluctuations during the monsoon season, ranging from 0.63 to 0.95 in NF, 0.62 to 0.95 in MF, and 0.64 to 0.96 in FF. This is followed by the pre-monsoon season, displaying BFI ranges of 0.57 to 0.83 in NF, 0.61 to 0.83 in MF, and 0.62 to 0.80 in FF. The winter season presents more stability in BFI, ranging from 0.81 to 0.97 in NF, 0.82 to 0.96 in MF, and 0.82 to 0.97 in FF. The post-monsoon season shows the most negligible variation, maintaining a consistent BFI range of 0.88 to 0.99 across all future periods. These findings underscore the impact of the SSP585 scenario on hydrological patterns, reflecting significant shifts in BFI across different seasons and regions.

Figure 20. Comparison of average monthly baseflow index, baseflow, and total streamflow at five gauge stations across the KRB under in baseline period and three different future periods (NN, MF, and FF) under the SSP585 scenario.

Figure 20. Comparison of average monthly baseflow index, baseflow, and total streamflow at five gauge stations across the KRB under in baseline period and three different future periods (NN, MF, and FF) under the SSP585 scenario.

Human activities and climatic factors, such as climate change and land use and cover alterations, can significantly influence streamflow and baseflow trends. BFI is subject to variability stemming from shifts in hydrological conditions, including changes in precipitation patterns, land usage, and the physical characteristics of watersheds [74,75]. Our study indicates that under the SSP245 and SSP585 scenarios, the minimum monthly BFI is projected to occur earlier than the historical period, typically in July. This shift could be attributed to increased streamflow from enhanced snowmelt, driven by rising minimum and maximum temperatures. Snowmelt contributes to streamflow predominantly as direct runoff, augmenting the proportion of direct runoff and consequently diminishing the BFI. Additionally, the higher elevation of gauge stations, particularly in the Rimna watershed, may contribute to its earlier projection of minimum monthly BFI compared to other watersheds. This analysis underscores the intricate interplay between climatic changes, land use patterns, and hydrological responses. More information on the baseflow separation analysis is in Tables S1 and S2.

Hence, climate change presents challenges and opportunities in the mountain, hilly and plains regions of the KRB. The KRB consists of various physiographic regions ranging from mountains to plains, and climate change introduces distinct challenges and opportunities in each region. The temperature increases are more pronounced at higher altitudes and diminish toward the plains, coupled with uncertain precipitation patterns, particularly a decrease in post-monsoon rainfall. These changes pose multifaceted impacts, such as glacier retreat affecting water availability, an increased risk of natural disasters such as floods and landslides, and shifts in agricultural zones impacting food security and biodiversity. Therefore, adaptation strategies must be region-specific. In the mountains, promoting climate-resilient agriculture and water conservation is vital. In the hilly regions, implementing soil and water conservation measures is necessary to prevent erosion and landslides. In the plains, enhancing flood management and irrigation efficiency is crucial. Building socioeconomic resilience across these regions involves integrating modern technologies with traditional practices, improving community-based monitoring and early warning systems, and fostering cross-sectoral approaches. This will ensure sustainable water resources management and agricultural practices, securing livelihoods against climate variability.

4. Conclusions

This study utilized a physically based, semi-distributed SWAT model for multi-site calibration across the KRB region. This approach was specifically chosen to assess accurately the impact of climate change on streamflow. To this end, we integrated outputs from CMIP6 GCMs into the SWAT model, providing a robust framework for our analysis. Our research evaluated two critical future scenarios: the high-emission (SSP585) and a more moderate, middle-of-the-road approach (SSP245). Additionally, these scenarios were examined across three distinct future time zones: the Near Future (NF, 2022–2047), the Mid Future (MF, 2048–2073), and the Far Future (FF, 2074–2100), offering a comprehensive view of potential outcomes. The significant findings of this study, which provide valuable insights into future water resource management and environmental planning, are summarized as follows: The top four highest-ranked GCMs based on performance metrics were INM-CM5-0, INM-CM4-8, MPI-ESMI-2-LR, and ACCESS-ESMI-5 for precipitation; INM-CM5-0, ACCESS-CM2, NorESM2-MM, and INM-CM5-0 for maximum temperature; and MPI-ESM1-2-LR, NorESM2-MM, ACCESS-CM2, and MRI-ESM2-0 for minimum temperature.

The Bernoulli Weibull—linear parametric transformation function and the smoothing spline were identified as the best bias correction methods for precipitation, maximum temperature, and minimum temperatures, respectively.

Annually, the study observed an increase in precipitation (+7 to +28%), a rise in maximum temperature (+0.018 to +0.118 °C/year), and a more pronounced increase in minimum temperature (+0.05 to +0.187 °C/year). Seasonally, under the SSP245 scenario, precipitation is expected to increase in the pre-monsoon (+17.4 to +59.2%), monsoon season (+10 to +15%), and winter seasons (+17.4 to +23.9%), and decrease in the post-monsoon period (−29.3 to −17.9%). Maximum temperature changes are more subtle, with a slight increase across seasons: +0.026 to +0.1 °C/year in pre-monsoon, +0.035 to +0.056 °C/year in monsoon, 0.007 to 0.017 °C/year in post-monsoon, and −0.005 to +0.07 °C/year in winter. Minimum temperature shows a consistent increase: +0.02 to +0.08 °C/year in pre-monsoon, +0.05 to +0.1 °C/year in monsoon, +0.034 to +0.058 °C/year in post-monsoon, and +0.08 to +0.12 °C/year in winter. Similarly, for precipitation under the SSP585 scenario, the study projects an increase in pre-monsoon (+19.9 to +57.3%), monsoon (+10.8 to +32.4%), and winter (+14.6 to +33.4%), and a decrease in post-monsoon (−42.7 to −31.2%). Maximum temperature is expected to change by −0.01 to +0.11 °C/year in pre-monsoon, +0.03 to +0.09 °C/year in monsoon, +0.03 to +0.1 °C/year in post-monsoon, and +0.04 to +0.17 °C/year in winter. For minimum temperature, increases are projected as follows: +0.03 to +0.17 °C/year in pre-monsoon, +0.06 to +0.21 °C/year in monsoon, +0.04 to +0.14 °C/year in post-monsoon, and +0.08 to +0.2 °C/year in winter. Overall, the study concludes that minimum temperatures are increasing faster than maximum temperatures, with temperatures rising across all stations and precipitation changes varying by season and location.

The study projects an increase in river discharge across all five gauge stations (Chisapani, Asaraghat, Benighat, Rimna, and Samaijighat) in the KRB of Nepal under both the SSP245 and SSP585 climate scenarios. This increase is expected in the near, mid, and far future periods. Annually, the projected increase in river discharge within the KRB ranges from 27.12% to 54.88% under the SSP245 scenario and from 45.4% to 93.3% under the SSP585 scenario, varying across different gauge stations. Seasonally, under SSP245, river discharge is predicted to change by +38.2 to +89.8% in pre-monsoon, +34.2 to +72.5% in monsoon, +28.8 to +32.9% in post-monsoon, and +9.2 to +16.4% in winter. Under the SSP585 scenario, the changes are even more pronounced: streamflow is expected to increase by +42.7 to +98% in pre-monsoon, +59.9 to +123.6% in monsoon, +35.8 to 62.8% in post-monsoon, and +17.4 to +35.3% in winter. Historically, the observed precipitation peak at all gauge stations has lagged the peak of river discharge by about a month. This delay is attributed to the natural melting of snow in the Himalayan region, which sustains river discharge even after a decrease in precipitation. However, with the projected increase in maximum and minimum temperatures, the snow is melting faster and more significantly. This change is expected to eliminate the historical delay, causing future peaks of precipitation and river discharge to coincide within the same month. The study demonstrates that the BFI varies significantly across seasons and climate scenarios, with the lowest values typically occurring in the monsoon season and the highest in the post-monsoon and winter seasons. This pattern indicates that baseflow contributes less to streamflow during wet seasons and more during dry seasons, providing crucial insights for water management in addressing regional droughts and floods.

The results of this study offer significant insights into how the watershed reacts to climate change, affecting the river basin’s water balance, despite multiple assumptions and limitations inherent in the models used. Evaluating various climate change scenarios has produced various outcomes, instrumental in facilitating effective planning and management of water resources in the KRB. Future research directions should also explore the implications of land use changes on future water flow dynamics.