Complementary use of multi-model climate ensemble and Bayesian model averaging for projecting river hydrology in the Himalaya

Considering the sensitivity and importance of water resources in the Himalayan uplands, this study intended to assess the hydrological responses to climate change in the Jhelum basin. Representative concentration pathway (RCP)–based projections from six dynamically downscaled global circulation models (GCMs) were bias-corrected for developing the climatic projections over the twenty-first century. The uncertainty associated with GCM outputs was addressed by using multi-model ensemble projections developed through Bayesian model averaging (BMA) technique. The assessment reveals that compared to the baseline (1980–2010) values, the annual mean maximum temperature in the basin will rise by 0.41–2.31 °C and 0.63–4.82 °C, and the mean minimum temperature will increase by 1.39–2.37 °C and 2.14–4.34 °C under RCP4.5 and RCP8.5, respectively. While precipitation is expected to decrease by 7.2–4.57% and 4.75–2.47% under RCP4.5 and RCP8.5, correspondingly. BMA ensemble projections were coupled with the Soil and Water Assessment Tool (SWAT) to simulate the future hydrological scenarios of the drainage basin. With the changing climate, the discharge of rivers in the Jhelum basin is expected to witness reductions by about 23–37% for RCP4.5 and 19–46% for RCP8.5. Moreover, the water yield of the basin may also exhibit decreases of 17–25% for RCP4.5 and 18–42% for RCP8.5. The projected scenarios are likely to cause water stress, affect the availability of water for diverse uses, and trigger transboundary water-sharing-related conflicts. The impact of climate change on discharge demands early attention for the formulation of mitigation and adaptive measures at the regional level and beyond.


Introduction
Changing climatic regimes and their impact on water resources have attracted considerable attention in the contemporary world because of their tremendous environmental and socioeconomic implications (Jasper et al. 2004;Luo et al. 2019;Reshmidevi et al. 2018). Mounting evidence exists about the impact of the increase in global mean temperatures on the regional water budget (Huntington 2006). Moreover, changing precipitation patterns and evapotranspiration rates projected under the climate change scenarios will further alter the hydrological systems (Kundzewicz et al. 2008;Raneesh & Santosh 2011). Therefore, while devising effective water management policies and adaptation strategies to offset climate change-induced stresses, it becomes quintessential to evaluate the hydrological responses and sensitivity of the riverine systems to the changing climate (Bhatta et al. 2019;Reshmidevi et al. 2018;Rodrigues et al. 2020).
Multiple studies have assessed the hydrological imprints of changing climate across the different river basins of the world, e.g., Liu et al. (2011) in the Yellow River basin, Tan et al. (2017) in the Kelantan River basin, Vetter et al. (2017) in 12 major sub-basins of the world, and Yu et al. (2018) in Yangtze River basin. Global circulation models (GCMs) are the popular tools in use for projecting climate change under different emission scenarios. However, GCMs provide information on a coarse scale with large grid sizes and insufficiently capture the regional heterogeneity that impacts the hydrological processes on the basin scale (Hakala et al. 2019). Regional climate models (RCMs) with a finer grid Responsible Editor: Marcus Schulz size have been developed to improve the outputs of GCMs. Although RCMs perform better in simulating the local climate, they inherit significant error/bias from GCMs (Chen et al. 2013). As a result, using the RCM outputs directly into the hydrological models can yield significant deviations in the model outputs with respect to the observed data (de Oliveira et al. 2017). As a result, prior to any hydrological simulations, bias adjustment of the raw RCM results is critical (Luo et al. 2019). "Coupled Model Intercomparison Project Phase 5(CMIP5)" GCMs provide the outputs based upon the representative concentration pathways (RCPs) that describe the net radiative forcing levels until 2100 ranging from a mitigating pathway (RCP2.6), medium stabilizing pathways (RCP4.5/RCP6), and extreme climatic pathway (RCP8.5) (Meinshausen et al. 2011;Taylor et al. 2012).
GCM output carries a substantial amount of uncertainty with it and must be taken into cognizance while carrying out any impact assessment studies. The cascade of uncertainty is nestled in the range of socioeconomic or developmental pathways that the world may take in the future; incomplete knowledge or representation of the climatic system and its interactions; and the structural or parametric differences in the GCMs. As such, the single model's results must be interpreted cautiously (Reshmidevi et al. 2018). To this end, ensemble projections have been proposed on account of their reliability and uncertainty assessment (Khan et al. 2021). In comparison to the outputs from separate models, the multimodel ensemble has been demonstrated to be more effective in the representation and simulation of climatic variables (Gleckler et al. 2008;Knutti et al. 2010;Zhang & Huang 2013). The ensemble techniques can vary from the simple arithmetic ensemble mean wherein each model is weighted equally irrespective of its performance, to the weighted ensemble mean. Additionally, ensemble projections have been generated through the adoption of the median values from the projections of multiple climate models (Cloke et al. 2013;Yang et al. 2022). These techniques improve the simulation results; however, uncertainty in the estimates is not addressed. BMA is a weighted ensemble method in which the relative competence of models is determined and model weights are optimized so that BMA ensemble projections closely match the observed data (Massoud et al. 2020). In this way, higher weights are applied to more skillful models as compared to low-skillful models. Additionally, the use of BMA is justified by the fact that uncertainty is included in the estimate in a substantially more efficient manner than when using a simple mean or a median ensemble. BMA has been previously used for developing ensemble projections in studies like Huang, (2014), Khan et al., (2021), and Massoud et al., (2020).
Generally, hydrological models forced with climatic scenarios have been employed to simulate the possible implications on stream flows, catchment storage, and other water balance components. Different hydrological models that can replicate the hydrological properties and processes operating within a basin have already been developed . Soil and Water Assessment Tool (SWAT) is being widely used for evaluating the likely implications of changing climatic regimes on hydrological processes (Bajracharya et al. 2018;Bhatta et al. 2019;Rodrigues et al. 2020;Narsimlu et al. 2013;Touseef et al. 2021). Numerous studies have also employed the SWAT model for simulating the combined impact of changing land-use patterns and climatic scenarios with reasonable precision (Ahmed et al. 2022;Yonaba et al. 2021a). Furthermore, the SWAT model contains a snowmelt module relevant for modeling the basins where snowmelt is a significant contributor to the streamflow. In addition to this, the SWAT model has been successfully used for simulating the streamflow in and around the Jhelum basin, e.g., Saddique et al. (2019), Shah et al. (2020, and Ougahi et al. (2022). Hence, the applicability of the SWAT model was extended to the present study.
The potential implications of climate change are more noticeable in the mountainous regions like the Himalayas, particularly on water resources, where snowmelt and glacial ice melt are the major determinants of streamflow (Viviroli et al. 2007;Lutz et al. 2014;Immerzeel et al. 2013). Significant increases in the temperature along with declining precipitation patterns have been reported for the Jhelum basin (Ahsan et al. 2021(Ahsan et al. , 2022. Having established climate change signals, it becomes obvious that the water resources will be most affected in this region. The imprints of the observed climatic variability are already visible over the hydrological regime of the Jhelum basin. For example, climate change is putting a lot of strain on the aquatic systems in the region (Alam et al. 2020). Lone et al. (2022) reported significant decreasing trends in the streamflow, driven by observed climate variability. Hence, climatic changes are likely to further impact the water security of the region and challenge the existing management strategies. The hydrological consequences of the climatic changes can transmit to the overall prosperity of the region due to dependence of other sectors, viz., irrigation, domestic, industrial, tourism, and hydropower generation on the water resource base. Limited studies have been attempted so far in the Jhelum basin that project the hydrological responses to changing climate, e.g., Singh Jasrotia et al. (2021). However, the key limitations of this study are that it uses outputs of a single RCM and lacks any detailed assessment of water balance components. Hence, it is imperative to have a rigorous quantification of the hydrological responses to the projected changes in climatic regimes. For this purpose, the current study utilizes outputs from 6 GCMs under a medium stabilizing pathway (RCP4.5) and an extreme pathway (RCP8.5). The study primarily aims to develop robust climatic projections for the twenty-first century, using the multi-model ensemble based on the Bayesian model averaging framework. The study subsequently employs the hydrological model SWAT coupled with the climatic scenarios to project the changes in the streamflow and different components of water balance over the twenty-first century. No such comprehensive study has been attempted yet for the Jhelum basin; the present one can thus serve as a benchmark study while planning for the mitigation of hydrological implications ascribed to climate change and pave the way for the formulation of adaptive policy for sustainable water management.

Study area
The present study is attempted for the Jhelum basin ( Fig. 1) located in the lap of the Himalayas with an area of about 15,000 km 2 . The trunk stream, Jhelum constitutes an important tributary of the Indus River system forming an elongated bowl-shaped basin that originated from the collision between the Indian and Eurasian plates (Alam et al. 2015(Alam et al. , 2017. Based on stratigraphy and elevation, the basin can be divided into 3 major physiographic divisions, viz., mountainous uplands, extensive plateaus of lacustrine deposits (Karewas), and valley floor with the numerous streams draining into Jhelum (Bhat et al. 2018). The slope and elevation of the basin lie between 0 to 65° and 1450 to 5500 mamsl, causing significant microclimatic variations. The region's climatic regime is controlled by two distinct weather patterns, viz., Western disturbance in winter months and southwest monsoons in the summer half. The impact of the latter is somewhat limited, and the major proportion of the precipitation is received in the winter months, mainly in the form of snowfall. The winter precipitation feeds the glaciers and acts as a buffer in maintaining the streamflow of the river Jhelum during the dry periods. Snowmelt currently Fig. 1 Location of the Jhelum basin with the major streams, climatic stations, and the river gauges contributes about 29% during the spring season, 58% during the summer season, and 38% during the autumn season to the streamflow of river Jhelum (Lone et al. 2022). The annual mean precipitation in the valley is around 100.5 cm with about 71.6% received during Oct-May and 28.4% during the Jun-Sep (Dad et al. 2021). The climate of the region is moderate having a sharp seasonality and has been categorized as a Sub-Mediterranean type (Bagnolus and Meher-Homji 1959). The valley records an annual mean maximum temperature of 17.6 ± 0.8 °C and a minimum temperature of 5.4 ± 0.4 °C (Dad et al. 2021). However, the annual temperature variations of the region range from sub-freezing during the winter and up to as high as 35 °C during the summer. The two major types of soil found in the basin are eutric cambisols on the valley floor and lithosols in the peripheral uplands (Romshoo et al. 2020). The region supports a population of ~ 7 million persons with heavy reliance on water resources for sustenance.

Observed hydroclimatic data
Observed climate data  for temperature and precipitation variables at a daily scale was procured from the India Meteorological Department, Regional Office Srinagar, J&K, India. This data was utilized for the bias correction of the GCM data and also as input for the baseline run of the SWAT model. For the other SWAT input climatic variables, i.e., relative humidity, wind speed, and solar radiation, and "Climate Forecast System Reanalysis (CFSR)" data were used. CFSR data have been used and validated in the past for the Kashmir basin by Ahsan et al. (2021). The observed streamflow data for the 3 gauges on the trunk stream, viz., Sangam, Asham, and Ram Munshi Bagh, were used in calibrating and validating the SWAT model.

Climate scenario projections
The Coordinated Regional Climate Downscaling Experiment (CORDEX) was developed under the aegis of the World Climate Research Programme to stimulate climatic scenarios and improve regional impact assessment. COR-DEX over the south Asian domain provides data for multiple GCM-RCM combinations at a grid size of 0.44° (~ 50 km). CORDEX-SA has previous successful applications in the Himalayan regions (e.g., Dimri et al. 2018;Krishnan et al. 2019). A total of 6 GCMs (Table 1) downscaled dynamically using the RegCM4 RCM (Giorgi et al. 2012) were used in the projection of the climatic scenarios under RCP4.5 and RCP8.5. The projected span (2006-2099) was split into 3 periods, viz., 2011-2040 (2020s), 2041-2070 (2050s), and 2071-2099 (2080s). The relative changes in mean climatology over the different spans of the twenty-first century were computed from the baseline  climatology.

Bias correction of GCM/RCM data
Prior to using the GCM/RCM data for the projection of climatic changes or their hydrologic implications, it is necessary to apply the bias correction techniques because the data obtained from the climate models inherit systematic biases (Luo et al. 2019). A wide range of bias correction techniques have been devised (Teutschbein and Seibert 2012); and in the present study, variance scaling (VS) was used for biascorrecting the temperature data, while a hybrid of power transformation (PT) and local intensity scaling (LOCI) techniques were employed for bias correcting the precipitation data. These bias correction methods use monthly correction factors to adjust the daily climate data.
Variance scaling is an effective technique and has the advantage of correcting both the mean and variance (Teutschbein and Seibert 2012). VS uses Eq. 1 to correct the bias in temperature. Power transformation of the precipitation data also corrects both mean and variance in the raw data. However, the limitation of the PT method is that it does not correct the wet day frequency. Hence, before its use, the wet day frequency is adjusted using the local intensity scaling method. In LOCI (Eq. 2) initially, a threshold for wet day P thresh,m is obtained from the raw data to match the wet day frequency in raw and bias-corrected precipitation. This is followed by the calculation of a scaling factor S m = (Pobs,m,d|Pobs,m,d>0) (Praw ,m,d |P raw ,m,d >P thres ,m) to match the mean of raw and bias-corrected precipitation.
After adjusting the wet day frequency using the LOCI method, in the PT method, b m is initially estimated to mini- is computed to match the mean of LOCI corrected and observed precipitation. In the final step, Eq. 3 is used to correct the bias.
where T is the temperature, P is precipitation, is the mean, is the standard deviation, e.g., T bc ,m,d is the bias-corrected temperature for the dth day of mth month; P obs,m is the standard deviation of the observed precipitation series for mth month (Fang et al. 2015).
To check the skill and efficiency of the RCM simulations, RCM data (raw and bias-corrected) were compared with observational data on a monthly timescale using the efficiency indicators like coefficient of determination (R 2 ), percentage bias (PBIAS), Nash Sutcliffe efficiency (NSE), and Kling-Gupta efficiency (KGE). Moreover, Taylor diagrams (Taylor 2001) were employed to assess the competence of RCMs for both raw and bias-corrected data. R 2 given by Eq. 4 is a measure of goodness-of-fit among the observed and model data ranging between 0 and 1. R. 2 values nearing 1 imply a higher agreement among observed and model values. PBIAS (Eq. 5) measures the model data's mean tendency to be greater or lesser than actual observations. PBIAS has an optimum value of 0 with positive values depicting underestimation and negative values depicting overestimation in model simulations. NSE is a standardized statistic (Eq. 6) that measures deviations between the modeled and observed data normalized by the observational data variance (Takele et al. 2021). It values between -∞ to 1 and NSE = 1 depicts modeled data are equal to that of the observed data (Nash and Sutcliffe 1970). KGE (Eq. 7) is an integrated efficiency criterion evaluating correlation, bias, and variability of model output in comparison to the observed data (Kling al. 2012). Varying between -∞ and 1, the model performance can be deemed as poor (0 ≤ KGE ≤ 0.5), acceptable (0.50 ≤ KGE ≤ 0.75), good (0.75 ≤ KGE ≤ 0.9), or excellent (KGE ≥ 0.9) (Yonaba et al. 2021a) where O i are the observational data values and M i are the model simulated values.
where r is the correlation coefficient between observed and simulated values, β is the ratio of simulated and observed means, and γ is the ratio of simulated and observed coefficients of variation.

Bayesian model averaging (BMA) for multi-model ensemble
On account of parametric and structural differences, GCMs/ RCMs produce different outputs (Wallach et al. 2016). Hence, a single climate model's output may still be subject to uncertainty (Uusitalo et al. 2015). The bias-corrected RCM output with different GCM forcing was utilized to develop multi-model ensemble projections employing the Bayesian model averaging (BMA) technique. BMA is an advanced statistical technique for deriving the multi-model ensemble and combines the individual models using optimized weights; weights being proportional to the relative skills of models in reproducing observed data over a historical baseline. The regression model in BMA takes into account all conceivable combinations of variables (in this case, GCMs) before calculating a weighted average. The posterior model probability (PMP) used as a weight is a measure of the model's skill during training (Raftery et al. 2005). Let k be the number of available models; the maximum count of regression models is s = 2 k . Here in the study using 6 models, the number of regression models is s = 2 6 = 64. PMP of n-th model M n | D , in 2 k candidate regression models is computed using Bayes theorem as given in Eq. 8 where p D | M n is the observed data likelihood given n-th model, p M n is the n-th regression model's prior probability, and p(D) is a constant used for normalizing as given in Eq. 9: Different forms of priors are available in the literature; for this investigation, we selected a uniform prior distribution based on Khan et al., (2021), which gives equal weight to all possible regression models. For the likelihood function p D | M n , Gaussian likelihood (Vogel et al. 2008) was used following Khan et al. (2021). Equatin 10 gives the variable's conditional prediction PDF based on training data.
where p y | M n is the prediction probability density function (PDF) based on the n-th regression model and p M n |D is the associated posterior probability and is being used as weight (Khan et al. 2021). The BMA analysis and visualization were carried out with R software using the various packages viz., BMA, BMS, BAS, and ensembleBMA.

SWAT model setup
Hydrological modeling is a key aspect in evaluating the hydrological implications of climatic changes (Praskievicz and Chang 2009). The basin was modeled by employing the SWAT hydrological model to project the changes in streamflow and water budget components for the twenty-first century. The water balance in the model is governed by Eq. 11: where SW t denotes the water content of soil at t timestep, SW 0 denotes the initial water content of the soil, R day is daily precipitation, Q surf denotes the surface runoff, E a denotes the actual evapotranspiration, w seep denotes the percolation, and Q gw denotes the groundwater flow (Bajracharya et al. 2018). The units of all components are in millimeters (mm).
Snowmelt was incorporated in the runoff calculations and is governed by Eq. 12: where SNOWmlt denotes the amount of daily snowmelt (mm), bmlt denotes the factor of daily melt (mm/day °C), snowcov denotes the fraction of snow-cover in the hydrological response unit (HRU) area, Tmax is the day's maximum temperature, Tsnow is the temperature of the snowpack (°C), and Tmlt denotes the snowmelt threshold temperature (°C) (Bhatta et al. 2019). SWAT uses various inputs, viz., climate data, soil map/ attributes, land use map/attributes, digital elevation model (DEM), and the observed streamflow. The climatic input variables used are given in "Observed hydroclimatic data". For the current study, the ALOS PALSAR Digital Elevation Model having a grid size of 12.5 m was used ( https:// search. asf. alaska. edu/#/). The land use/land cover data for the region was produced by the supervised classification of Landsat 8 OLI imagery acquired on 2013/06/03, using a maximum likelihood classification algorithm in ArcGIS 10.2. A total of 11 classes, viz., agriculture, built-up, barren, marshes, pasture, sparse forest, horticulture, water, scrub, dense forest, and glacier/ice were identified. The accuracy of the resultant map was assessed using the Kappa coefficient and had a value of 0.87. Values greater than 0.80 indicate robustness to the perfect amount of agreement (Yonaba et al. 2021b). For soil data, a 30-arc-second raster dataset from the Harmonized World Soil Database (HWSD) of the Food and Agriculture Organization (FAO) was employed.
DEM was employed as an input in the QSWAT interface in the QGIS environment to divide the catchment into subbasins. A threshold of 10,000 ha was used in the delineation of the drainage network and the analysis identified 41 sub-basins in the Jhelum basin. The smallest model unit is the hydrological response unit (HRU), which is defined by a unique fusion of soil type, slope class, and land use (Romagnoli et al. 2017). Since the seasonal snowmelt determines the basin's hydrology, elevation bands were created to represent the snowmelt as well as the orographic controls of precipitation and temperature. SWAT model provides a maximum of 10 elevation bands for a particular sub-basin and in the present analysis 5 elevation bands were set for the sub-basins having a significant attitudinal range whereas for sub-basins with less orographic differences only one elevation band was used, e.g., Bajracharya et al. (2018).

SWAT model calibration and validation
Multi-site calibration and validation were carried out for the 3 gauges on the trunk stream, viz., Sangam, Asham, and Ram Munshi Bagh. It has been demonstrated that multisite calibration considerably enhances the spatial patterns as depicted by the SWAT model (Gbohoui et al. 2021). Observed monthly streamflow data for 10 years (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) were used in calibration while as 5 years (2012-2016) were used in the validation of the model. Before the calibration of the model, sensitivity analysis was carried out to filter the most influential parameters determining the streamflow variations. The parameters were ranked as per their sensitivity towards the streamflow using the global sensitivity analysis embedded in the SUFI-2. It employs the p-stat and t-stat to determine the most sensitive parameters, having lower p values and higher t values (Bhatta et al. 2019). The SWAT-CUP considers the parameter uncertainty, model uncertainty, input error, and yields uncertainty measures like 95% prediction uncertainty (95PPU), r-factor, and p-factor. The r-factor (Eq. 13), which runs from 0 to ∞, is the mean width of the 95PPU divided by the standard deviation of the observational datasets. Similarly, the p-factor refers to the percentage of the observed data series that is encompassed by the 95PPU and ranges between 0 and 1. Satisfactory model performance is indicated by p-factors greater than 0.7 and r-factors less than 1.5 (Abbaspour et al. 2015).
where x t i ,97.5% s is the upper limit of the 95PPU band and x t i ,2.5% s is the lower limit of the 95PPU band at t timestep in the ith simulation. 0 is the observed standard deviation and n is the number of data points.
Statistical indicators such as R 2 , PBIAS, and NSE were also utilized to assess model performance. NSE was also used as an objective function for parameter optimization during the calibration of the model. Moriasi et al. (2007) and Almeida et al. (2018) have defined the thresholds of satisfactory and acceptable model performance for 0.50 < R 2 ≤ 0.60, ± 15 < PBIAS ≤ ± 25, and 0.36 < NSE ≤ 0.60.
Multi-model ensemble climatic projections were coupled with the calibrated model parameters to simulate the hydrological implications of changing climatic patterns in the Jhelum basin. After ascertaining the robustness of BMA projections, these were utilized as an input to project the streamflow at Ram Munshi Bagh, Asham, and Sangam gauges under RCP4.5 and RCP8.5. Moreover, the impact of climatic changes was assessed separately on the different components of catchment water balance, viz., evapotranspiration, snowmelt, surface runoff, and water yield. For estimating evapotranspiration (ET), the Hargreaves method (Eq. 14) was used (Hargreaves and Samani 1985). Due to the paucity of observed data, the SWAT simulated output for these components during the baseline run (2001-2010) was used as a reference to estimate the relative changes over the twenty-first century (Bajracharya et al. 2018).
where Tm, Tx, and Tn are the mean, maximum, and minimum air temperatures in °C, respectively. R a is the water equivalent of solar radiation in mm/day.

Bias correction of the GCM/RCM output
The competence of the bias correction techniques was evaluated by the comparison of RCM data (raw and BC) with the observational data on a monthly scale, using efficiency metrics, viz., R 2 , PBIAS, and NSE (Table 2). For Tmax and Tmin, R 2 values were high before bias correction but PBIAS and NSE got significantly improved. The raw Tmax (Tmin) displayed substantial biases, and bias correction minimized the PBIAS from − 60.6 to − 37.40 (− 135 to − 47.2) to 1.9-3.9 (2.3-6.5). Likewise, the NSE values for Tmax (Tmin) in raw RCM were of the order − 0.84 to In contrast to this, GCMs tend to resolve the precipitation poorly as indicated by lower values of efficiency metrics. However, the bias correction could improve the match between the simulated and the observational data (Table 2). Furthermore, Taylor diagrams (Taylor 2001) are used to display the RCM's performance (raw and BC) (Fig. 2). Considering the uncertainty in the single RCM output, multi-model ensemble projections were derived using BMA and were also compared with the baseline data. The BMA ensemble data matched more closely with the observed data than any of the individual RCMs, supported by the increased values of the efficiency indicators ( Table 2). The PBIAS in the BMA estimates of climatic variables reduced to 0 along with further improvements in R 2 , NSE, and KGE. Figure 3 depicts the performance of the BMA technique during the training period.

Climate change projections for the twenty-first century
The projections from the individual models and the BMA ensemble for the Jhelum basin are presented in Table 3. The results infer substantial warming in the annual mean maximum (Tmax) and minimum (Tmin) temperatures with varying rates of warming among the different GCMs. Multimodel ensemble projections show that under RCP4.5, the Tmax (Tmin) of the region will increase by 0.41 (1.39) °C during 2020s, 1.80 (2.12) °C during 2050s, and 2.31 (2.37) °C during 2080s. The rate of the warming for Tmax (Tmin) enhances under the RCP8.5 forcing levels by increments of 0.63 (2.14) °C during 2020s, 2.56 (3.12) °C during 2050s, and 4.82 (4.35) °C during 2080s. While comparing the magnitude of change simulated by the different models, it varied substantially; e.g., for Tmax, it varied from 0.81 °C (CNRM-CM5) to 1.42 °C (MPI-ESM-MR) for the 2020s, 1.59 °C (CNRM-CM5) to 3.01 °C (CSIRO-Mk3) for 2050s, and 2.17 °C (CNRM-CM5) to 3.66 °C (CSIRO-Mk3) for 2080s under RCP4.5 forcing levels. CNRM-CM5 showed consistently the least amount of change for Tmax and Tmin during the different spans of the twenty-first century under both RCPs. The station-scale BMA projections (Table 4) reveal higher warming rates in Tmax and Tmin for the Gulmarg, while the lowest warming rates were found for Srinagar (Tmax) and Pahalgam (Tmin). Moreover, ensemble estimates project enhanced night-time warming in the region over the twenty-first century pinpointed by the higher warming rates of Tmin as compared to the Tmax. The asymmetrical warming patterns of Tmax and Tmin exist even between models and also among the different spans of the twenty-first century, e.g., the CSIRO-Mk3 model projects higher rates of warming in Tmin than Tmax for the 2020s and vice versa during the 2050s and 2080s. The BMA projections for the Tmax and Tmin in the Jhelum basin are shown graphically in Fig. 4a and Fig. 4b, respectively.
Precipitation projections for the Jhelum basin show reductions over the twenty-first century vis-a-vis the baseline . The BMA estimates project decreases by 7.2-4.57% (RCP4.5) and 4.75-2.47% (RCP8.5) during the different periods of the twenty-first century. Even though precipitation tends to increase from the 2020s to the 2080s (Fig. 4c), there is a net decrease when compared to baseline values. The diminishing patterns of precipitation rhyme across all the individual models although with higher magnitudes than BMA estimates (Table 3), except IPSL-CM5A-LR which projects an increase of 9.6% in the precipitation of the Jhelum basin during the 2080s under RCP8.5. The station-scale BMA estimates (Table 4) show enhanced decreases for Kokernag (RCP4.5 and 8.5), Kupwara (RCP4.5 and 8.5), and Gulmarg (RCP8.5). In contrast to this, stations like Gulmarg (RCP4.5) and Qazigund (RCP 8.5) witness increases in the precipitation, while Srinagar shows almost no change (< 1%) over the twenty-first century when compared with the baseline precipitation.

SWAT model calibration and validation
Any hydrological model must be adequately calibrated and validated before adopting its outputs for further analysis. The SUFI-2 algorithm, embedded within the SWAT-CUP, was employed to calibrate and validate SWAT. Multi-site calibration was done for Sangam, Asham, and Ram Munshi Bagh gauges, located on the trunk stream using the observed mean monthly streamflow values for 2000-2010. The model was validated using 5 years (2011-2016) of streamflow data. Prior to that, sensitivity analysis was executed to limit the count of parameters. Table 5 shows the 11 most sensitive parameters as well as their p value and t-stat, important in the estimation of streamflow variability. The model conformity during its calibrating and validating process was assessed using 3 statistical indicators, viz., R 2 , PBIAS, and NSE. Furthermore, the uncertainty of the model outputs was evaluated using the r-factor and p-factor. The r-factor varied between 0.69 and 0.84 for the calibration period and 0.63 and 0.74 for the validation period. Similarly, the p-factor was found to be 0.39-0.47 and 0.43-0.54 for the calibration and validation phases. Overall, the uncertainty values reflect the satisfactory model performance; however, the r-factor in some instances was slightly less than the threshold value (0.7). Given the dearth of quality observed streamflow data for the region, such values can be taken as the acceptable limits of model performance. Table 6 shows the performance metrics of the SWAT model and the values reflect that model performance was well within the acceptable limits. The close agreement between actual and modeled streamflow (Fig. 5)   Fig. 2 Performance of GCMs and BMA ensemble during training period ; a mean maximum temperature, b mean minimum temperature, and c precipitation. BC = biascorrected GCM output, R = raw GCM output Environmental Science and Pollution Research (2023) 30:38898-38920 reveals that the validated model could be employed with considerable confidence and accuracy to model hydrologic responses to climatic changes.

Projection of streamflow under climate change scenarios
After the calibration and validation, SWAT was employed to translate the climatic change signals into streamflow responses using the multi-model ensemble climatic projections obtained using BMA. The relative changes in the streamflow during the different spans of the twenty-first century were compared with the baseline (2001-2010) values for 3 gauges (Sangam, Asham, and Ram Munshi Bagh) located along the course of the trunk stream (Table 7). In the wake of rising temperatures and declining precipitation patterns, the Jhelum basin is going to witness substantial reductions in streamflow over the twenty-first century.

Projection of water budget components under climate change scenarios
The variations in the streamflow are governed and linked with the different components of water balance like precipitation, evapotranspiration, surface runoff, snowmelt, etc. After having the precipitation projections for the twentyfirst century, the current study also intended to assess the relative changes in other water balance components, viz., actual evapotranspiration (ET), potential evapotranspiration (PET), snowmelt, surface runoff (SurQ), and water yield (Wyield) under climate change scenarios (Table 8). These components were simulated using the SWAT model for the baseline period (2001-2010) from the observational climate data and BMA ensemble climatic projections were employed for future simulations. Due to a lack of observed data, the SWAT simulated output for these components during the baseline run (2001-2010) was used as reference data (Bajracharya et al. 2018). The percentage changes in the water budget components throughout various periods in the tnwety-first century, under both RCPs, were calculated by comparing them to baseline values. Figure 7 shows the projected changes in the Jhelum basin's average annual hydrological components over different periods in the future. ET was estimated from the Hargreaves method (Hargreaves and Samani 1985), and it witnessed increases over the twenty-first century in the Jhelum basin. Actual ET shows increases of order 18% (20%) during 2020s, 28% (30%) during 2050s, and 31% (40%) during 2080s under RCP4.5 (RCP8.5). The changes in the ET for different subbasins over the twenty-first century under both RCPs are shown in Fig. 8. ET projections show a progressive increase for the future, determined by levels of radiative forcing and time, viz., near, mid, or end century. This increase in the ET is consistent with the projected increase in the temperatures driven by the enhanced radiative forcing levels. Due to the continuously rising temperatures, the ET values show a consistent rise over the twenty-first century. ET constitutes a major abstraction in the water balance which reduces the net water availability for the runoff and is a measure of the irrigation water requirements.
Snowmelt contribution in the Jhelum basin witnessed significant decreases in the future under the changing climatic regimes. The snowmelt shows substantial reductions by about 53% (50%) during 2020s, 58% (62%) during 2050s, and 60% (72%) during 2080s under RCP4.5 (RCP8.5). Moreover, Fig. 9 reveals that most of the decreases will be observed in the sub-basins situated in the lower elevation. The decrease in snowmelt contribution may be attributed to the declining precipitation trends over the twenty-first century. Jhelum basin is projected to witness decreases in the overall precipitation, and this would mean a reduction in the winter season precipitation too that is in the form of snowfall. Significant decreasing trends are already ominous in the winter season precipitation of the Jhelum basin (Ahsan et al. 2021). Furthermore, the warming temperatures would also reduce the precipitation amount falling as snow and there will be more liquid precipitation. Water yield is a measure of the cumulative water delivered to streamflow by sub-basins/HRUs and includes baseflow, lateral flow, and surface runoff (Andrade et al. 2021;Jain et al. 2017). The average water yield within the Jhelum basin witnessed a decline over the twenty-first century. For the RCP4.5 (RCP8.5), the water yield witnessed reductions by 17.7% (18.4%) during 2020s, 24% (32%) during 2050s, and 25% (42%) during 2080s. Figure 10 shows the spatial variations in the projected water yield of different sub-basins for different spans of the future. The projected decreases in the water yield can be attributed to the enhancement of the ET levels and decreases in the snowmelt contribution over the twenty-first century. The progressive decline in the water yield during the twenty-first century infers that water availability will get reduced significantly and will aggravate the water scarcity and stresses in the basin.

Discussion
The water resources of the mountainous basins are highly exposed to climatic changes wherein the regional hydrology is dependent on the glacial ice melt and seasonal snowmelt characteristics. The present study is attempted in the Jhelum basin-nestled in the complex folds of the North-Western Himalayas and forms important headwaters of the Indus River system. The imprints of global climate change are well established from the observational datasets in the Jhelum basin and its impacts are already being felt (Ahsan et al. 2021). Few studies (e.g., Ahsan et al. 2022) have also attempted to estimate the magnitude of climatic changes over the twenty-first century. However, these studies are limited in scope because the projections are based on a single GCM. It has been found that GCM outputs inherit considerable uncertainty in them and thus use of a single GCM can yield inaccurate results (Khan et al. 2021). For modeling the hydrological consequences of climatic changes, the development of rigorous and robust climate change projections is a primary requirement. To this end, 6 GCMs dynamically downscaled by the RegCM4 RCM from the CORDEX-SA experiments under a medium stabilizing forcing level (RCP4.5) and an extreme forcing level (RCP8.5) were collated for the present study. The use of RCM output is preferred for the basin-scale assessment than the GCM due to their higher resolution and the local sub-grid processes are parametrized better (Seager & Vecchi 2010). However, the raw RCM outputs may still contain significant biases, particularly over mountainous areas (Fowler et al. 2007). As a result, bias-correcting the RCM outputs is critical, and in this work, the temperature was bias-corrected using the VS approach, and precipitation was bias-corrected using a mix of LOCI and PT techniques. The bias correction significantly improved the performance indicators and yielded an overall better match between the observational data and the RCM simulation results during the baseline . Furthermore, to narrow the uncertainty, multi-model ensemble climate change projections were developed using the Bayesian model averaging (BMA) technique. There is consensus among the climate modeling community that by using the multi-model ensemble climatic projections, the uncertainty can significantly reduce, given the fact that individual models may simulate different aspects of the climate system well and the errors get canceled or reduced in the process (Brient 2019). Moreover, it has been reported that by using multiple simulations from a given RCM for different GCM forcings and emission scenarios, the accuracy of climatic projections can be further improved (Räisänen et al. 2004;van den Hurk et al. 2005;Salathé 2005). The climate change projections from all the GCMs and BMA estimates show unequivocal warming for the Jhelum basin over the twenty-first century. The temperature increments are progressive over time and vary directly with levels of radiative forcing levels. BMA estimates project warming of Tmax by 0.41-2.37 (0.63-4.82) °C during the various spans of the future under RCP4.5 (RCP8.5). Tmin is getting warmer at a relatively higher rate vis-a-vis Tmax and increases are of the order 1.39-2.37 °C under RCP4.5 and 2.14-4.35 °C under RCP8.5. Precipitation in the Jhelum Table 7 Projected changes (%) in streamflow of Jhelum under multi-model ensemble climatic projections for the twenty-first century basin witnessed reductions by varying amounts across different GCMs. These changing climatic regimes are going to trigger several implications for key sectors, viz., water resources, agriculture, energy consumption, and human health (Ahsan et al. 2021). The possible implications on the streamflow and water budget of the basin were assessed by coupling the climate change scenarios into the SWAT model. SWAT model performed quite well in simulating the basin hydrology as inferred by high values of the model skill metrics. The streamflow of the river Jhelum is projected to witness substantial reductions over the twenty-first century under both RCPs. The decreases aggravate over time and forcing levels with the highest decrease of about 40-46% for the 2080s under RCP8.5. These reductions in streamflow   Bajracharya et al. (2018), and Reshmidevi et al. (2018). Snowmelt is another important component of regional hydrology and sustains the streamflow in the Jhelum basin during the dry months (Jeelani et al. 2012). Under the changing climate, snowmelt witnessed a reduction of 50-70% over the twenty-first century. The decrease in the snowmelt may be attributed to the coupled impact of warming temperatures and diminishing precipitation. Higher temperatures particularly during the winter would cause more precipitation in the form of rain than snow and would ultimately decrease the snowmelt contributions. The decreasing snow and increasing rain proportion in the precipitation are already being observed in the basin as reported by Romshoo et al. (2018). This will also lead to the enhanced ablation of the glacial reserves in the basin. The glacial recession has amplified significantly in  the region during the last few decades (Murtaza & Romshoo 2016). Another important aspect of the snowmelt characteristic driven by the increasing temperatures is the earlier melting of the snowpacks. This impact is discernible during the 2080s under RCP8.5 where the peak streamflow month is showing a shift from May to March (Fig. 6c). This shifting of peak streamflow has been reported in numerous other snow-fed basins of the world and tends to create issues for water resource management due to the shifting of peak streamflow period away from peak demand season (Barnett et al. 2005). The combined impacts of the ET and snowmelt are yielding a decrease in the water yield of the basin and hence streamflow. The current study evaluated the effects of climate change on the streamflow and water balance components using static land use. However, the regional water balance could potentially vary as a result of changes in land use in both the present and future states. Over the recent past, there have been rampant land use changes in the basin with a considerable reduction in water bodies and a progressive increase in the built-up area and horticultural area (Alam et al. 2020). The effects of climate change on streamflow may get exacerbated or get better as a result of these changes. The integration of dynamic land use in the SWAT model can further narrow the band of uncertainty and make the inferences more robust (Yonaba et al. 2021a). Hence, this still stands as the study's primary limitation, and future research should attempt to model the combined effects of land use and climate change scenarios. Also, under the scenario of diminishing surface water availability, groundwater would be a key resource in sustaining the future water demand (Sheikh et al. 2022). There is also a significant requirement for modeling groundwater availability vis-à-vis climate change for integrated water resource management.

Conclusions
The discharge from the rivers in the Jhelum basin is of tremendous importance from the perspective of ecology, environmental health, and water security. The river waters are generally tapped for irrigation, hydropower generation, domestic consumption, and recreation purposes. Hence, any changes in hydrological conditions induced by climatic changes will entail a profound influence on the environment and society. With the application of multiple downscaled GCMs and BMA, the present study formulated robust climate change projections with due account of the inherited GCM uncertainties. The hydrological repercussions of climate change in the Jhelum basin for the twenty-first century were assessed with the SWAT model. The changes are projected for the annual mean maximum and minimum temperatures, precipitation percentages, and river discharge. Moreover, changes in the catchment water balance were simulated under changing climatic regimes. This analysis found that the maximum temperature may rise to 4.82 °C and the minimum temperature may go up by 4.35 °C coupled with a reduction in the precipitation to the maximum of 7.2% in the given period. Consequent to the changes in temperature and precipitation conditions, the discharge of the rivers is likely to fall by 17 to 42%. The significant drivers of diminishing streamflow are expected to be an increase in the evapotranspiration levels and a decrease in the snowmelt contributions. Temperature enhancements, the decline in precipitation, and the huge decrease in the seasonal flow of rivers may result in dwindling aquatic ecosystems, loss of biodiversity, reduced productivity, and conflicts in watersharing agreements, e.g., Indus Water Treaty. The deliverables of the study may therefore serve as a guide for foreseeing future challenges and developing policies to mitigate and adapt to the changing climate and the flow regimes of the river basins in the Himalaya.