Assessing the impacts of climate change on drought-prone regions in Bhima sub-basin (India) using the Standard Precipitation Index

In Bhima sub-basin, the water sector is at high demand and in critical stress due to rapid urbanization. The past few decades have witnessed extreme events and seasonal shifts due to anthropogenic activity triggered climate change. Thus, to evaluate the variability of extreme events, assessing the historical and future trends of precipitation in climate change scenarios is vital for developing comprehensive mitigation and adaptation strategies. This paper examines the drought-prone regions by studying spatio-temporal variation of drought scenarios using the Standard Precipitation Index (SPI). The change factor method is used to downscale precipitation data from general circulation model (GCM) outputs under four Representative Concentration Pathway (RCP) scenarios to project future downscaled precipitation, to be input to examine the drought for 12 months. GCM and scenario uncertainty in climate change impact assessments are examined using box-whisker plots. Temporal variation for 12-month SPI shows significant changes over RCP scenarios. For the beginning of the period, 2021 precipitation is scanty for RCP 2.6 and 4.5 scenarios. Mild to moderate and severe drought events for the RCP 2.6 scenario are more predominant. Severe drought events under the RCP 6.0 scenario dominate over others. Lastly, the inconsistent pattern of drought events for RCP 8.5 is reported.


INTRODUCTION
Climate change in India will significantly impact agricultural production and water will be a major constraint in the future for food production (Aggarwal 2008). Developing countries like India are more vulnerable to the impacts of climate change and its implications on rain-fed agricultural areas have been affecting rural populations owing to the large dependency on agriculture for livelihoods. Agriculture primarily depends on monsoon, and changing precipitation patterns and intensities along with duration will probably influence the balance in water supply and demand mechanisms, thus affecting the hydrological cycle leading to frequent natural disasters (Gleick et al. 2010). Global warming is changing our climate, resulting in dire consequences. In the past few decades this has worsened as extreme events and seasonal shifts have bcome a reality. The consequences of global warming may not be uniform everywhere but its adverse effects will certainly affect life-supporting systems. According to Bates et al. (2010), the projected increase in greenhouse gas (GHG) emissions will intensify the water cycle and cause extreme events such as droughts and floods.
One of the most important sources of information for impact assessment studies is the output from general circulation models (GCMs) of the climate system (Robock et al. 1993). Owing to anthropogenic activities and variability in nature, it has been a challenge for researchers and policymakers to assess natural resources (more importantly water) with minimum uncertainty under perturbed climate conditions. Quite a number of studies in the recent past report the use of GCMs to assess the impacts of climate change at the regional level using statistical downscaling techniques (Woldemeskel et al. 2016;Sharma et al. 2018;Kang et al. 2019;Rehana 2019). GCMs can be used with higher confidence and less uncertainty using systematic evaluation of model performance (Nyunt et al. 2016).
The climate change impact assessment on the regional spatial scale (finer resolution) is hampered significantly by a considerable amount of uncertainty (Kundzewicz et al. 2018). The inter-modal uncertainty and inter-scenario uncertainty, referred to hereafter as GCM uncertainty and scenario uncertainty, respectively, have influenced the GCM outputs derived using statistical downscaling models. According to an Intergovernmental Panel on Climate Change IPCC report (2007), lack of scientific knowledge on climate processes and initial conditions results in GCM uncertainty. On the other hand, scenario uncertainty considers GHG emission scenarios under IPCC and causes unpredicted forecasting on the socio-economic and human intervention in policymaking. Thus, the uncertainty accumulates at various levels to propagate the large ranges at the regional level (Rowell 2006;Wilby & Harris 2006;Ghosh & Mujumdar 2007;Lindner et al. 2014;Fang et al. 2018).
The water sector is under high demand as well as under critical stress owing to rapid urbanization and industry needs. Large numbers of dams across the main river and tributaries of Krishna river (Bhima is a sub-basin of Krishna river basin) have changed the natural flow regime to such an extent that, almost zero flow conditions prevailed in the lower catchment area for a longer period of time (Venot et al. 2008). Hence, it is critical to further pursue a study on drought trends to examine the severity and its impact on food security and managing the water resource (Mahajan & Dodamani 2015). The mean and maximum temperature in the Krishna basin is rising significantly during all the seasons (Deshpande et al. 2016). Evidently, droughts in the Krishna basin have been occurring since 1901 with regular frequency and prolonged periods of ten to 15 years. However, in the 20th century, droughts lasted for a smaller duration of two to five years and were distributed widely and mainly concentrated in the upper Krishna and Bhima river basin . More recently, studies by Kulkarni et al. (2014) and Chanapathi & Thatikonda (2020) carried out work on Krishna basin using hydrological models to assess the impact of climate change on water resources. Approximately 64% of the population in Bhima sub-basin depends on agriculture as a major source of income and has been suffering due to frequent droughts in Maharashtra (upper Bhima subbasin). The Upper Bhima river basin is facing both periodic and prolonged water shortages due to exhaustive irrigation development (Garg et al. 2012). Udmale et al. (2014) highlighted several key challenges related to environmental degradation and ecosystem, putting this basin at potential risk regarding climate change. Surinaidu et al. (2013) reported a significant decline in groundwater (by À6 m), which for decades has endangered the basaltic aquifers and impacted agricultural activities. Developing countries basically lack an integrated approach to monitoring and effective evolution of agricultural projects. Thus, Biswas (1987) recommends the involvement of local authorities in the monitoring and evaluation process of watershed management to effectively mitigate the frequent droughts in the Upper Krishna basin in Maharashtra (Bhima sub-basin). Hence, robust and resilient strategies for managing the natural resources and climate change mitigation measures for the Bhima subbasin need to be sensitized to the local administration through policy documents. Watershed management policies set the guidelines and directives for the optimal utilization of resources. For instance, water policy guides the allocation priorities, namely, agriculture, industry, domestic, etc., as well as for awards, incentives and penalties towards conservation and discharge of water resources. Therefore, this paper on the impact assessment of climate change is vital in the framework of watershed management policies. rainfall and recorded average annual rainfall of 2,731.11 mm (Biggs et al. 2007). Lower Bhima had 678 mm of average annual rainfall meaning there was scanty rainfall in this region (Karnataka and Telangana). The temperature in the upper region (Bhimashanker) most of the times drops to 5-6°C in the lower reaches of Bhima (Sholapur district of Maharashtra) and sometimes rises up to 46°C. The Indian Meteorological Department (IMD) Pune reported a variation in annual mean temperature of 25°C in the upper region to 27.2°C in the lower region during 1969-2009. Literature examining the impacts of changing climate on water resources and watershed management for the Bhima subbasin reported that approximately 64% of the population depends on agriculture as a major source of income and has been suffering due to frequent droughts in Maharashtra (Upper Bhima sub-basin). The Upper Bhima river basin is facing both periodic and prolonged water shortages due to exhaustive irrigation development (Garg et al. 2012). Faster growing urbanization and the needs of stakeholders have disturbed the natural regime to hamper the rainfall distribution for the past decades causing an unbalanced supply and demand mechanism. Uneven distribution and uncertainty in rainfall pattern across this river basin has resulted in the stakeholders of this region depending on available groundwater resources, making the basin a waterscarce region. Greater parts of Bhima sub-basin covering Andhra Pradesh, Karnataka, Maharashtra and Telagana are facing massive rainfall deficits, drought scenarios and crop failures. Hence, examining the historical changes in the Bhima sub-basin region is, indeed, an important exercise and so historical trend analysis using the Mann-Kendall test is performed.
Historical trend analysis Mann (1945) proposed a non-parametric test based upon the application of Kendall's (Kendall 1975) test considering the randomness alongside the time for correlation and it is well-known to the research community as the Mann-Kendall test (MK test).
The mathematical equations for calculating Mann-Kendall statistics S, var and standardized test statistics Z are as follows: where, if S .0 then later observations in the time series tend to be larger than those that appear earlier in the time series, while the reverse is true if S ,0. The variance of S, (i.e., var) is given by where t varies over the set of tied ranks and f t is the number of times (i.e., frequency) that the rank t appears.
The MK test uses the following test statistic: where se ¼ the square root of var. If there is no monotonic trend (the null hypothesis), then for time series with more than ten elements, z ∼ N (0, 1), i.e., z has a standard normal distribution. The present study detects and assesses the significant trends in the spatial and temporal variations of historical data for the Bhima basin. For this purpose, the spatial variation of changes in rainfall IMD 0.25°Â0.25°gridded data for the period of 1969 to 2016 is an input for generating the annual mean rainfall over the entire Bhima basin. Figures 2 and 3 illustrate the annual mean rainfall and the maximum daily rainfall observed during 1969-2016, respectively. MK trend test 'tau' value ( Figure 4) explains the rainfall trend for the annual rainfall series of the IMD gridded data over the basin. The tau value of the MK test is quite similar and as good as the correlation coefficient with its value ranging from À1 to þ1. A positive (þ) tau value shows the increasing trend of the series of data (rainfall in this case) and its counterpart negative (À) tau value represents the decreasing trend. Figure 5 illustrates the average annual rainfall over the basin (with tau value). In a similar context, the trend in the annual maximum daily rainfall is analysed and the tau value of the gridded data with the trend is shown in Figures 6 and 7, respectively.
Upon examining the spatial plots and analysing the trends in annual rainfall over the basin, it is clear that there is neither a fixed increase nor a decrease of annual rainfall over the basin. However, the grids towards the Western Ghats region (Maharashtra state) lying in the Upper Bhima (K5) show increasing trends (positive tau values) in rainfall ( Figure 6). The annual    annual maximum rainfall reaches 2,300 mm and declines subsequently, such that, during 2005, the Bhima sub-basin experienced flash floods on account of extreme events of rainfall causing tangible and intangible effects on the people of this subbasin. However, the extreme rain events during 2005 are attributed to an increase in the fluctuations of the monsoon westerly winds, due to increased warming in the Arabian Sea (Kulkarni et al. 2020). The average slope of yearly time series data as per the linear equation is À0.2086 resulting in a 0.044% decline of annual rainfall for every decade. As the tau values are more predominant in Lower Bhima this may result in drought conditions in the regions of Lower Bhima.

Downscaling of precipitation
Statistical downscaling methods are generally used in climate change impact assessment studies for downscaling climate projections provided by the GCMs. GCMs are robust mathematical tools to estimate the probable changes in climate through the GHG concentrations to respond to future climate systems. Table 1 outlines the GCMs selected for downscaling the precipitation (rainfall in India is a major form of precipitation and thus, precipitation and rainfall in India are synonymous). The selection of these five GCMs is upon higher rankings having daily simulation for the future in four different RCPs on the ranking component considered by Raju et al. (2016) for the Godavari and Krishna basins of India (Bhima sub-basin is the subbasin of Krishna). A study by Raju et al. (2016) on 45 Coupled Model Intercomparison Project 5 (CMIP5)-based GCMs ranked the GCMs based on several metrics. However, the inclusion of all these GCMs into the current research is quite challenging owing to the non-availability of complete data and challenge that lies in handling the missing data (if any) for considering all or very large numbers of GCMs, such that many of these GCMs are not completely equipped with run data as provided by IPCC. Besides, the exclusion of some GCMs may even narrow the range of uncertainty represented by the remaining models. However, multi-model ensembles have now become a standard procedure to quantify and assess agreement and uncertainties among multiple GCMs in projecting future climate with a fair number of selections in GCMs (Wang et al. 2018).
A well-known statistical downscaling technique, the change factor (CF) method, is employed to downscale the precipitation data from GCM outputs under four Representative Concentration Pathways (RCP) scenarios published under the Fifth Assessment Report of the IPCC (van Vuuren et al. 2011). Among the wide ranging downscaling techniques, the CF method is a commonly used approach by the research community. CF is easy to implement and assess the impacts of changing climate at regional scale-related applications (Boé et al. 2007) and enabling representation of multiple GCMs and emission scenarios with the least minimal computing effort. For impact assessment studies involving complications, the CF technique has been adopted by the research community owing to its simplicity and flexibility of application with direct scaling of local data in line with changes suggested by the GCM simulations most relied upon. More information on downscaling climate variables using the CF method can be found in the papers by Boé et al. (2007) and Anandhi et al. (2011).
CF is broadly classified as: i) additive change factor (ACF) method and ii) multiplicative change factor (MCF) method. ACF evaluates the numerical difference in GCM variable obtained from present climate simulation and from a future climate scenario considered at the same GCM grid location as shown in Equation (4). This method is typically used for temperature, solar radiation, wind speed and relative humidity climate variables, except precipitation, since precipitation will not have a negative range.
The MCF method is the same as ACF except that the ratio between the future and current GCM simulations is calculated, instead of an arithmetical difference, as shown in Equation (5): where CF add and CF mul are additive and multiplicative change factors, respectively: GCM f is a GCM future climate scenario and GCM b is the value from the GCM baseline considered for a temporal domain. The next step is to obtain the local scaled future values (LSf add,i and LSf mul,i ) by relating the change factors obtained from Equations (1) and (2).
where LOb i are the observed values of climate variables considered at the ith time step on an individual meteorological station, or the averaged meteorological time series for a catchment for the designated temporal domain. LSf add,i and LSf mul,i are the values of future scenarios of climate variables obtained using additive and multiplicative formulation of the CF method. The methodology to estimate future scenarios using CF is best described pictorially in Figure 8. An output from a single GCM for a particular scenario has predominant uncertainties and should not account for predictions. Rather, to provide superior qualitative and quantitative climate change information it is essential to ascertain the uncertainties in the predicted results. This may be attained by the ensemble methods by minimizing the associated uncertainties in the multi-prediction approach. A multi-model ensemble average will help in addressing the underlying uncertainties by employing a set of climate models and scenarios (Knutti et al. 2010;Warner et al. 2010) by taking a simple arithmetic average. The IPCC has identified scenarios and, as per the IPCC's Fifth Assessment Report (AR5) these scenarios are called Representative Concentration Pathways (RCPs) with four pathways: RCP 8.5, RCP 6.0, RCP 4.5 and RCP 2.6. These RCP pathways have been adopted by the climate research community, terrestrial ecosystems and emission inventory experts, who have put forward comprehensive data sets with high spatial resolution until 2100 to designate different climate futures depending upon GHG emissions. The multi-model ensemble of precipitation projections for the four RCPs 2.6, 4.5, 6.0 and 8.5 validated the highest capability for precipitation downscaling (Homsi et al. 2019).
The mean annual rainfalls over the basin for the four RCP scenarios are shown in Figures 9-12. It is inferred from these spatial maps and histograms that the percentage (%) change in mean annual rainfall is positive and is increasing for most of the grids lying both in Upper Bhima and Lower Bhima region for all the four RCP scenarios. The percentage change in mean annual rainfall varies from 2% to 30% for the grids in the basin under RCP scenarios RCP 2.6, 4.5 and 6.0. Further, for RCP 8.5 scenario, this variation is slightly higher with the range extending to 35% and, moreover, during 2061-2100 almost all the grids in the basin under RCP 8.5 scenario depict extremity changes (30-35% variation) in the record of mean annual rainfall. However, this may be due to rising radiative forcing pathway RCP 8.5 scenarios leading to very high GHG emissions (IPCC 2013). A histogram, shown below each of the spatial maps, is used to summarize climate data. Histogram plot illustrates a graphic interpretation of numerical data by presenting the number of data points that fall within a quantified range of values. Red line in the histogram indicates mean percentage change.

Standard precipitation index (SPI)
For monitoring drought and guiding early warning, Steinemann (2003) defined drought indicator as a variable to identify and assess drought conditions. The levels of severity in climate can be categorized into 'mild', 'moderate', 'severe', or 'extreme'

drought. Generally used drought indicators include Standard Precipitation Index (SPI), Palmer Drought Severity Index (PDSI), Crop Moisture Index (CMI), Surface Water Supply Index (SWSI) and Reclamation Drought Index (RDI) (WMO & GWP 2016).
Most of the aforementioned indicators of drought need numerous input data such as rainfall, temperature, soil moisture, snowpack, reservoir storage, etc. However, SPI proposed by McKee et al. (1993) is used as a drought index using a single indicator 'precipitation' as an input variable. SPI has been a popular method owing to its 'ease of use'. SPI in comparison with other indices is an easy method for calculation and it is effective in giving early drought warning and drought damage control (Wang et al. 2015).
To assess the impact of drought in the past, the World Meteorological Organization (WMO) had recommended SPI to monitor and account for drought conditions (Hayes et al. 1999). This method makes use of the downscaled GCM data for forecasting future scenarios. It is the most popular among the other drought indices (Karabulut 2015) and also recognized widely for illustrating meteorological droughts (Hayes et al. 1999).
A recent study on monitoring of drought by Zhao et al. (2018) using Gravity Recovery and Climate Experiment (GRACE) data is compared and contrasted with SPI. Machine learning tools adopt the basics of SPI mechanism to predict standardized streamflow index for hydrological drought analysis (Shamshirband et al. 2020).
To determine the SPI values, the long-term rainfall data are fitted with a probability distribution. The current study intended to use the Gamma distribution to fit the long-term rainfall record similar to the methods adopted by Tsakiris et al. (2007) and Sönmez et al. (2005).
The standard procedure for SPI is as follows: (1) Gamma distribution is fitted to the time series of non-zero rainfall (X) for all the period 't' of interest (e.g., 3 months, 6 months, 9 months, 12 months, etc.). Avoid the overlapping between two time periods.
The Gamma probability distribution function (PDF) [g(X)] is given as: where α.0 and β.0 are shape and scale parameters, respectively. Fitting a Gamma distribution involves computation of parameters α and β. Further, shape and scale parameters α and β estimates are given by: where, and n is the number of periods of non-zero rainfall in the time series.
(2) Calculate the value of cumulative distribution function CDF [G(X)] conforming to each value of non-zero rainfall (X): G(X) can be computed using incomplete Gamma function.
(3) Estimate the zero rainfall (or, no-rainfall) probability (q) from the historical time series.
The value of CDF [H(X )] for an explicit rainfall (X) is then given by: (4) The value of SPI for the rainfall is computed as the value of a standard normal deviate corresponding to the value of CDF [H(X)]. McKee et al. (1993) used the categories as shown in Table 2 to define drought intensities resulting from the SPI.

Quantification of uncertainty
In order to address both GCM and scenario uncertainties an attempt is made in this section to interpret the outputs obtained from GCMs with four RCP scenarios in assessing the severity of climate change. It would be difficult to present the uncertainty study for each IMD grid in the basin and, thus, basin average value is considered for this research. Likely changes (based on the five GCMs ensemble average) in the aforementioned climate variables for future projected time series ranging from 2021 to 2100 and compared with that of their mean values over the historical period 1986-2005 are presented. Thus, for  Figure 13 shows box plots of the GCM and scenario uncertainty in climate change impact on projections of precipitation obtained over the evaluation period. Projections obtained from five GCMs and four RCP scenarios (ensemble average) are indicated as 'þ' (increase) and 'À' (decrease) sign with respect to average value of historical period and tabulated as maximum, mean and minimum percentage change of respective climate variables (Table 3). The historical value of the climate data is denoted by the horizontal line. The variability of the projection represented by the box in the box plot and median of the box plots indicate the deviation of the data.
Inferring from Table 3 and Figure 13, the mean percentage change for precipitation with the values 7.64%, 9.11%, 12.41% and 16.33% shows an increasing tendency along with the time series. The box plot for 2021-2040 and 2041-2060 shows the symmetric distribution of data, while the projections for 2061-2080 and 2081-2100 has skewness with the later one being more unsymmetrical and showing relatively more uncertainty than earlier ones.

Meteorological drought
Precipitation, being the major indicator for most of the themes on applications to link the availability of water, is the key input data and indicator of drought showing the different characteristics of a drought. The future projections of precipitation for the  period 2021-2100 involved CMIP5 scenarios along with RCPs as accepted in the AR5 released in 2014 by the IPCC. Time series of long-term records (2021-2100) of downscaled monthly precipitation data in the Bhima sub-basin is used as input to the SPI program. The salient features of the places under examination for analysis of 'dry' period are listed in Table 4 along with the grid location of places selected in Bhima sub-basin and with district/administrative boundaries' details in Figure 14.  Severity of future drought using SPI estimated from the downscaled GCM output The impact studies on meteorological changes, more importantly ascertaining drought, require quantitative assessment. Thus, SPI as drought index is used to project the future meteorological condition under different RCP scenarios to assess the impact on Bhima basin. To assess its ability to predict the drought conditions, monthly precipitation time series data for Bhima basin (average basin conditions) covering the period of 2021-2100 are considered. The SPI is capable of computing drought intensity over any specified interval, yet the most influential feature of this method is its inherent ability to simultaneously evaluate drought over a set of timescales, e.g., one, three, six, nine, 12, etc., popularly used by the research community. This research work considers SPI 12-month time scale using the algorithm presented in Gimbel et al. (2015) with modification using the 'precintcon' package to analyse the precipitation intensity, concentration and anomaly. SPI based on 12 months compares the precipitation data for the consecutive 12 months' recorded data within the same 12 consecutive months for all the available historical data. Since these particular timescales' data are the cumulative result of smaller periods, they may be above or below normal, making the longer SPI fall towards zero except for a typical dry or wet trend. Long-term precipitation patterns are reflected by these SPIs. These timescales are linked with the applications on stream flows, groundwater phenomena and the reservoir levels considered at the longer duration.
Temporal variation for 12-month SPI presented in Figure 15 shows significant changes over RCP scenarios. For the beginning of the period, 2021-2040 precipitation appears to be scanty for RCP 2.6 and RCP 4.5 scenarios. Mild to moderate and severe drought events for RCP 2.6 scenarios are more predominant as is evident from the summary of the report presented in Table 5. Extreme drought condition with 17.5% occurrence for the basin during 2021-2040 has more impact in the mid-region of the basin under G5, G6 and G7 grid point locations (Figure 16(a)). The upper region of Bhima (K5) under G1 and G2 grid locations has nearnormal conditions during this period. When we shift to the next time slot for 2041-2060, the occurrence of extreme drought events drops marginally. For the 2061-2080 time period, there is an increase in the extreme event of drought with 21.5% of the basin falling under the extreme drought category, and the grid locations G1 and G2 which were under near-normal conditions have been falling under this category as is evident from spatial plots. This event is continuous for 2081-2100 with 28% of the basin falling under extreme drought events. Near normal conditions are most likely to prevail during 2021-2100 for the RCP 4.5 scenario with the magnitude ranging from 26 to 35% occurrence for the Bhima sub-basin. The variations for different time intervals are shown in the spatial maps ( Figure 16(b)). The extreme drought events for RCP 4.5 tend to decline from the 2021 time period. However, during 2061-2080 there are zero occurrence extreme drought events (Table 5). Moderate and severe drought events have predominant impact during 2021-2040 in the lower region at grid locations G8, G9 and G10 and upper region at G1 and G2 grids during 2081-2100. Mild to moderate and severe drought conditions outweigh the other conditions under RCP 6.0. Most parts of upper Bhima are affected with these conditions (Figure 16(c)), with extreme drought events during 2027, 2047, 2058 and 2066 (as seen in the temporal variation plot, Figure 15). Extreme drought events occur with 20.5% intensity in the basin during 2021-2040 and mostly affect the grid locations G5, G6, G7 and G9. Spatial plots of RCP 8.5 scenarios shown in Figure 16(d) have varied occurrence of dry periods across the different time periods. The inconsistent pattern of occurrence of drought events over the time period shows 18% of the basin with extreme drought during 2041-2060 followed by 14%, 12% and 8.5% of its occurrence during 2021-2040, 2061-2080 and 2081-2100, respectively. It is probable that the emission scenarios of GHG considered for RCP 8.5 may likely resemble an inconsistent pattern of drought events.
Assessing the severity of future meteorological drought using the SPI estimated from the downscaled GCM output is yet another important objective of the research. To assess its ability to predict drought conditions, monthly precipitation time series data for Bhima basin (average basin conditions) covering the period of 2021-2100 are considered. Twelve-month SPI (i.e., SPI-12), the occurrence of mild to moderate drought events under RCP 2.6 and 6.0 scenarios is predominant for the basin, while severe drought events recorded for RCP 6.0 scenario dominates over other RCPs scenarios. However, the extreme drought events occur commonly and extensively under RCP 2.6 and 6.0 scenarios. The results under RCP 8.5 scenario show abnormal and extreme conditions in which GHG mitigation policies are not implemented and GHGs are emitted according to current trends. The current work on meteorological drought assessment will provide valuable information for managers and policymakers to investigate the impacts of climate change on drought conditions. The proposed work in this paper for developing the strategies and mitigation plan for regional assessment will serve as a decision-making tool for effective watershed