Machine-learning- and deep-learning-based streamflow prediction in a hilly catchment for future scenarios using CMIP6 GCM data

Abstract

The alteration in river flow patterns, particularly those that originate in the Himalaya, has been caused by the increased temperature and rainfall variability brought on by climate change. Due to the impending intensification of extreme climate events, as predicted by the Intergovernmental Panel on Climate Change (IPCC) in its Sixth Assessment Report, it is more essential than ever to predict changes in streamflow for future periods. Despite the fact that some research has utilised machine-learning- and deep-learning-based models to predict streamflow patterns in response to climate change, very few studies have been undertaken for a mountainous catchment, with the number of studies for the western Himalaya being minimal. This study investigates the capability of five different machine learning (ML) models and one deep learning (DL) model, namely the Gaussian linear regression model (GLM), Gaussian generalised additive model (GAM), multivariate adaptive regression splines (MARSs), artificial neural network (ANN), random forest (RF), and 1D convolutional neural network (1D-CNN), in streamflow prediction over the Sutlej River basin in the western Himalaya during the periods 2041–2070 (2050s) and 2071–2100 (2080s). Bias-corrected data downscaled at a grid resolution of 0.25${}^{\circ }$ $\times$ 0.25${}^{\circ }$ from six general circulation models (GCMs) of the Coupled Model Intercomparison Project Phase 6 GCM framework under two greenhouse gas (GHG) trajectories (SSP245 and SSP585) were used for this purpose. Four different rainfall scenarios ($R_{\mathrm{0}}$, $R_{\mathrm{1}}$, $R_{\mathrm{2}}$, and $R_{\mathrm{3}}$) were applied to the models trained with daily data (1979–2009) at Kasol (the outlet of the basin) in order to better understand how catchment size and the geo-hydromorphological aspects of the basin affect runoff. The predictive power of each model was assessed using six statistical measures, i.e. the coefficient of determination ($R^{\mathrm{2}}$), the ratio of the root mean square error to the standard deviation of the measured data (RSR), the mean absolute error (MAE), the Kling–Gupta efficiency (KGE), the Nash–Sutcliffe efficiency (NSE), and the percent bias (PBIAS). The RF model with rainfall scenario $R_{\mathrm{3}}$, which outperformed other models during the training ($R^{\mathrm{2}}$ $=$ 0.90; RSR $=$ 0.32; KGE $=$ 0.87; NSE $=$ 0.87; PBIAS $=$ 0.03) and testing ($R^{\mathrm{2}}$ $=$ 0.78; RSR $=$ 0.47; KGE $=$ 0.82; NSE $=$ 0.71; PBIAS $=$ $-$0.31) period, therefore was chosen to simulate streamflow in the Sutlej River in the 2050s and 2080s under the SSP245 and SSP585 scenarios. Bias correction was further applied to the projected daily streamflow in order to generate a reliable times series of the discharge. The mean ensemble of the model results shows that the mean annual streamflow of the Sutlej River is expected to rise between 2050s and 2080s by 0.79 % to 1.43 % for SSP585 and by 0.87 % to 1.10 % for SSP245. In addition, streamflow will increase during the monsoon (9.70 % to 11.41 % and 11.64 % to 12.70 %) in the 2050s and 2080s under both emission scenarios, but it will decrease during the pre-monsoon ($-$10.36 % to $-$6.12 % and $-$10.0 % to $-$9.13 %), post-monsoon ($-$1.23 % to $-$0.22 % and $-$5.59 % to $-$2.83 %), and during the winter ($-$21.87 % to$-$21.52 % and $-$21.87 % to $-$21.11 %). This variability in streamflow is highly correlated with the pattern of precipitation and temperature predicted by CMIP6 GCMs for future emission scenarios and with physical processes operating within the catchment. Predicted declines in the Sutlej River streamflow over the pre-monsoon (April to June) and winter (December to March) seasons might have a significant impact on agriculture downstream of the river, which is already having problems due to water restrictions at this time of year. The present study will therefore assist in strategy planning to ensure the sustainable use of water resources downstream by acquiring knowledge of the nature and causes of unpredictable streamflow patterns.

1 Introduction

Human-induced global warming has altered the patterns of rainfall worldwide (Goswami et al., 2006; Trenberth, 2011) and has also increased the risks of extreme events such as the droughts and floods (Easterling et al., 2000; Trenberth et al., 2015; Otto et al., 2017). It has impacted hydrology of many river basins globally, including the variation in streamflow (Gerten et al., 2008; Nepal and Shrestha, 2015; Singh et al., 2015a; Ali et al., 2018; Lutz et al., 2019; Singh et al., 2022). A study of the long-term (1948–2004) streamflow (discharge) data of the 200 largest rivers of the globe showed considerable changes in their annual discharge; however, the results were statistically significant only for 64 rivers (Dai et al., 2009). Out of these, 45 were marked with decreasing trends, and the remaining 19 showed increasing trends in their annual discharge. Similar decreasing and increasing trends in discharge of the rivers were also reported at the regional scale, i.e. Asia (Kundzewicz et al., 2009; Krysanova et al., 2015), Europe (Stahl et al., 2010; Stahl et al., 2012), and North America (Pasquini and Depetris, 2007). Moreover, it has been established that the effects of rainfall variation and extreme events on annual discharge are likely strong compared with other drivers (Kundzewicz et al., 2009; Miller et al., 2012; Van der Wiel et al., 2019). Zhao et al. (2021) examined how precipitation, evapotranspiration, and the timing of snowmelt impacted runoff in the Kaidu River basin in China. They discovered that, as global warming increased, the timing of snowmelt became less significant, while the influence of precipitation increased comparatively. A projected rise of $\sim$ 2 to 5 ${}^{\circ }$C in the mean annual global temperature by 2100 under higher greenhouse gas emission scenarios, as predicted by the general circulation models (GCMs; Gao et al., 2017), will considerably affect the rainfall pattern (intensity and amount) and may alter hydrological cycles (Oki and Kanae, 2006; Haddeland et al., 2014). This could subsequently impact the availability of water resources and present challenges for their management, since a rise in the demand of water is also predicted (Lutz et al., 2019). Therefore, it is indispensable to know the underlying hydrological dynamics occurring within a basin in the context of climate change for effective management and sustainable use of the water resources.

The underlying hydrological processes controlling rainfall–runoff generation in a basin can be understood with the use of a hydrological model which is based on complex mathematical equations and theoretical laws governing physical processes in the basin (Kirchner, 2006; Singh et al., 2019). It simulates/or predicts the response of the basin to climatological forcings such as rainfall (Sood and Smakhtin, 2015) and generates a synthetic time series of hydrological data that can be used by water managers and scientists for varied applications ranging from water budgeting and partitioning (Conan et al., 2003; Schreiner-McGraw and Ajami, 2020) to inundation mapping and modelling (Mahato et al., 2022). A hydrological model is supposed to not only have a good predictive power but also the ability to capture relationships among the forcing factors and catchment response so that an accurate estimate of the rainfall–runoff could be made (Shortridge et al., 2016). However, until now, there has been no hydrological model that can simulate basin behaviour universally well against all the hydrological challenges inflicted by climate change and human interventions (Yang et al., 2019). As a result, many hydrological models have been devised, considering the functioning and robustness of models for explaining the underlying complexity in quantifying the basin-scale response to the small-scale spatial complexity of physical processes (Shortridge et al., 2016; Herath et al., 2021). Broadly, these can be grouped into two categories, i.e. physical or process-based models and empirical or data-driven models (Yang et al., 2019; Kabir et al., 2020). The latter category of models uses a mathematical relationship established between runoff and affecting factors in the basin for deriving the runoff (Adnan et al., 2019).

It is purported that the data-driven model, despite the inherited limitations in the physical interpretability of the processes, has outperformed the physical models in terms of the prediction accuracy in many hydrological applications (Shortridge et al., 2016; Adnan et al., 2019; Kabir et al., 2020; Herath et al., 2021). Also, data-driven models are preferred over the physical models for rainfall–runoff modelling and/or streamflow prediction modelling due to the limited requirements of data as input, where the data limitation is the major challenge (Beven, 2011). These models, in the past, were heavily criticised on the grounds of being incompetent to model the non-linear behaviour of streamflow (Yang et al., 2019). But recent developments in computational intelligence, in the areas of machine learning (ML) and deep learning (DL) in particular, have greatly expanded the capabilities of empirical modelling (Adnan et al., 2020; Fu et al., 2020; Rahimzad et al., 2021; Ghobadi and Kang, 2022). This has resulted in the development of many non-linear models such as the artificial neural network (ANN), random forest (RF), support vector regression (SVR), and long short-term memory (LSTM) models, which can capture and model the non-stationarity of the rainfall–runoff relationships (Yaseen et al., 2015; Shortridge et al., 2016; Adnan et al., 2019; Yang et al., 2019; Xiang et al., 2020). Yang et al. (2019) applied three machine learning models, namely ANN, SVR, and RF, to predict monthly streamflow over the Qingliu River basin in China, under changing environmental conditions between 1989 and 2010, and compared their results with the six process-based hydrological models. They concluded that the ML model performed better than the process-based model not only in terms of prediction accuracy but also in terms of flexibility when it came to including other runoff-effect factors into the model. Similar outcomes for Lake Tana and the adjacent rivers in Ethiopia were also reported by Shortridge et al. (2016), where ML models demonstrated noticeably lower streamflow prediction errors than the physical models developed for the region. However, they inferred that linear machine learning models, such as the multivariate adaptive regression splines (MARSs) and generalised additive model (GAM), were sensitive to extreme climate events, so the degree of uncertainty in their predictions needed to be carefully considered.

The limitations of such data-driven models can be overcome by adopting more advanced ML and DL models (Xiang et al., 2020). Rasouli et al. (2012) compared the performance of the multilinear regression (MLR) model with the Bayesian neural network (BNN), SVR, and Gaussian process (GP) in terms of the daily streamflow prediction for the Stave River, a mountainous basin, in British Columbia, Canada, and found that the BNN model performed better than others. According to Hussain and Khan (2020), the supervised-learning model RF outperformed the multilayer perceptron (MLP) and SVR, in terms of accuracy, while predicting monthly streamflow for the Hunza River in Pakistan by 33.6 % and 17.85 %, respectively. Recently, deep neural network (DNN), convolutional neural network (CNN), and LSTM models, which are based on deep learning, have seen a surge in the number of streamflow prediction applications due to their abilities to handle complex stochastic datasets and abstract the internal physical mechanism (Fu et al., 2020; Ghobadi and Kang, 2022). Based on statistical performance evaluation criteria, Rahimzad et al. (2021) found that the LSTM outperformed the LR, SVR, and MLP models in daily streamflow prediction over the Kentucky River basin in the USA. However, Van et al. (2020) showed that CNN outperformed LSTM in streamflow modelling in the Vietnamese Mekong Delta by a small margin. Comparing data-driven models to a given problem yields a range of results for distinct geographical and climatic conditions (Hagen et al., 2021). Adnan et al. (2020) examined the predictive accuracy of the optimally pruned extreme learning machine (OP-ELM), least-square support vector machine (LSSVM), MARSs, and model tree (M5Tree) models in order to estimate monthly streamflow in the Swat River basin (Hindukush Himalaya), Pakistan. They came to the conclusion that the LSSVM and MARSs are the most effective at forecasting streamflow. In contrast, Hussain et al. (2020) discovered that ELM outperformed 1D-CNN while forecasting streamflow on three timescales, i.e. daily, weekly, and monthly in the Gilgit River, Pakistan. This suggests that it is challenging to find a data-driven model that is effective across all application domains and scales (Yaseen et al., 2015; Fu et al., 2020).

The use of machine-learning- and deep-learning-based models for streamflow simulations within catchments is generally limited to observable periods and the resulting forecasts (Eng and Wolock, 2022). There are very limited studies worldwide in which these models were applied for predicting the long-term streamflow for future periods in the context of climate change (Das and Nanduri, 2018; Thapa et al., 2021; Adib and Harun, 2022). This can be attributed to the challenges associated with data assimilation brought on by the use of coarse-resolution-scenario data obtained from general circulation models (GCMs), which limits their direct application in regional impact assessment (Hagen et al., 2021; Adib and Harun, 2022). Das and Nanduri (2018) integrated relevance vector machine (RVM) and support vector machine (SVM) models with the Coupled Model Intercomparison Project Phase (CMIP5) GCMs to project monthly monsoon streamflow across the Wainganga basin (India) for monsoon season. Adib and Harun (2022) studied the variations in the monthly streamflow pattern of the Kurau River (Malaysia) from 2021 to 2080 by coupling two ML models (RF and SVR) with the Coupled Model Intercomparison Project Phase (CMIP6) GCMs. Despite the significance potential of the ML and DL models in streamflow prediction, relevant studies assessing the application of these models for streamflow prediction under future scenarios over the mountainous basins are limited due to non-availability of long-term data (Xenarios et al., 2019; Adnan et al., 2020). Thapa et al. (2021) used a combination of the LSTM model and the CMIP5 GCM scenarios to estimate streamflow patterns in the Langtang basin of the central Himalaya. Their analyses revealed a notable increase in streamflow as a result of the predicted increase in precipitation. The projections from Coupled Model Intercomparison Project Phase 3 (CMIP3) GCMs and CMIP5 GCMs inherit limitations in the simulation of extreme precipitation (Kim et al., 2020), which are the principal drivers for the runoff generation in the catchment. This causes large uncertainty in streamflow predictions (Wang et al., 2021). Uncertainty in streamflow prediction can be minimised by using scenarios from the CMIP6 GCMs, which are likely to be more realistic than previous generations, i.e. CMIP3 GCMs and CMIP5 GCMs, given their significant improvement in simulating rainfall and temperature for historical records (Chen et al., 2020; Gusain et al., 2020; Kim et al., 2020). Therefore, projected changes in streamflow patterns derived from the CMIP6 GCM scenarios would give a better understanding of the catchment's future hydrological regime than previous ones. To the authors' knowledge, no work has been published about a mountainous basin that integrates ML and/or DL models with CMIP6 GCMs scenarios to predict changes in streamflow patterns for future periods. Hence, it is important to test whether machine learning approaches can be effectively used over a mountainous river basin to predict streamflow using hydrometeorological variables and CMIP6 GCM scenarios as the input data.

With a catchment area of 56 874 km${}^{\mathrm{2}}$ (up to Bhakra Dam), the Sutlej also pronounced as “Satluj”, is an important river in the western Himalaya and runs through diverse climatic zones. The flow in the upper and middle catchment is primarily impacted by glacier melt and snowmelt induced by a seasonal temperature shift and the preceding winter precipitation, while the lower section of the catchment area is mostly regulated by rainfall, both in the winter and during the monsoon season (Singh and Jain, 2002; Archer, 2003; Miller et al., 2012). Based on data from the period 1986–1996, Singh and Jain (2002) estimated the mean yearly contribution of snowmelt and glacier melt and rainfall to the Sutlej River to be 59 % and 41 %, respectively. However, the discharge in the river peaks is directly related to the peak in rainfall during the monsoon (Lutz et al., 2014). Recent studies on this basin has raised concerns about the implications of climatic changes on streamflow, since a warming climate has brought changes in the amount and spatiotemporal distribution of precipitation (Singh et al., 2014, 2015b). To date, previous research has only used process-based hydrological models and scenarios from CMIP3 GCMs and CMIP5 GCMs when examining the effects of climate change (past and future) on streamflow patterns in the region (Singh and Jain, 2002; Singh et al., 2015a; Ali et al., 2018; Shukla et al., 2021), which leaves a gap in the use of machine and deep learning models and scenarios from the latest CMIP6 GCMs. This study, for the very first time, examines the potential of five ML models and one DL model, namely the Gaussian linear regression model (GLM), Gaussian generalised additive model (GAM), MARSs, ANN, RF, and 1D-CNN, in streamflow prediction over the middle Sutlej River basin (rainfall-dominated zone) in western Himalaya using different Shared Socioeconomic Pathway (SSP) scenarios from CMIP6 GCMs. The pattern of variations in the Sutlej River's monthly, seasonal, and annual streamflow is assessed for the future periods of 2041–2070 (2050s) and 2071–2100 (2080s), with respect to the reference period of 1979–2009, under SSP245 and SSP585. The findings of the study will help to develop a better plan for the operation of hydroelectric power projects and water resources management in the catchment.

2 Study area

The selected study area is a sub-catchment within the Satluj basin (Fig. 1), with an area of 2457 km${}^{\mathrm{2}}$. Topographically, it is very rugged (0–80${}^{\circ }$) and is dominated mostly by forests (56.20 %), grassland (26.4 %), agricultural lands (17.1 %), and glaciers and snow cover (0.3 %; Singh et al., 2015a). The presence of mountain barriers in the sub-basin's north, the large variation in altitudes (500–5000 m), and the aspect all contribute to the region's diverse climate. It varies from a hot and moist tropical climate in the lower valleys to a cool temperate climate at about 2000 m and tends towards an alpine climate as the altitude increases beyond 2000 m. The mean annual discharge (averaged over the period of 1979–2009) of the river gauged at Kasol was 12 469.43 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$. There is a large inter-diurnal and monthly variation in the pattern of the river discharge. The minimum and maximum daily discharge recorded at Kasol was 64.30 and 2891 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, respectively. The early months of the year, i.e. starting from January up to March, are characterised by low streamflow. After this, a continuous and rapid rise in flow occurs, with the maximum being in the month of July ($\sim$ 22 %–23 %). Then, it again starts decreasing, and the flow reaches its minimum in the month of December (2 %–3 %). The details of the sub-catchment are summarised in Table 1.

Figure 1

The location of the sub-catchment within Sutlej River basin. The three hydrometeorological stations (Kasol, Sunni, and Rampur), from which this study employed observed data for the years 1979 to 2009, are also shown.

[Figure omitted. See PDF]

Table 1

Characteristics of the study catchment over the evaluation period of 1979–2009.

Parameters	Details
Details of the sub-catchment
Drainage area of the sub-catchment (km${}^{\mathrm{2}}$)	2457 km${}^{\mathrm{2}}$
Altitude	$\sim$ 500–5000 m
Slope	0–80${}^{\circ }$
Geology	Granite, Jutogh formation, and Chail, Salkhala, and Hemanta formations
Soil	Dystric Cambisols, dystric Regosols, and eutric Fluvisols
Streamflow measured at the outlet (Kasol) of the sub-catchment
Average annual streamflow	411.2 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$
Minimum streamflow (daily)	64.3 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$
Maximum streamflow (daily)	2891 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$
Standard deviation (SD) of the annual streamflow	1750.7 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$
Coefficient of variation (CV) of the annual streamflow	0.1 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$
Rainfall integrated over the sub-catchment
Average annual rainfall	1001.3 mm
Average monsoon rainfall (July–September)	403.0 mm
Average winter rainfall (December–March)	277.3 mm
Temperature integrated over the sub-catchment
Average annual maximum temperature ($T_{max}$)	28.3 ${}^{\circ }$C
Average annual minimum temperature ($T_{min}$)	13.9 ${}^{\circ }$C

The sub-basin is bestowed with large hydropower potential. There are three major hydroelectric power projects, namely the Sunni Dam Hydro Electric Project of 1080 MW,Rampur hydroelectric power plant (RHEP) of 412 MW, and Nathpa Jhakri India Hydroelectric Power Plant (NJHEP) of 1500 MW. The sub-basin is climatologically sensitive and, at present, facing challenges due to climate change and human interventions (Singh et al., 2015b, c). Changes in future climate will alter the patterns of flow in the river and could further affect the water resources and hydroelectric power production (Singh et al., 2014).

3 Description of the data and methods

The methodology involved in predicting streamflow for the period 2041–2100 in the Sutlej River is included in Sect. 3.1, which is a collection of hydrometeorological data, Sect. 3.2, which is a selection of machine and deep learning models, Sect. 3.3, which is a performance evaluation of the developed models, and Sect. 3.4, which is a bias correction in the streamflow projection. These are described in detail under following sub-headings.

3.1 Hydrometeorological data

The daily rainfall, temperature ($T_{max}$ and $T_{min}$), relative humidity, solar radiation, wind speed, and discharge data used to study the performance of the different machine and deep learning models on streamflow modelling were collected for 31 years, i.e. 1979–2009. Rainfall, temperature, and discharge data were obtained from the Bhakra Beas Management Board (BBMB), while relative humidity, solar radiation, and wind data were extracted from the Climate Forecast System Reanalysis (CFSR) Global Weather Data (http://globalweather.tamu.edu/, last access: 4 October 2020). These data were collected for three hydrometeorological stations, namely Kasol, Sunni, and Rampur (Fig. 1).

The downscaled outputs from the CMIP6 GCMs, the latest generation of climate models, were used for streamflow prediction in the future (2050s and 2080s). This framework of CMIP6 GCMs was run to simulate future climate under four Shared Socioeconomic Pathway (SSP) scenarios, which are designed to explain potential future greenhouse gas (GHG) emissions under various global socioeconomic shifts that could occur by 2100 (Riahi et al., 2017; Karan et. al, 2022). Even when using downscaled outputs, however, the regional climate change projections inherit biases from the GCM boundary conditions (Jose and Dwarakish, 2022), which were corrected in the dataset detailed in Mishra et al. (2020) for South Asia. They used the empirical quantile mapping (EQM) method to remove bias in the downscaled data. This dataset provides bias-corrected downscaled climate change projections for 13 CMIP6 GCMs and four GHG emission scenarios (SSP126, SSP245, SSP370, and SSP585), and the latter are briefly summarised in Riahi et al. (2017). Climate projections from CMIP6 GCMs that have been generated under the SSP245 and SSP585 scenarios are used in this study. SSP245, a medium scenario, represents the average pathway of future GHG emissions, with a radiative forcing of 4.5 W m${}^{-\mathrm{2}}$ by the year 2100, while SSP585 is the scenario with an upper limit in the range of scenarios with a radiative forcing of 8.5 W m${}^{-\mathrm{2}}$ by the end of this century (O'Neill et al., 2016). The data are available at a daily timescale and horizontal spatial resolution of 0.25${}^{\circ }$ $\times$ 0.25${}^{\circ }$. Seven grids of the downscaled CMIP6 GCM data cover the study area. The temperature ($T_{max}$ and $T_{min}$) data were adjusted for topographical bias by separating the study area into a number of homogenous elevation bands spaced by at an interval of 1000 m and applying a temperature lapse rate of 6.5 ${}^{\circ }$C per 1000 m within each grid. A digital elevation model (DEM) of 30 m spatial resolution derived from CartoSat-1 stereo data (http://www.bhuvan.nrsc.gov.in​​​​​​​, last access: 20 November 2020) was used for this purpose. The values of rainfall and temperature at each grid were then averaged over the catchment, using the Thiessen polygon method, in order to provide daily rainfall data integrated at the catchment scale for assessing changes in the future climate with respect to the observed period, i.e. 1979–2009.

Figure 2

Taylor diagram showing comparative skills of 13CMIP6 GCMs in simulating climatic variables (rainfall, $T_{max}$, and $T_{min}$) over the Sutlej sub-basin during reference period (1979–2009). The degree of the correlation coefficient ($r$) between the observed and CMIP6 GCMs, centred root mean square error (CRMSE), and departure of the models' standard deviation (SD) from the observed data (dashed black arc line) are shown in panel (a) for rainfall, panel (b) for $T_{max}$, and panel (c) for $T_{min}$. The units of SD for rainfall and temperature are in centimetres and degrees Celsius, respectively.

[Figure omitted. See PDF]

Furthermore, the ranking of CMIP6 GCMs was done to find the most appropriate models that can generate the most plausible scenarios of future climate in the catchment and, ultimately, be employed in the streamflow projection. A Taylor diagram (Taylor, 2001), a robust graphical plot, is widely used to rank GCMs due to its effectiveness in determining the relative strengths of the competing models and in evaluating overall performance as a model evolves (Abbasian et al., 2019; Ghimire et al., 2021). It integrates three statistical metrics, including the degree of correlation ($r$), centred root mean square error (CRMSE), and ratio of spatial standard deviation (SD). Combining these metrics allows the determination of the degree of pattern correspondence and explains how exactly a model represents the observed climate (Taylor, 2001). Therefore, the performance of 13 CMIP6 GCMs in modelling climatic variables (rainfall, $T_{max}$, and $T_{min}$) in the Sutlej sub-basin was compared to the observed data (1979–2009) using a Taylor diagram (Fig. 2a–c). The models were then ranked as a result of this comparison. A high positive correlation ($r$ $=$ 0.84 to 0.96) and low CRMSE ($<$ 3 ${}^{\circ }$C) error were found in all 13 CMIP6 GCMs for temperature ($T_{max}$ and $T_{min}$; Fig. 2b–c). Additionally, it was found that the models' standard deviations, which ranged from 5.60 to 6.03 ${}^{\circ }$C for $T_{max}$ and 6.34 to 6.63 ${}^{\circ }$C for $T_{min}$, were close to the SD of the observed data (6.01 and 6.07 ${}^{\circ }$C). These results imply that all CMIP6 GCMs may be able to predict the most likely future temperature over the catchment.

However, not all CMIP6 GCMs showed the high degree of similarity in predicting rainfall; in fact, 2 (CanESM5 and NorESM2-LR) of the 13 models revealed a negative correlation (Fig. 2a). In the pool of 13 CMIP6 GCMs, only six models showed relatively higher correlation ($r$ $\ge$ 0.56), smaller CRMSE ($<$ 12 cm) errors, and a high similarity to the standard deviation of the observed data (13.2 cm). They were the (1) Earth Consortium Earth 3 Veg model (EC-Earth3-Veg), (2) Russian Institute for Numerical Mathematics climate model version 4.8 (INM-CM4-8), (3) Russian Institute for Numerical Mathematics climate model version 5.0 (INM-CM5-0), (4) Max Planck Institute for Meteorology Earth System Model version 1.2 with higher resolution (MPI-ESM1-2-HR), (5) Max Planck Institute for Meteorology Earth System Model version 1.2 with lower resolution (MPI-ESM1-2-LR), and (6) Norwegian Earth System Model version 2 with medium resolution (NorESM2-MR). Furthermore, within these models, the highest and lowest correlations between observed and simulated rainfall were found for the INM-CM4-8 ($r$ $=$ 0.69) and NorESM2-MR ($r$ $=$ 0.56), respectively. These six CMIP6 GCMs were finally selected to examine future patterns in streamflow for the periods 2050s and 2080s in the Sutlej River basin, as they had also shown high performance in simulating temperatures ($r$ $=$ 0.90 to 0.96).

3.2 Selection of machine learning and deep learning models for streamflow modelling

In this study, five machine learning models and one deep learning model, namely GLM, GAM, MARSs, ANN, RF and 1D convolution neural network (1D-CNN), were selected, and their performances in predicting streamflow in Sutlej River were compared. These are regression-based models which capture the relationship between the predictors (dependent variables) and the predictand (independent variables) and provide the value of the output variables (Adnan et al., 2019; Kabir et al., 2020). The models were trained with daily observed data recorded during 1979–2009 at Kasol (the gauging site) and simulated historical projections of CMIP6 GCMs. The climatic projections of the grid corresponding to the Kasol station were taken into consideration as the input from the CMIP6 GCMs. However, prior to building the models, all of the data were normalised, using standard normalisation techniques, to standardise the features on a common scale. Furthermore, the entire dataset was split into training and testing datasets, since a cross-validation method was adopted in this study. The training dataset (80 %) was used for fitting the models, whereas a testing dataset was used for checking model accuracy (20 %). Under the cross-validation method, the process was repeated until every part of the allocated data was used in testing (Kabir et al., 2020). Six different program codes were written in the Python language for ANN, GAM, GLM, MARSs, RF, and 1D-CNN simulations. Out of these six selected models, GLM, GAM, and MARSs are linear models, whereas other three i.e. ANN, RF, and 1D-CNN, are non-linear in nature (Shortridge et al., 2016; Yang et al., 2019; Herath et al., 2021). Additionally, excluding GLM, all of the remaining models are based on a non-parametric regression approach, where the functional relationship between the predictor and predictand is not predetermined but can be adjusted to capture the unusual or unexpected features of the data (Shortridge et al., 2016). A detailed description of these models can be found elsewhere (Shortridge et al., 2016; Adnan et al., 2019; Yang et al., 2019; Kabir et al., 2020; Ghimire et al., 2021; Herath et al., 2021; Shu et al.,2021).

Since the 1D-CNN model is based on weight sharing, it needs fewer training parameters than other models (Kiranyaz et al., 2021). It has mainly three layers, i.e. a convolution layer, pooling layer, and fully connected layer. The primary job of the convolution layer is to non-linearly map input data into a set of feature maps or a series of feature vectors. When working as a visual cortical perceptron, the filter kernels are convoluted with the input data of their receptive fields. The convolution results with biases are then passed on to the activation function to create feature maps. The pooling layer, which comes after each convolution layer, primarily serves to reduce the dimension of the feature maps and maintain the invariance of characteristic scale. The fully connected layer uses a completely connected single layer perceptron to combine the feature maps that were acquired by the prior convolution and pooling layers in order to build a higher-level feature (Kiranyaz et al., 2021). In this study, one convolution layer with 64 filters, a kernel of size 2, and a rectified linear activation function (ReLU) was employed. This was followed by max pooling layer with pool size of 2 and the faltterm layer. After that, two fully connected layers are applied with ReLU and a linear activation function, respectively. However, for optimisation, the adaptive moment estimation (Adam) algorithm was applied (Ghimire et al., 2021; Shu et al., 2021). Six variables, namely rainfall, $T_{max}$, $T_{min}$, relative humidity, solar radiation, and wind speed were used as input for developing the models. Additionally, to understand the control of the catchment size and geo-hydromorphological characteristics of the basin in generating runoff, these models were simulated under the following four rainfall scenarios: rainfall on the same day ($R_{\mathrm{0}}$), rainfall lagged by 1 d ($R_{\mathrm{1}}$), rainfall lagged by 2 d ($R_{\mathrm{2}}$), and rainfall lagged by 3 d ($R_{\mathrm{3}}$). The remaining meteorological parameters were held constant during the processes.

Table 2

The information on hyperparameters used for estimating model parameters.

Model name	Hyperparameter	Values
Artificial neural network (ANN)	build_fn warm_start random_state optimizer loss metrics batch_size validation_batch_size verbose callbacks validation_split shuffle run_eagerly epochs	value $=$ build_regressor value $=$ False value $=$ None value $=$ rmsprop value $=$ None value $=$ None value $=$ 64 value $=$ None value $=$ 1 value $=$ None value $=$ 0.0 value $=$ True value $=$ False value $=$ 500
Generalised additive model (GAM)	formula family data weights subset na.action offset method optimizer control scale select knots sp min.sp H gamma fit paraPen G drop.unused.levels drop.intercept discrete	value $=$ None value $=$ gaussian() value $=$ list() value $=$ Null value $=$ Null value $=$ Null value $=$ “GCV.Cp” value $=$ c(“outer”, “newton”) value $=$ list(), value $=$ 0 value $=$ False value $=$ Null value $=$ Null value $=$ Null value – Null, value $=$ 1 value $=$ True value $=$ Null value $=$ Null value $=$ True value $=$ Null value $=$ False
Generalised linear model (GLM)	endog exog family offset exposure freq_weights var_weights missing	value $=$ 1D value $=$ 1D value $=$ sm.families.Gaussian(sm.families.links.log()) value $=$ None value $=$ None value $=$ None value $=$ None value $=$ str
Multivariate adaptive regression splines (MARSs)	max_terms max_degree allow_missing penalty endspan_alpha endspan minspan_alpha minspan thresh	value $=$ 20 value $=$ 3 value $=$ False value $=$ 3.0 value $=$ 0.005 value $=$ $-$1 value $=$ 0.005 value $=$ $-$1 value $=$ 0.001

Table 2

Continued.

Model name	Hyperparameter	Values
	zero_tol min_search_points check_every allow_linear use_fast fast_K fast_h smooth enable_pruning feature_importance_type feature_importance_type	value $=$ 1 $\times$ 10${}^{-\mathrm{12}}$ value $=$ 100 value $=$ $-$1 value $=$ True value $=$ False value $=$ 5 value $=$ 1 value $=$ False value $=$ True value $=$ None value $=$ 0
Random forest (RF)	n_estimators criterion max_depth min_samples_split min_samples_leaf min_weight_fraction_leaf max_features max_leaf_nodes min_impurity_decrease	value $=$ 500 value $=$ “squared_error” value $=$ None value $=$ 2 value $=$ 5 value $=$ 0.0 value $=$ auto value $=$ None value $=$ 0.0
1D convolution neural network (1D-CNN)	Conv1D_filter Conv1D_kernel_size Conv1D_pool_size Learning rate Epoc Batch size loss	value $=$ 64 value $=$ 2 value $=$ 2 value $=$ 0.0001 value $=$ 30 value $=$ 280 value $=$ MSE

3.3 Model performance evaluation

It has been found that overfitting in a model may lead to large errors in out-of-sample predictions (Hastie et al., 2009). Therefore, it has been evaded by establishing model parameters for GLM, GAM, MARSs, ANN, and RF through automated hyperparameter tuning methods. In total, 500 bootstrapped resamples of the training dataset were generated for each parameter value to be assessed. Table 2 presents the information on the specific parameters evaluated for each model.

Figure 3

Evaluation of the models' (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) performance in simulating streamflow under rainfall scenarios (a) $R_{\mathrm{0}}$, (b) $R_{\mathrm{1}}$, (c) $R_{\mathrm{2}}$, and (d) $R_{\mathrm{3}}$ at Kasol during the training phase using six statistical metrics ($R^{\mathrm{2}}$, KGE, NSE, RSR, MAE, and PBIAS).

[Figure omitted. See PDF]

Figure 4

Evaluation of the models' (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) performance in simulating streamflow under rainfall scenarios (a) $R_{\mathrm{0}}$, (b) $R_{\mathrm{1}}$, (c) $R_{\mathrm{2}}$, and (d) $R_{\mathrm{3}}$ at Kasol during the testing phase using six statistical metrics ($R^{\mathrm{2}}$, KGE, NSE, RSR, MAE, and PBIAS).

[Figure omitted. See PDF]

The accuracy with which the simulated flow matches the observed flow during the training (calibration) and testing (validation) phases determines whether a hydrological model is appropriate for a given application (Refsgaard, 1997). Several methods, including quantitative statistics and graphical methods, have been developed in the past for assessing the accuracy of model predictions (Legates and McCabe, 1999). Moriasi et al. (2007) grouped these methods into three categories, namely standard regression, dimensionless, and error index, depending on how well each method explains the relationship between observed and simulated values, compares the relative performance of the models, and quantifies the deviation in the units of the data of interest. Moreover, it has been established from previous studies that a single metric is inadequate to evaluate a model's performance; hence, multiple metrics should be used (Adnan et al., 2020). Therefore, in this study, the prediction accuracy of different models was compared using six statistical measures, out of which one was a standard regression (coefficient of determination, $R^{\mathrm{2}}$), two of which were dimensionless (Kling–Gupta efficiency, KGE, and Nash–Sutcliffe efficiency, NSE), and the remaining three were the error index (ratio of the root mean square error to the standard deviation of the measured data, RSR, the mean absolute error, MAE, and the percent bias, PBIAS). These metrics are defined below by Eqs. (2)–(7): $$\begin{array}{}\text{1}& {\displaystyle }\begin{array}{rl}& R^{\mathrm{2}}=\left({\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left(Q_i-\stackrel{\mathrm{‾}}{Q}\right)\left(P_i-\stackrel{\mathrm{‾}}{P}\right)}{\sqrt{\sum _{i=\mathrm{1}}^n{\left(Q_i-\stackrel{\mathrm{‾}}{Q}\right)}^{\mathrm{2}}}\times \sqrt{{\left(P_i-\stackrel{\mathrm{‾}}{P}\right)}^{\mathrm{2}}}}}\right)\\ & \left(\text{range: 0 to 1}\right)(\mathrm{Van}\mathrm{Liew}\mathrm{et}\mathrm{al}.,\mathrm{2003})\end{array}\text{2}& {\displaystyle }\begin{array}{rl}& \mathrm{KGE}=\mathrm{1}-\sqrt{{\left(r-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{p}}}{{\mathit{\sigma }}_{\mathrm{ob}}}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{P_i}{Q_i}}-\mathrm{1}\right)}^{\mathrm{2}}}\\ & \left(\text{range: 0 to 1}\right)(\mathrm{Gupta}\mathrm{et}\mathrm{al}.,\mathrm{1999})\end{array}\text{3}& {\displaystyle }\begin{array}{rl}& \mathrm{NSE}=\mathrm{1}-\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(Q_i-P_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(Q_i-\stackrel{\mathrm{‾}}{Q}\right)}^{\mathrm{2}}}}\right]\\ & (\text{range:}-\mathrm{\infty }\mathrm{to}\mathrm{1})(\mathrm{Nash}\mathrm{and}\mathrm{Sutcliffe},\mathrm{1970})\end{array}\end{array}$$

$$\begin{array}{}\text{4}& {\displaystyle }\begin{array}{rl}& \mathrm{RSR}={\displaystyle \frac{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left(Q_i-P_i\right)}^{\mathrm{2}}}}{{\mathit{\sigma }}_{\mathrm{ob}}}}\\ & \left(\text{range: 0 to }\mathrm{\infty }\right)(\mathrm{Murphy}\mathrm{et}\mathrm{al}.,\mathrm{2004})\end{array}\text{5}& {\displaystyle }\begin{array}{rl}& \mathrm{MAE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left|P_i-O_i\right|}{n}}\\ & \left(\text{range: 0 to }\mathrm{\infty }\right)(\mathrm{Adnan}\mathrm{et}\mathrm{al}.,\mathrm{2020})\end{array}\text{6}& {\displaystyle }\begin{array}{rl}& \mathrm{PBIAS}=\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^n\left(Q_i-P_i\right)}{\sum _{i=\mathrm{1}}^n\left(Q_i\right)}}\right]\times \mathrm{100}\\ & (\text{range: }-\mathrm{100}\mathit{\%}\mathrm{to}\mathrm{100}\mathit{\%})(\mathrm{Gupta}\mathrm{et}\mathrm{al}.,\mathrm{1999}),\end{array}\end{array}$$ where $P_i$ are the predicted values, and $Q_i$ are the observed values. $n$ accounts for the number of samples, $\stackrel{\mathrm{‾}}{Q}$ represents the mean of the observed data, and $\stackrel{\mathrm{‾}}{P}$ is the mean of the predicted data. However, $r$ is the Pearson's correlation coefficient, whereas ${\mathit{\sigma }}_{\mathrm{ob}}$ and ${\mathit{\sigma }}_{\mathrm{p}}$ refer to the standard deviation of the observed and predicted values, respectively.

Figure 5

Comparison of observed and simulated streamflow for all six models (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) under rainfall scenarios $R_{\mathrm{0}}$.

[Figure omitted. See PDF]

$R^{\mathrm{2}}$ evaluates the percentage of the variation in the measured data that can be explained by the model, whereas the NSE estimates the relative size of the residual variance in relation to the variance in the measured data (Nash and Sutcliffe, 1970; Van-Liew et al., 2003). According to Mazrooei et al. (2021), the NSE is sensitive to extreme flows; as a result, the KGE is also used to evaluate a model's performance, while taking extreme flows into account (Adib and Harun, 2022). Other metrics, like RSR, MAE, and PBIAS, shed light on the overall inaccuracies in the projected flow relative to the observed. The values of $R^{\mathrm{2}}$, KGE, and NSE should all be 1 in an ideal model, whereas RSR, MAE, and PBIAS values should be 0 (Nash and Sutcliffe, 1970; Van-Liew et al., 2003; Gupta et al., 1999; Adnan et al., 2020). Moriasi et al. (2007) developed a guideline for interpreting the results of these metrics and a rank for the hydrological models based on a thorough review of the available literature. They found that a model can be classified as very good, good, satisfactory, or unsatisfactory if its NSE value is between 0.75 and 1, 0.65 and 0.75, 0.50 and 0.65, or less than 0.50, respectively. Similarly, $R^{\mathrm{2}}$ values between 0.6 and 0.7 are considered satisfactory, 0.85 and 1 are very good, and numbers below 0.5 are unsatisfactory (Van-Liew et al., 2003). However, for RSR, numbers above 0.7 are considered to be poor, whereas values between 0 and 0.5 are considered to be in the very good range. Thus, the lower the RSR value, the better the model. This is also true for PBIAS and MAE, where lower values are favourable. According to Moriasi et al. (2007), PBIAS values of less than $\pm$10 % are considered to be highly acceptable, while values of more than $\pm$25 % are considered to be unsatisfactory. The negative number indicates that the model has overestimated its bias, whereas the positive value indicates that the model has underestimated its bias (Gupta et al., 1999).

Figure 6

A comparison of observed and simulated streamflow for all six models (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) under rainfall scenarios $R_{\mathrm{1}}$.

[Figure omitted. See PDF]

3.4 Bias correction

Uncertainty in streamflow prediction may be caused by the GCM shortcomings (e.g. coarse spatial resolution, simplified physics and thermodynamic processes, numerical methods, or poor knowledge of climate system dynamics) in accurately replicating natural climate variability (Sperna Weiland et al., 2010). As a result, its quantification and correction are critical for generating a future time series of streamflow that is reliable and recommended for devising water resource management plans in the catchment. This study used the bias correction method proposed in Hawkins et al. (2013) to correct the uncertainty (bias) between observed and CMIP6-GCM-predicted streamflow. The mathematical expression for this formula is given below:

7$$Q_{\mathrm{bc}}={\stackrel{\mathrm{‾}}{Q}}_{\mathrm{ob}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{ob}}}{{\mathit{\sigma }}_{\mathrm{p}}}}\left(Q_{\mathrm{future}}-{\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}\right),$$ where $Q_{\mathrm{bc}}$ and $Q_{\mathrm{future}}$ are the bias-corrected and raw daily discharge for future simulation, respectively.

${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{ob}}$ and ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{p}}$ are the mean discharge of observed and historical simulation for the reference period (1979–2009), respectively. ${\mathit{\sigma }}_{\mathrm{o}}$ and ${\mathit{\sigma }}_{\mathrm{p}}$ are the standard deviation in the observed and historical simulation for the reference period, respectively. This method captures the variability in both the observation and GCM simulations (Hawkins et al., 2013), which is the interest of this study.

Figure 7

Comparison of observed and simulated streamflow for all six models (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) under rainfall scenarios $R_{\mathrm{2}}$.

[Figure omitted. See PDF]

4 Results

4.1 Streamflow simulation and evaluation of model performance

The simulation (1979–2009) results generated under different rainfall scenarios ($R_{\mathrm{0}}$, $R_{\mathrm{1}}$, $R_{\mathrm{2}}$, and $R_{\mathrm{3}}$) on a daily timescale for all six models (GLM, GAM, MARSs, ANN, RF, and 1D-CNN) during training and testing are shown in Figs. 3 and 4, respectively. The model performed slightly better during training than during testing periods. $R^{\mathrm{2}}$, NSE, and KGE values across models ranged from 0.69 to 0.90, 0.52 to 0.87, and 0.69 to 0.91 and from 0.69 to 0.81, 0.49 to 0.74, and 0.68 to 0.82 during training and testing, respectively. Likewise, it was found that RSR, MAE, and PBIAS varied from 0.31 to 0.55, from 71.95 to 123.25 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and from $-$2.11 % to $+$4.31 % during training, as well as from 0.56 to 0.46, from 123.06 to 106.64 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and from $-$3.74 % to $+$2.21 % during testing, respectively. Non-linear models (ANN, 1D-CNN, and RF) outperformed linear models (GAM and GLM) in runoff prediction under all rainfall scenarios ($R_{\mathrm{0}}$, $R_{\mathrm{1}}$, $R_{\mathrm{2}}$, and $R_{\mathrm{3}}$), with the exception of MARSs, which produced results that were more or less comparable with those of the ANN model. Figures 3–4 show that both models (RF and 1D-CNN) satisfy the performance requirements outlined by Moriasi et al. (2007), but RF slightly outperformed CNN in terms of the error index. $R^{\mathrm{2}}$, NSE, KGE, RSR, MAE, and PBIAS values for the RF model during the training ranged from 0.88 to 0.90, 0.85 to 0.87, 0.86 to 0.87, 0.32 to 0.34, 71.95 to 77.49 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and $+$0.03 % to $+$0.13 %, respectively. For the 1D-CNN, however, it varied from 0.87 to 0.89, 0.85 to 0.87, 0.90 to 0.91, 0.34 to 0.35, 80.29 to 83.14 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and $-$1.25 % to $+$0.13 %. A similar pattern with slightly lower values was revealed during testing for the both models. This implies that RF can effectively capture non-linear interactions and can provide insights about actual watershed functions (Shortridge et al., 2016). On the other hand, GLM showed the poorest results. $R^{\mathrm{2}}$, NSE, KGE, RSR, MAE, and PBIAS values for the GLM during the training varied from 0.69 to 0.71, 0.52 to 0.56, 0.71 to 0.72, 0.54 to 0.55, 134.80 to 140.56 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and $+$2.63 % to $+$2.73 %, respectively. During testing, they varied between 0.69 and 0.71, 0.49 and 0.54, 0.68 and 0.70, 0.54 and 0.56, 134.35 and 141.26 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, and $+$1 % and $+$1.31 %, respectively. Furthermore, it was observed that the models with rainfall scenario $R_{\mathrm{3}}$ revealed reasonably better results in comparison to the $R_{\mathrm{0}}$, $R_{\mathrm{1}}$, and $R_{\mathrm{2}}$ scenarios, indicating the delayed contribution of rainfall–runoff to the river.

Figure 8

Comparison of observed and simulated streamflow for all six models (ANN, GAM, GLM, MARSs, RF, and 1D-CNN) under rainfall scenarios $R_{\mathrm{3}}$.

[Figure omitted. See PDF]

Figures 5, 6, 7, and 8 show a comparison of the observed and simulated streamflow under rainfall scenarios of $R_{\mathrm{0}}$, $R_{\mathrm{1}}$, $R_{\mathrm{2}}$, and $R_{\mathrm{3}}$ for all the models at Kasol, which is the outlet of the basin. As observed in Figs. 5–8, RF was able to follow the curve better compared to the other models. It is also deduced from the comparison of scatterplots wherein a relatively smaller deviation in the observed and estimated discharge of streamflow was found for the RF model. GLM performed the worst out of the six models with respect to the time variation graphs. A limitation faced by all the six models was the simulation of peak values. The models slightly underperformed at the prediction of higher values of streamflow. These findings led to the ultimate decision to use the RF model with rainfall scenario $R_{\mathrm{3}}$ to predict streamflow in the Sutlej River in the future (2050s and 2080s) under the SSP245 and SSP585 scenarios.

Figure 9

Box plot comparing observed and CMIP6-GCM-simulated (mean ensemble of models) streamflow for various months of the year, as derived over the period 1979–2009. The line inside the box denotes the median values of the streamflow, while the upper and lower whiskers indicate the highest and minimum values, respectively.

[Figure omitted. See PDF]

Figure 10

Probability exceedance curves developed using 10 % of the highest and lowest flows from the observed and CMIP6-GCM-simulated (mean ensemble of models) flows over the time span of 1979–2009 for annual and seasonal (pre-monsoon and monsoon) flows.

[Figure omitted. See PDF]

Figure 11

Projected change in the mean monthly rainfall in the sub-basin using different CMIP6 GCMs under the SSP245 and SSP585 scenarios in the (a, c) 2050s​​​​​​​ and (b, d) 2080s.

[Figure omitted. See PDF]

4.2 Comparison of streamflow simulated with observed and CMIP6 GCMs data

The uncertainty between observed and CMIP6-GCM-predicted streamflow during the reference period (1979–2009) was investigated by comparing the streamflow simulated by the RF model with the observed and CMIP6 GCM data. A large difference in streamflow patterns was seen in the box plot of observed and CMIP6-GCM-simulated discharge (Fig. 9) derived for various months of the year, particularly from June through September (monsoon season), when a pattern of intense daily rainfall was observed over the catchment. Additionally, it was discovered through the analysis of the probability exceedance curves generated using 10 % of the time series' highest flows that, despite the streamflow in the two datasets being comparable throughout the pre-monsoon season (Fig. 10c), they differ noticeably for high flows during the annual (Fig. 10a) and monsoon season (Fig. 10c). Similar trends were seen in the comparison of the probability exceedance curves for low flows during the monsoon season, although there was strong agreement for annual (Fig. 10b) and pre-monsoon measurements (Fig. 10d). This may be due to the fact that orography has a considerable impact on the regional Indian summer monsoon (ISM) climate, making it challenging for climate models to predict daily monsoonal rainfall accurately across the Himalaya (Turner and Annamalai, 2012; Niu et al., 2015; Choudhary et al., 2022). The regional climate model (RCM) based on CMIP5 GCMs was used by Sanjay et al. (2017) to study the pattern of change in precipitation and temperature over the Hindukush Himalaya region. As a condition of the model's inability to accurately represent complicated feedback mechanisms, the results revealed large uncertainty in the summer and winter precipitation over the northwestern Himalaya. This is also supported by the study of Kadel et al. (2018). They evaluated the performance of 38 CMIP5 GCMs in simulating rainfall over the central Himalaya and came to the conclusion that the majority of the models studied performed poorly when it came to reproducing the spatial distribution of monsoonal rainfall. Although the most recent study by Gusain et al. (2020) in India reported that an ISM simulation using CMIP6 GCMs over CMIP5 GCMs had significantly improved, there are discrepancies between the models, and this indicated uncertainty in the predictions. Lalande et al. (2021) examined the abilities of 26 CMIP6 GCMs to simulate the rate of precipitation across the Himalayan region and concluded that the models consistently overestimated the rate of precipitation by 31 % to 281 %. Additionally, a cold bias in temperature estimation was also reported. Therefore, bias correction, as described in Sect. 3.4, was applied to the projected streamflow for the future periods (2050s and 2080s) under all scenarios and for all six models in order to provide accurate times series of the discharge.

Figure 12

Projected change in the mean seasonal rainfall in the sub-basin using different CMIP6 GCMs under the SSP245 and SSP585 scenarios in the (a, c) 2050s and (b, d) 2080s.

[Figure omitted. See PDF]

Figure 13

Projected changes in the (a) mean annual rainfall, (b) $T_{max}$, and (c) $T_{min}$ in the sub-basin using different CMIP6 GCMs under the SSP245 and SSP585 scenarios in the 2050s and 2080s.

[Figure omitted. See PDF]

4.3 Projected change in rainfall and temperatures in 2050s and 2080s under SSP245 and SSP585

Figure 11 shows how the catchment's mean monthly rainfall is expected to change under SSP245 and SSP585 in the 2050s and 2080s compared to the reference period (1979–2009). Within the months and for the CMIP6 GCMs, a sizeable shift in the rainfall pattern is seen. With the exception of March, June, and September, the mean ensemble of the models generally predicts a rise in rainfall throughout the year in the 2050s and 2080s under all scenarios. The models also show a significant variation in the seasonal and yearly rainfall patterns expected for the catchment in the 2050s and 2080s under various emission scenarios. However, based on the mean ensemble of the models, it is predicted that seasonal (Fig. 12) and annual (Fig. 13a) rainfall will generally increase in the 2050s and 2080s under SSP245 and SSP585. Pre-monsoon, monsoon, post-monsoon, and winter rainfall in the 2050s will increase by 8.75 % to 8.85 %, 10 % to 20.80 %, 85 % to 91.91 %, and 12.48 % to 14.16 %, respectively, under SSP245 and SSP585. However, under SSP245 and SSP585 in the 2080s, it will rise by 7.69 % to 17.50 %, 21.52 % to 41.43 %, 56.16 % to 89.66 %, and 22.48 % to 12.43 %, respectively. Under both scenarios in the 2050s and 2080s, the pre-monsoon and post-monsoon will have the lowest and highest percentage increases in rainfall, respectively. The monsoon season is, however, anticipated to have the greatest rise in terms of quantity ($\sim$ 40–167 mm). The predicted range for the increase in mean annual rainfall is 13.85 % to 18.61 % in the 2050s and 17.91 % to 34.31 % in the 2080s. It is observed that the predicted pattern of change in rainfall across the sub-basin under various SSPs is consistent in terms of the direction of change with other studies conducted over the Sutlej and Himalaya region. Lalande et al. (2021) reported an overall increase in the mean annual precipitation over the Himalayan region based on 10 CMIP6 GCMs. According to their analysis, the mean ensemble of model precipitation is predicted to increase by 8.6 % to 25.4 % in 2081–2100 under SSP245 and SSP585. The same study also showed an increase in the region's winter (November to April) and ISM (June to September) rainfall. This contradicts past studies that showed a trend towards declining ISM rainfall after the 1950s (Sabin et al., 2020). They postulated that the region's higher winter rainfall would have been caused by the strengthening of the western disturbances; however, the intensification of the ISM is responsible for the region's enhanced summer rainfall.

The analysis of the CMIP6 GCM projections leads to the conclusion that, for all months and seasons in the 2050s and 2080s, maximum (excluding April and pre-monsoon in the 2050s under SSP245) and minimum temperatures will rise under both scenarios (Figs. 14a–d and 15a–d). Similarly, increase in mean annual $T_{min}$ and $T_{max}$ is also predicted in the 2050s and 2080s under all scenarios (Fig. 13b and c). The increase will be relatively higher for the $T_{min}$ compared to the $T_{max}$. This is also reported by Singh et al. (2015c). The increase in rainfall and temperature is typically higher under SSP585 than SSP245 in both eras (2050s and 2080s), as expected, due to a larger increase in radiative forcing brought on by increased greenhouse gas emissions.

Figure 14

Projected change in the mean seasonal maximum temperature ($T_{max}$) in the sub-basin using different CMIP6 GCMs under the SSP245 and SSP585 scenarios in the (a, c) 2050s and (b, d) 2080s.

[Figure omitted. See PDF]

Figure 15

Projected changes in the mean seasonal minimum temperature ($T_{min}$) in the sub-basin using different CMIP6 GCMs under the SSP245 and SSP585 scenarios in the (a, c) 2050s and (b, d) 2080s.

[Figure omitted. See PDF]

4.4 Assessment of the change in streamflow in 2050s and 2080s under SSP245 and SSP585

The Sutlej River's mean monthly streamflow change, compared to the reference period's observed flow (1979–2009), is shown in Fig. 16 under scenarios SSP245 and SSP585 for the future periods (2050s and 2080s). According to both scenarios and all six models, the Sutlej River's streamflow will decrease between January ($-$33.80 % to $-$14.38 %), February ($-$32.40 % to $-$14.15 %), March ($-$23.55 % to $-$0.84 %), November ($-$21.06 % to $-$5.14 %), and December ($-$29.88 % to $-$18.38 %) in the 2050s and 2080s. Moreover, except for MPI-ESM1-2-HR and MPI-ESM1-2-LR, which show an increase in streamflow in the 2080s under the higher-emission scenario, all of the CMIP6 GCMs indicate a decrease in the river's discharge in June ($-$20.24 % to $-$0.57 %) under SSP245 and SSP585 for both the periods. Similarly, excluding EC-Earth3-Veg (under SSP245 in 2050s) and INM-CM5-0 (under SSP245 in the 2050s and 2080s and under SSP585 in the 2050s), all of the CMIP6 GCMs indicate a decrease in the river's discharge in May ($-$25 % to $-$2.85 %) during the study period. In contrast, under SSP245 and SSP585 in the 2050s and 2080s, all of the CMIP6 GCMs predict a rise in the river's discharge in April (20.24 % to $-$0.57 %; excluding SSP585 in the 2080s), August (16.84 % to 5.28 %), and September (55.27 % to 4.35 %). But no clear pattern of streamflow change is seen for the remaining months (July and October) of the year, making results difficult to generalise because the projected decrease and/or increase in streamflow over the months is inconsistent among models under various emission scenarios in the 2050s and 2080s. The variations in climate-variable projections caused by differing spatial resolutions and parameterisation levels in the climate models may be the cause of these discrepancies in streamflow estimates (Sperna Weiland et al., 2010; Singh et al., 2015a). According to Murphy et al. (2004), the average of an ensemble of GCMs cancels out the errors in each individual model, and as more models are used, the ensemble uncertainty decreases. Therefore, in order to reduce the uncertainty in the projection of streamflow related to individual CMIP6 GCMs, the streamflow pattern of the Sutlej River was analysed by also using the mean ensemble of all six GCMs.

Figure 16

Predicted change in monthly streamflow pattern of the Sutlej River with respect to the reference period (1979–2009) in (a, b) 2050s and (c, d) 2080s under SSP245 and SSP585 scenarios for different CMIP6 GCMs.

[Figure omitted. See PDF]

Figure 17

Comparison of monthly observed (1979–2009) and projected discharge of the multimodel ensembles for the 2050s and 2080s under the SSP245 and SSP585 scenarios.

[Figure omitted. See PDF]

The mean ensemble of the models predicts that the Sutlej River's mean monthly streamflow (excluding April) will decrease under both scenarios from November ($-$18.45 % to $-$17.17 %) to June ($-$10.90 % to $-$8.06 %) between the 2050s and 2080s (Fig. 17). The river flow, which would have been expected to increase in April under both scenarios in 2050s, will also decline in the 2080s for the higher emission scenarios (SPP585). The maximum and minimum streamflow declines are predicted to occur in the 2050s under SSP245 for the months of December ($-$24.25 %) and May ($-$7.77 %), respectively. In comparison to SPP245, the decline will generally be slightly higher under SSP585 in 2050s, and for the 2080s, the projected decrease in streamflow will not show much difference under both the scenarios. Opposite to this, the mean ensemble of the models predicts that the Sutlej River's flow will increase from July (5.50 % to 5.91 %) to October (3.01 % to 11.42 %) in the 2050s and 2080s under both scenarios. The maximum and minimum streamflow increases are predicted to occur in the 2080s under SSP245 for the months of September (25.82 %) and July (5.50 %), respectively. In all scenarios, the increase will be slightly greater in the 2080s than it will be in the 2050s. When compared to SPP245, it will be higher for SSP585 in the scenarios.

Figure 18

Predicted change in seasonal streamflow pattern of the Sutlej River with respect to the reference period (1979–2009) in the (a, c) 2050s and (b, d) 2080s under the SSP245 and SSP585 scenarios for different GCMs.

[Figure omitted. See PDF]

The projected change in the seasonal streamflow of the Sutlej River in 2050s and 2080s is shown in Fig. 18. The 2050s and 2080s would see an increase in streamflow during the monsoon (4.46 % to 16.14 %) and a decrease during the pre-monsoon ($-$17.40 % to $-$0.51 %) and winter ($-$28.81 % to $-$12.42 %) for all six CMIP6 GCMs, with the exception of INM-CM5-0 in the 2050s under SSP245 and MPI-ESM1-2-HR and MPI-ESM1-2-LR in the 2080s under SPP585, which indicate an increase rather than a decrease in streamflow during the pre-monsoon. The predicted streamflow for the post-monsoon season, however, does not show a consistent pattern of change across time within the models under SSP245 and SSP585 scenarios. But there is a high probability, based on the mean ensembles of models projections, that streamflow will also decline during the post-monsoon in 2050s ($-$1.23 % to $-$0.22 %) and 2080s ($-$5.59 % to $-$2.83 %) under all scenarios. Similarly, the predicted decline for pre-monsoon and winter will be between $-$10.36 % and $-$6.12 % and $-$21.87 % and $-$21.52 % under SSP245 and between $-$10.0 % and $-$9.13 % and $-$21.87 % and $-$21.11 % under SSP585, respectively. With the exception of winter, when there are no significant differences in the projected streamflow, the decline will be slightly larger in the 2080s than in the 2050s in all scenarios. In addition, the results of the mean ensemble of the models indicate that the Sutlej River's flow will increase during the monsoon under both scenarios, from 9.70 % to 11.41 % in the 2050s and 11.64 % to 12.70 % in the 2080s.

Similarly, Fig. 19 lists the projected change in the mean annual streamflow for the Sutlej River in 2050s and 2080s with respect to the reference period (1979–2009) under different emission scenarios. Although the nature of the direction of change within models varies, the mean ensemble of the models reveals a persistent increasing pattern in the streamflow for all scenarios in 2050s and 2080s. The Sutlej River's annual streamflow will rise between 2050 and 2080 by 0.79 % to 1.43 % for SSP585 and 0.87 % to 1.10 % for SSP245, according to the mean ensemble of the models. The rise is expected to be higher in the 2080s as compared to 2050s under SSP585.

Figure 19

Predicted change in the mean annual streamflow of the Sutlej River with respect to the reference period (1979–2009) in the 2050s and 2080s under SSP245 and SSP585 scenarios for different GCMs.

[Figure omitted. See PDF]

5 Discussion

This study reveals an increase in the Sutlej River's mean annual and monsoonal streamflow in the 2050s and 2080s, which is in contrast to earlier studies (Singh et al., 2014; Ali et al., 2018) that reported a reduction based on long-term investigation of station data over historical era. The pattern of rainfall and temperature predicted by CMIP6 GCMs for the future periods under the SSP245 and SSP585 emission scenarios, in addition to the physical processes occurring within the basin, have contributed to this increase in the Sutlej River's streamflow. For instance, it is speculated that the projected increase in the mean streamflow during the monsoon season under both scenarios in the 2050s and 2080 for all models is related to the projected percentage increase in the rainfall amount over the catchment and the melting of glaciers brought on by the increased maximum and minimum temperatures. On the one hand, this increase in river streamflow and its propensity to raise the silt load may have an impact on both the capacity of reservoirs and the hydropower potential of hydroelectric facilities situated in the sub-basin and downstream of it. On the other hand, despite the predicted increase in rainfall throughout the pre-monsoon, post-monsoon, and winter seasons, the anticipated decrease in streamflow of the Sutlej River during pre-monsoon, post-monsoon, and winter may be explained by the projected rise in temperatures, which may have led to increased evaporation from the surface. Similar conclusions were reached by Adib and Harun (2022), who studied the Kurau River in Malaysia and predicted a drop in streamflow during the months of January, April, and October, despite receiving more rainfall. Moreover, during winter and post-monsoon, most of the precipitation in upper part of the catchment occurs in the form of snowfall, which has a minimal effect on the runoff generation in the catchment. Additionally, the large increase in monsoonal streamflow predicted during study periods is what led to the projected increase in the Sutlej River's mean annual flow. Predicted decreases in the Sutlej River streamflow over the pre-monsoon (April to June) and winter (December to March) seasons may have a significant impact on agriculture and hydropower generation downstream of the river, which is already struggling due to water shortages at this time of year. Ali et al. (2018) predicted that the hydroelectric production from the Nathpa Jhakri and Bhakra Nangal hydropower projects will decline during May to June in the future due to a projected decline in the streamflow of the Sutlej River.

The projected streamflow patterns for the Sutlej River under SSP245 and SSP585 scenarios in the 2050s and 2080s show similar tendencies, but with differing magnitudes, that have been found by past researchers using process-based hydrological models. For instance, Singh et al. (2015a) used the SWAT (Soil Water Assessment Tool) model, a semi-distributed hydrological model, to simulate streamflow for future periods using two CMIP3 GCMs (CGCM3 and HadCM3), and they discovered that the Sutlej River's mean annual streamflow would increase in the range of 0.6 % to 7.8 % for the future periods (2050s and 2080s). Similar to this, using the variable infiltration capacity (VIC) and SWAT models, respectively, Ali et al. (2018) and Shukla et al. (2021) estimated increases in the Sutlej River's mean annual streamflow for the 2050s and 2080s under RCP4.5 and RCP8.5. The study of Shukla et al. (2021) estimated that, under RCP4.5 and RCP8.5, the mean streamflow of the river would increase by 14 % and 21 % (at Rampur), respectively, in the 2080s. The previous studies observed that a substantially higher increase in projected streamflow may be attributable to the overestimation by the CMIP3 GCMs and CMIP5 GCMs of monsoonal precipitation over the Himalayan region (Choudhary et al., 2022; Sanjay et al., 2017; Gusain et al., 2020; Lalande et al., 2021). Similar to this, the results of Singh et al. (2015a), Ali et al. (2018), and Shukla et al. (2021) corroborated the expected decrease in streamflow during pre-monsoon and winter in addition to a rise during monsoon. This suggests that the RF model can accurately predict runoff and analyse the effects of climate change, while capturing the non-linearity of a hilly catchment.

6 Conclusion

This study compared the performance of the five machine learning models (GLM, GAM, MARSs, ANN, and RF) and one deep learning model (1D-CNN), which were further divided into linear (MARSs, ANN, and RF) and non-linear (ANN, 1D-CNN, and RF) models, in simulating rainfall–runoff responses over the hilly Sutlej River basin in order to determine the best model for predicting streamflow response to future climate change in the 2050s and 2080s under SSP245 and SSP585 using CMIP6 GCM data. The important findings of the study are summarised below.

• In general, non-linear models (ANN, 1D-CNN, and RF) outperformed linear models (GAM, GLM, and MARSs) in runoff prediction under all rainfall scenarios ($R_{\mathrm{0}}$, $R_{\mathrm{1}}$, $R_{\mathrm{2}}$, and $R_{\mathrm{3}}$). Among all the models, RF and 1D-CNN were identified as being the best models as per the model evaluation criteria. However, RF outperformed CNN in terms of error index (MAE and PBIAS), and as a result, it was used to investigate impact of future climate change on the Sutlej River pattern in the 2050s and 2080s under the SSP245 and SSP585 emission scenarios.

• The developed RF model slightly underperformed at the prediction of higher values of streamflow during training and testing. This implies that it is less effective at predicting flash floods that are caused by intense rainfall in the catchment. However, it was determined that the results produced by RF were comparable to process-based hydrological models for a long-term change study in streamflow pattern.

• Significant variations in the streamflow pattern were observed throughout the periods of months, seasons, and years and for the CMIP6 GCMs. The differences in the spatial resolution and parameterisation levels of CMIP6 GCMs, which caused a noticeable change in the projected amounts of temperature and precipitation during the study periods, may serve as an illustration of these variances in streamflow prediction. The Sutlej River's mean annual streamflow, based on the mean ensemble of models, is predicted to rise between the years 2050 and 2080 by 0.79 % to 1.43 % for SSP585 and by 0.87 % to 1.10 % for SSP245. Additionally, under both emission scenarios, streamflow will decrease during the pre- and post-monsoon ($-$1.23 % to $-$0.22 % and $-$5.59 % to $-$2.83 %) and during the winter ($-$21.87 % to $-$21.52 % and $-$21.87 % to $-$21.11 %) but increase during the monsoon (9.70 % to 11.41 % and 11.64 % to 12.70 %) in the 2050s and 2080s.

• The increase in the Sutlej River's streamflow (annual and monsoon) is due to both physical processes that occur within the basin and rainfall and temperature patterns that are predicted by CMIP6 GCMs for future time periods under the SSP245 and SSP585 emission scenarios. On the one hand, the projected rise in mean streamflow during the monsoon season is associated with both the projected percentage increase in rainfall over the catchment and the melting of glaciers brought on by the increasing maximum and minimum temperatures. On the other hand, the predicted increase in temperatures, which may have led to increased evaporation from the surface, may be used to explain the anticipated reduction in streamflow of the Sutlej River during pre-monsoon, post-monsoon, and winter.

• Additionally, the projected changes in the mean annual and seasonal streamflow of the river are consistent with earlier research done using process-based physical hydrological models. Thus, the outcomes of the overall study indicate that the RF model is efficient for simulating streamflow in the Himalayan catchment and that water availability during monsoon will rise as a result of an increase in catchment precipitation, which would eventually lead to an increased sediment load and affect hydropower generation. However, the predicted reduction in streamflow during pre-monsoon, post-monsoon, and winter will put stress on agriculture and hydropower generation downstream of the river, which is already struggling due to water shortages at this time of year. The administrators of local water resources and the government organisations in charge of maintaining reservoirs downriver may find these details on streamflow patterns to be of great use.

Code availability

The codes developed for this study can be made available to the readers on reasonable request to the corresponding author.

Data availability

The observed station data are confidential, and the authors do not have permission to share the data.

Author contributions

DS and SL conceptualised the problems, supervised the entire research activity from its inception to its completion, contributed to the data collection, processing, and interpretation, and wrote the research paper. MV and RS contributed to the development of the model, generation of figures, and analysis of data. PC and DC contributed to the data analysis and interpretation.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Hydrological response to climatic and cryospheric changes in high-mountain regions”. It is not associated with a conference.

Acknowledgements

The authors acknowledge the National Natural Science Foundation of China (NFSC; grant no. 42171129), for funding this research work, and the Bhakra Beas Management Board (BBMB), India, for the hydrometeorological data used in this study. We thank Pratibha Shrivastava, a fourth-year computer science student (B. Tech) at NIT Raipur, for her help in building the model.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant no. 42171129).

Review statement

This paper was edited by Yue-Ping Xu and reviewed by two anonymous referees.