Modelling hydrological responses under climate change using machine learning algorithms – semi-arid river basin of peninsular India

Hydrological models are promising tools to study various water resources engineering problems such as flood prediction and design, drought assessment, water quantity and quality assessment and hydrological responses under climate variability (Sood & Smakhtin 2015). The most important classification of hydrological models is empirical (data-driven), conceptual and physically based models (Aghakouchak & Habib 2010). Empirical models have evolved as data-driven using various regression (Zhang et al. 2020), machine learning algorithms such as Artificial Neural Network (ANN) (Hu et al. 2018), least square-Support Vector Regression (SVR) (Bharti et al. 2017), etc. Such empirical models are primarily based on observations without considering the various hydrological processes involved (Jaiswal et al. 2020), whereas the physically-based models work on the principles of conservation of mass, momentum and energy to consider various hydrological processes along with field measured parameters and other basin characteristics (Aghakouchak & Habib 2010). The conceptual models are based on the simplified mathematical conceptualization of a hydrological system by considering various components at catchment scales and use model parameters in simulating the discharges (Jaiswal et al. 2020). Physically-based models are complex and require fine spatio-temporal data, whereas purely data-driven empirical models are limited completely to observations. Compared to physical and empirical based hydrological models, the conceptual models have gained much interest in the research community due to their capability to include various hydrological processes (e.g. runoff, evapotranspiration, surface storage, etc.) and ease of implementation with available data (Ammann et al. 2020). In this context, various conceptual models have evolved in the past, such as HYMOD (Parra et al. 2018), HBV model (Vormoor et al. 2018), etc. Model parameters or calibration factors are most prominent in conceptual hydrological models affecting the accuracy of streamflow estimates (Arnold et al. 2012). In this context, auto-calibration approaches for hydrological models have progressed in recent years, for example, the Soil Water Assessment Tool – Calibration Uncertainty Programs (SWAT-CUP) (Ha et al. 2017), and the Shuffled Complex Evolution Algorithm (SCE-UA) which carries the limitation of multiple optimal parameters and corresponding accuracies (Pohlert et al. 2007). Calibration of hydrological model parameters using remotely sensed variables has also gained interest in recent years, based on satellite evapotranspiration data (Immerzeel & Droogers 2008; Rientjes et al. 2013). To this end, various catchment processes have been conceptualized in hydrological modelling by including model parameters whose values can be determined through calibration (Mianabadi et al. 2019). In this context, calibration factors were implemented for various processes such as closure of water balance (Jarsjö et al. 2008), radiation (Choudhury 1999) and vegetation dynamics (Donohue et al. 2007). Conceptual models developed by including such calibration factors on various processes have been used to study the hydrological responses under climate, land use and water changes as a function of various hydroclimatic variables and observed historical streamflows (e.g. Asokan et al. 2010).

Hydrological variables such as precipitation, evapotranspiration and streamflows have shown pronounced changes under climate signals and therefore the model parameters also should account for such changes (Srinivasa Raju & Nagesh Kumar 2018). Most of the conventional hydrological models apply such calibration parameters which were established based on observed historical streamflow data and other hydrological variables in the hydrological climate change impact assessment (e.g. Meenu et al. 2013). Such models integrate hydrological models and statistical downscaling models to study the hydrological responses based on the historically established model parameters (Hengade et al. 2018; Saraf & Regulwar 2018). Inclusion of time-invariant model parameters which are independent of the hydro-climatic variables and estimated solely based on the observed data are limited to capture the temporal variability of water-energy balance variables under climate change (Rehana et al. 2020a). The present study proposes to establish a relationship between the calibration factors and various hydroclimatic variables and uncalibrated runoff. Such an established functional form of calibration factors can be easily implementable for future scenarios under projected hydroclimatic variables. In this context, developing a calibration approach relating precipitation, evapotranspiration and uncalibrated streamflows can improve the model performance and understanding of hydrological processes (López et al. 2017). Therefore, estimating model parameters dynamically accounting for the expected hydrological variability under climate change can be promising. Furthermore, such approaches can be further extended with climate change projections to study the hydrological response of river basins.

The study proposed a modelling framework by relating the conceptual hydrological model parameter with various hydrological variables of precipitation, evapotranspiration and simulated runoff. Such established relationships will ease the estimation of model calibration factors for future scenarios with climate change projections obtained from statistical downscaling models. The major challenges in implementing such a framework are limited data and poor understanding of complex relationships between the calibration factors and various hydrological process-based variables. To overcome such limitations, the present study proposed an approach to model the parameters of conceptual models using various hydrological variables by adopting machine learning algorithms. The study emphasized the use of advanced machine learning algorithms for developing such complex relationships between model parameters and hydrological variables, which can form as a basis for studying the hydrological responses under climate variability. The proposed modelling framework has been demonstrated on a semi-arid river basin, Krishna River Basin (KRB), India.

Krishna River Basin (KRB) is the fifth largest river system in India and second largest river basin of peninsular India with a total catchment area of 258,948 km² (Figure 1). The KRB, 73°17′–81°9′E and 13°10′–19°22′N, has heavy precipitation over the Western Ghats with decreasing precipitation towards the upper and lower parts of the basin. The KRB originates in the Western Ghats of India and travels for about 1,400 km across the four states of Karnataka, Maharashtra, Andhra Pradesh and Telangana before joining the Bay of Bengal at Hamsaladeevi in Andhra Pradesh near the east coast of India. The major tributaries of the KRB include Ghataprabha, Malaprabha, Bhima, Tunga-Bhadra and Musi. About 44% of the KRB lies in Karnataka, with about 26, 15 and 15% of the basin in Maharashtra, Telangana and Andhra Pradesh respectively. The annual average precipitation over the basin is about 784 mm, with 90% of the precipitation occurring during the south-west monsoon from June to October (http://india-wris.nrsc.gov.in/wrpinfo/?title = Krishna). Most of the basin has a prevailing semi-arid climate and the states covering the basin (Maharashtra, Karnataka and Telangana) are majorly drought-prone (http://india-wris.nrsc.gov.in). The major water use of the KRB is for irrigation with about 61.9 billion cubic meters per year (BCM/year) along with domestic and industrial water uses as 1.6 and 3.2 BCM/year (Rooijen et al. 2009).

In this study, hydrological modelling of the KRB was performed using the distributed hydrological model, PCRaster-water flow module using POLFLOW code (Wit 2001; Jarsjö et al. 2008) (http://pcraster.geo.uu.nl/pcraster/4.2.0/documentation/pcraster_manual/sphinx). The KRB was delineated and a local drainage direction (ldd) map was created using the topographical data (Figure 2). The basin was discretized with a rectangular grid (1 × 1 km) and for each grid various hydrological processes can be provided as input. The Digital Elevation Map (DEM) data used was at a resolution of 30 arc-seconds (approximately 1 km) from the Global 30 Arc-Second Elevation (GTOPO30) data set provided by the US Geological Survey. Spatially distributed data of daily precipitation (P) and temperature (T) for the KRB were from the India Meteorological Department (IMD). The gridded daily precipitation data at a spatial resolution of 0.25 × 0.25° (Pai et al. 2014) and temperature at a resolution of 1 × 1° (Srivastava et al. 2009) were cropped to the basin for a common time period from 1951 to 2014. The temperature, which was available at 1 × 1° resolution, was interpolated to 0.25 × 0.25° using the inverse distance weighting method. The discharge data was obtained from Krishna & Godavari Basin Organization (KGBO), Central Water Commission (CWC), Hyderabad, Government of India (www.kgbo-cwc.ap.nic.in) for about 25 discharge locations for the period 1966–2015 (Figure 1(b)).

It can be noted that the study employed a conceptual catchment scale hydrological model on an annual time scale, neglecting the storage component of water in terms of groundwater in the water-balance equation. For long-term soil water storage changes the groundwater recharge influenced by geological features and human interaction water withdrawals are considered to be negligible (Wu et al. 2017). However, it can be noted that the storage component of water cannot be neglected at any temporal scale less than annual and the resulting phenomenon and the hydrological responses may vary due to the inclusion of groundwater storages being taken into account (Gunkel & Lange 2017).

The annual scale calibration factors (⁠⁠) were estimated and considered as predictand variable and uncalibrated runoff (P-AET) resulting from the PCRaster model, along with P and AET which were considered as predictor variables in training and testing of machine learning models. The present study used two machine learning models, Ensemble Regression (ER) and Random Forest (RF) to predict the model parameters.

Ensemble Regression (ER) methods are machine learning paradigms in which multiple methods which are often referred to as weak learners are trained to solve a problem and are combined to obtain better results (Friedman 2001). The ER model works on the hypothesis that a diverse set of models can make better predictions in comparison with an individual model. The ER models have gained interest in hydrological model assessments in recent years (Sajedi-Hosseini et al. 2018).

The input training dataset of ‘N’ points {X, Y} = ⁠, where is the set of predictors and as the observed predictand value at the ith timestep, will be considered. Initially, all predictors are given equal weighting coefficients and an initial model (⁠⁠) will be developed to predict values of the form ⁠.

In recent years, data mining algorithms such as Random Forest (RF) have been widely applied in various studies related to water resources engineering such as hydrological (Moore et al. 1991), soil moisture forecasting (Prasad et al. 2018) and water quality assessment (Bui et al. 2020). RFs are an ensemble of decision trees, which split the data space similar to the way humans and each node represents a predictor variable (Breiman 2001). The data splitting is generally done based on inequality condition, Gini index, which is a measure of the performance of the split and measures how diverse the data is until it reaches a terminal node (Breiman 2017). The RF tree size is determined based on the number of nodes, which serve as a basis to minimize the variance of each split. RFs can perform regression as well as classification using various decision trees based on bootstrap aggregation (bagging), which is a random data re-sampling methodology with replacement. The bootstrapping process was used to sample S data points with N variables by replacement. A decision tree was trained on each of these random samples. Bagging is the process of making overall predictions from average predictions of various bootstrapped samples trained on individual decision trees, making averaging predictions from each decision tree. The steps involved in designing the RF models are as follows:

• 1.

  Draw a bootstrap samples S* of size N from the training data

• 2.

  An RF tree will be developed to the bootstrapped data recursively for each terminal node of the tree.

• • • Selecting S variables at random from N variables

  • • Choosing the best split-point among the S

  • • Split the node into two daughter nodes

  • • Continuing the above steps until the size of the minimum node is reached.

• 3.

  The ensemble of trees can be estimated as

• 4.
  Prediction at a new point, x, can be made using the regression formulation as

  

  (13)

The RF model has been identified as a robust model due to its ability to capture the complex relationships between features and labels efficiently, accounting for the uncorrelated decision trees in the prediction.

The trained and tested machine learning algorithms will be further used to predict the calibration factors which can be used to estimate the calibrated runoff at the catchment scale. To quantify the performance of models and to understand the dependence between calibration factors and hydrological variables, Root Mean Square Error (RMSE), Absolute Percentage Bias (APB), Mean Absolute Error (MAE), Nash–Sutcliffe Efficiency (NSE) (Nash & Sutcliffe 1970), Kling-Gupta Efficiency (KGE) (Knoben et al. 2019) and Likelihood were selected as performance measures.

Statistical downscaling models are the state-of-the-art climate change projections prediction models relating to large-scale climate variables with surface hydrological variables (e.g. precipitation) using statistical methods (Wilby & Dawson 2013). The present study employed a multisite statistical downscaling model to predict the climate change projections of precipitation and temperature (Figure 3). The basic formulation of the downscaling model includes data preprocessing to remove systematic bias in the modelled and actual climate observations (bias correction), data reduction method (principal component analysis, PCA), predictand variable states estimation (K-means clustering), fitting algorithm to relate predictors and predictand states (CART), and transfer function (Support Vector Regression). National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) climate data from 1948 to the present at a resolution of 2.5 × 2.5° were considered as large-scale observed climate variables. The MIROC-ESM derived from the Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology, Japan, was considered as the GCM output. Spatial resolution mismatch between NCEP and GCM data was resolved by applying an inverse distance weighting interpolation. The GCM predictor data undergoes bias correction based on quantile mapping method by comparing the cumulative distribution functions (CDFs) of NCEP and GCM predictor data for the historical and future scenarios to remove the systematic bias associated with the climate outputs (Li et al. 2010). Both bias-corrected GCM and reanalysis climate data undergoes standardization, which involves subtraction of mean and division by standard deviation estimated with all the data points of the time series. After data processing, statistically significant climate variables (predictors) were considered to predict the precipitation over various grid points of the basin. The predictor variables which have shown significant correlation coefficients with precipitation were identified as surface air temperature, wind speed, humidity, etc., which were also found as potential predictors in other studies over India (Rehana & Mujumdar 2012; Salvi et al. 2013). NCEP/NCAR reanalysis data for surface air temperature, mean sea level pressure, specific humidity at 500 mb pressure level, zonal and meridional wind velocity at the surface level is extracted for the latitudes in the range of 12.5–20°N and the longitudes in the range of 72.5–82.5°E, surrounding the entire Krishna Basin.

The study applied a PCRaster distributed hydrological model for 348 grids encompassing the KRB at 0.25 × 0.25° resolution at the annual scale. The original water balance equation of (P-AET) was estimated for each grid from 1965 to 2014 to estimate the precipitation surplus. The precipitation surplus will be converted to runoff by using the area of each grid cell (Equation (7)). The basin outlet uncalibrated runoff will be estimated using the accuthresholdflux function of PCRaster using a local drainage network created based on DEM data. The observed and uncalibrated runoff along with precipitation and AET will be used to estimate the basin averaged model parameter, X_cal on an annual scale. The model parameters estimated based on the closure of water balance by comparing the observed and uncalibrated runoff will be considered as the predictand. The hydrological variables of P, AET and uncalibrated runoff (R_cal) are considered as the predictor variables in modelling of calibration factors using machine learning algorithms.

The calibration factors obtained based on the closure of water balance (Equation (10)) were considered as observed ones and we have used P, AET (Equation (5)) and uncalibrated runoff resulting from PCRaster (Equation (8)) as independent variables in the modelling of dynamic calibration factors using ER and RF models. The ER and RF models were trained with observed calibration factors obtained from Equation (10) as a dependent variable and P, AET and uncalibrated runoff obtained from PCRaster as independent variables for the 1966–2000 data. The trained ER and RF models have been tested with the independent data of P, AET and uncalibrated runoff for the period 2000–2014. The trained and tested ER and RF models were compared in terms of various performance measures, as given in Table 1. The error between the modelled and predicted values of calibration factors were found to be less for the RF model compared to ER. The RMSE and MAE values for the RF (ER) model was 0.04 (0.02) and 0.14 (0.09), respectively. Furthermore, the efficiency performance measures such as NSE and KGE have also shown that the RF model is robust to predict the calibration factors compared to the ER model with limited data points in training (34) and testing (15). Also, the RF model has been identified as a better fit compared to ER in terms of similarity measure (likelihood) and bias (APB), as given in Table 1 and Figure 4.

It can be noted that the calibration factors in Figure 4(a) and 4(b) are in the range of 0.9–1.5, which represents the multiplying factor to be applied on the AET. These are the range of values which have to be applied on the estimated AET based on Equation (5) to implement Equation (6) in the estimation of actual water availability in terms of runoff. The resulting runoff obtained after application of the calibration factor on AET will represent the calibrated runoff which is comparable with the observed runoff on the catchment scale. Figure 4(c) and 4(d) show the observed and simulated runoff values with ER and RF models for the training and testing periods. Due to the robustness of the RF model in modelling the calibration factors, the present study used the RF model to predict the calibration factors for the current and future scenarios. The trained and tested RF model has been used with projections of precipitation and temperature to estimate the projected runoff response over the basin. Figure 5 shows the comparison of NSE coefficients estimated between observed and NCEP and GCM simulations for the training period of 1951–1989 and the testing period of 1990–2005 of the statistical downscaling model for precipitation and temperature. The basin averaged NSE values estimated between observed and NCEP data for precipitation for the training and testing period were 0.65 and 0.44 respectively, whereas the NSE values estimated between observed and MIROC simulations of precipitation for the training and testing period were 0.57 and 0.32 respectively. Temperature predictions were more convincing compared to precipitation due to less variability in data. The basin averaged NSE values estimated between observed and NCEP data for temperature for the training and testing period were estimated as 0.98 and 0.93 respectively, whereas the NSE values estimated between observed and MIROC simulations of temperature for the training and testing period were 0.94 and 0.86 respectively. The observed basin averaged annual precipitation for the training and validation was estimated as 763.06 and 812.04 mm respectively, whereas the NCEP (MIROC) simulations for training and testing period were obtained as 663.68 (650.57) and 611.83 (648.26) mm, respectively. The observed mean annual temperature for the training (testing) period was noted as 26.26 (26.54) °C, with NCEP simulated as 26.25 (26.42) °C and MIROC simulated as 26.27 (26.20) °C. The future projections of precipitation and temperatures for the period of 2021–2080 were estimated with MIROC GCM under the Representative Concentration Pathway (RCP) 4.5 representing atmospheric radiation at 4.5 Wm^–2 at the end of 2100. Both precipitation and temperatures were predicted to increase over the KRB for the period of 2021–2080. The annual average precipitation increase was predicted to be about 13.12% with a temperature increase of 0.6 °C for the period 2061–2080 compared to the observed period of 1951–1989 with MIROC GCM outputs. Overall, the precipitation and temperature are predicted to increase over the KRB, affecting the water-energy balance variables (Rehana et al. 2020a). Such increase of precipitation and temperature projections based on GCM outputs was found to be comparable with the earlier research findings over the KRB based on Regional Circulation Model (RCM) based Coordinated Regional Downscaling Experiment (CORDEX) projections (e.g. Rehana et al. 2020a, 2020b). Their study noted an approximate increase of 2.19% in precipitation and an approximate increase of 1.29 °C in temperature for the period 2021–2040 and 2041–2060 respectively compared to the observed period of 1966–2003 with the RCP 4.5 scenario with various CORDEX model outputs.

The projections of precipitation and temperature have been used as input to the hydrological modelling framework as shown in Figure 2 to estimate the runoff for the current and future scenarios. The annual AET for observed period of 1951–1989 and 1990–2005 was estimated as 340.22 and 495.96 mm respectively. The estimated AET with MIROC GCM projections of precipitation and temperature for the periods 2020–2040, 2041–2060 and 2061–2080 are 580.63, 613.75 and 648.49 mm respectively (Table 2). Overall, the AET has been predicted to increase under an increase of temperature and the consequent increase of PET over the basin. The decrease of basin outlet runoff for the periods 2021–2040, 2041–2060 and 2061–2080 compared to the observed period of 1951–1989 was noted as 87.73 (15%), 86.17 (14.69%), 58.77 (10.09%) m³/s respectively over the KRB (Table 2). The increase of AET over the basin was found to be more pronounced compared to the increase of precipitation over the KRB using GCM climate change projections.

The present study proposed a modelling framework to study the hydrological response under climate change over the KRB. A conceptual hydrological model was adopted to study the runoff response under variability of P and AET. A time-variant calibration approach was developed to estimate the model parameters relating to precipitation, evapotranspiration and uncalibrated runoff using two machine learning algorithms, Ensemble Regression and Random Forest model. The Random Forest algorithm was identified as a robust method to model the calibration factors with limited data available for training and testing with precipitation, evapotranspiration and uncalibrated runoff. The developed model was further used to study the runoff response under climate change variability of precipitation and temperature. A statistical downscaling model based on K-means clustering, CART and SVR was developed to predict the future scenarios of precipitation and temperature. Both precipitation and temperature were predicted to increase by 13.12% and 0.6 °C respectively for the period 2061–2080 compared to the observed period of 1951–1989 with MIROC GCM outputs. The increase of temperature has predicted an increase of PET for the future scenarios. It can be noted that the AET estimates are defined by P and PET changes. The increase of temperature has resulted in an increase of PET and consequent increase of AET over KRB for the future scenarios. It can be noted that the projected increase of AET is more pronounced under an increase of PET compared to the increase of P, resulting in a decrease of runoff for the future scenarios over the KRB. The projected decrease of basin outlet runoff for the periods 2021–2040, 2041–2060 and 2061–2080 compared to the observed period 1951–1989 was noted as 15, 14.69 and 10.09% respectively. Overall, the basin outlet runoff was predicted to decrease (13.26%) for future scenarios under climate change due to an increase in temperature (0.6 °C) compared to the precipitation increase (13.12%) over the KRB. As the conceptual model applied is based on the basic hydrological variables of P and AET, the resulting runoff outlets were predicted to decrease under the projected increase of AET. Overall, severe water shortage was predicted under climate change due to the projected increase of AET compared to the precipitation over KRB.

The proposed modelling framework has considered one hydrological model with one particular GCM simulation under a particular concentration pathway scenario for demonstration purposes. It should be noted that hydrological response can be expected to vary with various hydrological models, calibration parameters, statistical downscaling models, various concentration path scenarios, and various machine learning algorithms. A detailed study may be evident under the possible range of uncertainties arising at each stage of the modelling framework, which can be a potential research problem.

The present study considered a calibration-free formulation of Budyko's hypothesis to estimate the annual AET on the catchment scale. The study intended to demonstrate how a hydrological calibration factor can be modelled based on various hydro-climatological variables and the same can be utilized for the climate change impact assessment. For this purpose, the catchment scale single calibration factor based conceptual hydrological model was adopted. However, the proposed methodology can be implemented with any such standard hydrological model, accounting for the number of calibration factors involved and the suitable data-driven model to be adopted. Furthermore, the present modelling framework can be extended by introducing multiple calibration factors accounting for the land use, soil, etc. (Jarsjö et al. 2008) in the hydrological modelling, which is a potential area of research. The main emphasis of the study was to develop a modelling framework to study the hydrological responses under climate change using a conceptual hydrological model with an integration of a calibration framework based on machine learning algorithms. The proposed modelling framework can be extendable for other river basins with sufficient data availability and validation by applying various machine learning algorithms. With the use of most readily available and operational hydrometeorological variables such as precipitation, temperature and runoff, the proposed model can be promising to study the catchment scale hydrological responses under climate change. However, validation of the modeled results is most critical in the implementation for other catchment studies. Specifically, the use of calibration-free Budyko formulation has to be validated with the water-balance or most dependable and state-of-the-art satellite based AET estimates. Furthermore, the study limited for annual water balances by neglecting the storage changes which can be further extended with the inclusion of appropriate storage components to study the water balances other than on annual scale. Due to the limitation over the availability of various meteorological variables in the estimation of PET, the present study adopted the Thornthwaite model, which can be extendable to the standard Penman model. Overall, the proposed methodology of inclusion of time-variant calibration factors in the conceptual hydrological model accounting for various hydrological variables is more reliable for catchment-scale water management under climate change.

The research work presented in the paper was funded by Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India through Start-up Grant for Young Scientists (YSS) Project no. YSS/2015/002111 to Dr S. Rehana. We sincerely acknowledge the India Meteorological Department (IMD) for providing the gridded observations of precipitation and temperatures and details of the data are given in the website at https://iridl.ldeo.columbia.edu/SOURCES/.IMD/index.html?Set-Language = en. The authors sincerely thank the editors and anonymous reviewers for reviewing the manuscript and providing insightful comments.

Data cannot be made publicly available; readers should contact the corresponding author for details.