Understanding the Climate Change and Land Use Impact on Streamflow in the Present and Future under CMIP6 Climate Scenarios for the Parvara Mula Basin, India

Abstract

Understanding the likely impacts of climate change (CC) and Land Use Land Cover (LULC) on water resources (WR) is critical for a water basin’s mitigation. The present study intends to quantify the impact of (CC) and (LULC) on the streamflow (SF) of the Parvara Mula Basin (PMB) using SWAT. The SWAT model was calibrated and validated using the SWAT Calibration Uncertainty Program (SWAT-CUP) for the two time periods (2003–2007 and 2013–2016) and (2008–2010 and 2017–2018), respectively. To evaluate the model’s performance, statistical matrices such as R², NSE, PBIAS, and RSR were computed for both the calibrated and validated periods. For both these periods, the calibrated and validated results of the model were found to be very good. In this study, three bias-corrected CMIP6 GCMs (ACCESS-CM2, BCC-CSM2-MR, and CanESM5) under three scenarios (ssp245, ssp370, and ssp585) have been adopted by assuming no change in the existing LULC (2018). The results obtained from the SWAT simulation at the end of the century show that there will be an increase in streamflow (SF) by 44.75% to 53.72%, 45.80% to 77.31%, and 48.51% to 83.12% according to ACCESS-CM2, BCC-CSM2-MR, and CanESM5, respectively. A mean ensemble model was created to determine the net change in streamflow under different scenarios for different future time projections. The results obtained from the mean ensembled model also reveal an increase in the SF for the near future (2020–2040), mid future (2041–2070), and far future (2071–2100) to be 64.19%, 47.33%, and 70.59%, respectively. Finally, based on the obtained results, it was concluded that the CanESM5 model produces better results than the ACCESS-CM2 and BCC-CSM2-MR models. As a result, the streamflow evaluated with this model can be used for the PMB’s future water management strategies. Thus, this study’s findings may be helpful in developing water management strategies and preventing the pessimistic effect of CC in the PMB.

1. Introduction

Semi-arid areas are particularly vulnerable to human induced changes in the terrestrial ecosystem and environment [1]. These areas are habitat to 15% of the world’s population, who rely heavily on precipitation to supply their water needs for domestic, industrial, and agricultural uses [2]. Ref. [1] suggests that between 1990 and 2004, the entire semi-arid region rose by 7% globally. Food and water scarcity are forecast to occur in the dryland regions as a result of climate change (CC) and land use land cover (LULC) shifting above a critical level [3]. The global CO₂ concentration is rising at an alarming rate right now. As a result, since the start of the industrial revolution, the mean temperature has risen by 1 °C [4]. It is expected that global warming will continue to accelerate. Due to current climate change, a rise in the occurrence of drought events is predicted [3]. In particular, in the dry land regions, it will place enormous strain on greenery area cover, surface water properties, and agronomic production [5]. Arid or semi-arid regions cover over 35% of India’s land area, where 45% of the country’s population and livestock reside. Water supplies in these areas are dependent on rainfall brought in by the Indian summer monsoon, which provides almost 80% of unpolluted water resources. The climate and LULC of Central India have changed noticeably in recent decades. These changes have been identified as one of the primary causes of pre-monsoon precipitation failure due to air pollution, climate change factors, sea levels rising, and temperatures rising, etc., with many area seasons also changing slowly [6]. As a result, the occurrence of water scarcity rose by about half over the last two periods [7]. Groundwater levels have dropped due to decreased precipitation. There has been a rise in risky rainwater incidents by 10–30% and these are expected to continue [8]. These factors have had a significant impact on local rainwater availability. Thus, it is critical to examine the joint effect of climate change and LULC variations on the hydrological system in semi-arid river catchments for better water resource management [9,10]. The availability of new water is a source of apprehension for many people [11]. Climate factors such as changing rainfall patterns, increasing hotness, lowest temperature, possible evapotranspiration, and river flow have an impact on surface water resource planning and management [12]. Environment variation has a major effect on unpolluted water accessibility and variability as weather excesses such as heatwaves, droughts, and alterations in precipitation patterns become more common as an outcome of worldwide warming [13]. According to recent world reports, climate change is affecting the planet at an unprecedented rate. Rising temperatures, extreme weather events, and rising sea levels are just a few of the consequences of a changing climate. These changes are having a profound impact on land use, as the natural systems that support our planet are being pushed to their limits. One of the key drivers of changes in land use is agriculture. As temperatures rise, crops and livestock are struggling to survive in their traditional habitats, and farmers are being forced to adapt their practices to cope with changing conditions. This has led to a shift in the types of crops that are being grown, as well as the areas where farming is taking place. Another factor that is contributing to changes in land use is urbanization. As populations continue to grow, more and more people are moving into cities, which are expanding into previously rural areas. This has led to the destruction of natural habitats, as well as increased pollution and other environmental problems [10,11,13,14].

The effect of changing climate is particularly felt by developing nations whose economies depend heavily on rain feed agriculture [14]. As a huge part of India’s terrain is dry and semi-arid and depends on precipitation, the consequences of climate change directly affect the country’s economy. In order to incorporate various preventive measures, it is important to comprehend how climate change may affect future streamflow (SF) changes. Using a numerical hydrological model is the most accurate and thorough way to determine how changing climate may affect the SF of a river basin. Nevertheless, the choice of hydrological models depends on the individual and is influenced by factors such as data accessibility, site-specific issues, cost, and model accuracy [15]. According to the literature, the Soil Water Assessment Tool (SWAT) has been adopted by many researchers to determine the climate change impact on SF of a river [16,17,18]. Several studies have been completed using different climate scenarios and climate models to assess the combined effects of changing climate and LULC on the SF of various river basins across the globe. For instance, Ref. [19] used the Providing Regional Climates for Impacts Studies (PRECIS) weather and SWAT model to investigate the influence of forthcoming weather changes in the Upper Sind River Basin, India. They estimated that streamflow would rise by 16% and 93% in the mid and far future, respectively. Ref. [20] discovered that climate change had a greater impact on streamflow variability than LULC variation in the Brahmani-Baitarani basin in India. Ref. [21] used the SWAT model to evaluate the synergistic influence of weather patterns and LULC variation. They discovered that the RCP 4.5 scenario will result in a 50% rise in surface runoff, streamflow, and water supply in the Krishna River basin in India by the end of the century. Ref. [22] suggested the catchment’s LULC shift has increased streamflow in the Narmada River basin. Ref. [23] found that a 1.7% reduction in rainfall had a 7.8% impact on surface runoff in the Godavari basin. Ref. [24] assessed the predicted impact of global warming in the Betwa River basin by utilising the CMIP5 scenario. They discovered that streamflow is likely to rise by 4–29% and 12–48%, respectively, between 2040–2069 and 2070–2099. Most of the previous studies focused on single or at most two GCMs for predicting the future projections of precipitation, temperature, and SF in the basin. However, studies suggest that there is a significant level of uncertainty involved in a particular GCM prediction [25,26]. Hence ensemble models are advised in order to minimize such uncertainties [27]. In this study, in order to reduce uncertainty, a mean ensembled model is developed to determine the change in streamflow under different scenarios for different future time projections, i.e., near future (nf) (2020–2040), mid future (mf) (2041–2070), and far future (ff) (2071–2100).

To our knowledge, such research has not been performed to understand the impact of changing climate and LULC in the PMB. Over the past few decades, the PMB has undergone a number of changes, including LULC changes, an increase in the occurrence of extreme events, and a weakening of the summer monsoon rainfall [6,28,29]. Therefore, it is critical to understand the impact of climate change and LULC on the PMB under both the current and future climatic scenarios. Therefore, the main objective of this research is to create a strong foundation for long-term water resource forecasting and policy in the PMB by understanding how streamflow responds to climate change and LULC at a representative watershed size. The specific objectives are: (1) developing, calibrating, and validating the SWAT model for simulating the SF of the PMB; (2) investigating the effect of climate change on precipitation and temperature in the PMB based on multiple GCMs (ACCESS-CM2, BCC-CSM2-MR, and CanESM5) under different scenarios (ssp245, ssp370, and ssp585); (3) quantifying the potential impact of changing climate on the SF of the PMB; and (4) developing a mean ensemble model to determine the net change in SF under different scenarios for different future time projections, i.e., nf (2020–2040), mf (2041–2070), and ff (2071–2100). The findings of this study can be applied to long-term planning and policymaking, particularly in the management of water resources.

2. Study Area and Data Collection

2.1. Study Area

The PMB, which is situated in the Ahmednagar district of Maharashtra, India, was carefully chosen as the study area. The Pravara and the Mula are two independent Godavari tributaries. At an altitude of approximately 1295.4 metres above sea level, the Pravara River begins somewhat on the eastern side of the Sahyadri mountains, between Kulung and Ratangad. The Mula River is a main tributary of the Pravara River, which arises at a height of 1422 metres above MSL on the eastern edge of Harishchandragad.

The catchment size of the PMB is 5800 km² and is entirely in the state of Maharashtra, India. The basin serves the entire Ahmednagar district, which has a population of about 412,000 people. The PMB is located between 73°38′09.7312″ and 74°50′34.9810″ E longitudes and 19°02′23.1340″ and 19°45′18.9745″ N latitudes (Figure 1). It descends east from the Western Ghats and varies in elevation from 401 to 1596 metres above MSL. The precipitation increases from east to west, and it receives the majority of its precipitation during the southwest monsoon. It has an arid to semiarid climate, with average precipitation ranging between 500 and 1200 mm. In the basin, the monthly average temperature ranges from 14.38 °C to 38.60 °C. Figure 1 illustrates the location of the study area.

2.2. Data Collection

2.2.1. Hydro-Meteorological Data

Hydro-meteorological data such as rainfall, temperature, and discharge, along with the digital elevation model (DEM), soil map, slope map, and LULC map, form the key inputs for the SWAT model. Daily rainfall and temperature data were acquired for a 29-year period, from 1990 to 2018, from the India Meteorological Department (IMD). Daily discharge data for the gauging site Newasa, which is the outlet of the study area, for the same period of 29 years, from 1990 to 2018, has been obtained.

Table 1. Data Resolution and Source.

Data	Resolution	Source
DEM	30 m	https://earthexplorer.usgs.gov/ (accessed on 1 March 2023)
Soil Map	1000 m	https://data.apps.fao.org/map/catalog/srv/eng/catalog.search (accessed on 1 March 2023)
Slope Map	30 m	https://earthexplorer.usgs.gov/ (accessed on 1 March 2023)
LULC	30 m	https://earthexplorer.usgs.gov/ (accessed on 1 March 2023)
Rainfall	0.25 m	India Meteorological Department (IMD)
Temperature	1°	India Meteorological Department (IMD)
Discharge	-	Central water Commission (CWC)–Krishna Godavari Basin Organization

shows the resolution and source for each data set that has been used for model setup in SWAT.

2.2.2. GCM Data

CMIP6-GCM data were used for analysing the future SF for the period 2020–2100. For India, 13 GCMs are readily accessible at a resolution of 0.25 m × 0.25 m and are bias-corrected [30]. For the PMB, bias-corrected daily precipitation and temperature data were acquired from (https://zenodo.org/record/3874046#.YOQWg0kzZPa, accessed on 1 March 2023).

3. Methodology

The goal of this study is to assess the effect that changing climate has on the SF of the PMB. Since hydrology is a complex and uncertain process, models with various scenarios have been set up to assess the SWAT model’s appropriateness. The probable impact of climate change on hydrology is analysed in three different time horizons: nf (2020–2040), mf (2041–2070), and ff (2071–2100). Figure 2 illustrates the methodology used in this study. The following steps should be followed in order to achieve the set objectives:

• Collect and sort meteorological and hydrological data.
• Preparation of input maps such as DEM, soil, slope, and LULC.
• Watershed delineation using DEM involves fixing the basin outlet and setting a threshold area.
• HRU generation by overlaying the LULC, soil, and slope maps.
• Setting up the SWAT model with the input data.
• Evaluating the model and obtaining the output data.
• Calibrating the model using the SUFI-2 algorithm by finding the sensitive parameter.
• Checking the model’s performance for the calibration period; if okay, validating the model, or else adjusting the parameters again.
• Validation of the model under the fitted parameter values obtained during the calibration.
• Updating the SWAT model with future climatic data.
• Simulating the model with the future climatic data for impact assessment.

3.1. SWAT Model Setup

The PMB SWAT model simulates the SF at the rivers. The SWAT is a semi-distributed model which divides the catchment into a number of hydrological response units (HRUs) based on the similarity of the soil, LULC, and slope. It calculates surface runoff using the Soil Conservation Services (SCS) curve number method [31]. Lastly, it resolves the water balance equation (Equation (1)) to calculate the various hydrological processes in accordance with this formula:

$${\mathrm{SW}}_{\mathrm{t}}={\mathrm{SW}}_{\mathrm{o}}+{\sum }_{i=1}^t(\mathrm{Rday}-\mathrm{QSurf} - \mathrm{ETi}-\mathrm{WSeep} - \mathrm{Qgw})$$

(1)

where:

SW_t = Final Soil water content in (mm).

SW_o = Initial Soil water content in (mm).

t = Time in days.

$\mathrm{Rday}$ = Precipitation in (mm).

$\mathrm{QSurf}$ = Surface runoff in (mm).

$\mathrm{ETi}$ = Evapotranspiration in (mm).

$\mathrm{WSeep}$ = Percolation in (mm).

$\mathrm{Qgw}$ = Return flow in (mm).

SWAT uses meteorological data along with elevation, soil, and LULC as the input data. The model uses a soil map, a slope map, and a LULC map for the development of HRUs. The model was set up for two distinct time periods: 2000–2010 and 2011–2018. For the simulation of the years 2000–2010 and 2011–2018, respectively, we used the LULC map for the years 1990 and 2018 as a baseline. For both models, we selected the first two years, 2000–2010 and 2011–2018, as a warm-up period. The model is then executed to simulate SF.

3.1.1. Watershed Delineation

In this study, the Arc SWAT 2012 version was adopted for simulations. A DEM with a 30 m resolution was adopted to delineate the watershed. The DEM was used to study the drainage patterns and calculate sub basin features. The watershed delineation of the PMB has generated five sub basins, as shown in Figure 3a. In the Figure 3a, the purple colour line divides the basin into five sub basins, the blue colour line represents the river reach, and the black colour line represents the watershed boundary of the PMB.

3.1.2. Digital Elevation Model (DEM)

A DEM, which describes the elevation of every point in a given area at a definite spatial pixel size, was used to characterise terrain. A 30 m × 30 m resolution DEM was downloaded from the USGS earth explorer website. The DEM was used to represent the basin and examine the land surface regions, and it was also used to extract the stream network properties. Before extracting data, the DEM was projected to a coordinate system (WGS 1984 UTM Zone 43N). As seen in Figure 3b, the elevation values of the PMB range from 465 m to 1596 m.

3.1.3. Slope Map

The terrain’s slope has an impact on SF and infiltration. The flow increases with slope steepness. The slope, expressed as a percentage rise, ranges from 0% to 73.33%. Depending on the slope value, the slope map of the research area was classified into three groups, as shown in Figure 3c. DEM was used to create the study area’s slope map (30 m × 30 m resolution).

3.1.4. Soil Map

The global soil map is provided by the Food and Agriculture Organization (FAO) on a scale of 1:5,000,000. Figure 3d illustrates the soil map of the study area. Details of the hydrologic group and type of soil are listed in Table S1.

3.1.5. Land Use Land Cover

The satellite images of Landsat-8 obtained from the USGS Earth Explorer website were taken for LULC classification. The study area of the PMB was divided into five different classes: water bodies, built-up areas, vegetation land, agricultural land, and barren land. In ArcGIS, supervised classification of maximum likelihood was used to create the LULC map for the years 2010 and 2018, which are shown below in Figure 3e,f. Each LULC class area has been computed and listed in

Table 2. Representing areas for each class of LULC (2010) and LULC (2018).

LULC (2010)				LULC (2018)
Sr. No	Class	Area (km2)	% Area	Sr. No	Class	Area (km2)	% Area
1	Water Bodies	32.42	0.97	1	Water Bodies	66.56	1.19
2	Vegetation	125.19	2.25	2	Vegetation	16.66	0.30
3	Agricultural Land	1895.73	34	3	Agricultural Land	2046.46	36.70
4	Built-up Area	165.81	2.97	4	Built-up Area	554.64	9.95
5	Barren land	3355.97	60.19	5	Barren land	2890.78	51.85

.

3.1.6. Climate Data

For the analysis of future SF, CMIP6–GCM data were utilized. Out of the thirteen readily available GCMs, three were selected for the PMB based on the literature [32,33]. Furthermore, three scenarios were utilized, which were SSP245, SSP370, and SSP585. The GCMs used in this study are presented along with their sources in

Table 3. GCMs and their sources.

Sr. No	GCM	Source
1	ACCESS-CM2	Australian Community Climate and Earth System Simulator Coupled Model
2	BCC-CSM2-MR	Beijing Climate Centre Climate System Model
3	CanESM5	Canadian Earth System Model

. The details about the scenarios used in this study are explained in Table S2.

3.1.7. Hydrologic Response Unit (HRU) Generation

The HRU represents all surface features, including LULC, soil, and slope, as one unified group. If the surface properties are studied individually, the extent of the model’s complexity will increase. Therefore, it is desirable to merge all the attributes into a cohesive unit that has a lumped value. In the present study, the watershed was delineated with a threshold area of about 13,000 ha, which divided the basin into five sub basins. It was than overlain by LULC, soil, and slope, with each having a threshold value of 5%, respectively. As a result of which a total of 194 (HRUs) were generated. The LULC classification in terms of percentage area is illustrated in Table 3 and Table 4. Since it is hard to predict the precise future LULC, therefore for the present study to analyse the climate change impact, LULC was assumed to be unaffected for the future.

3.2. Model Simulation, Calibration, and Validation

3.2.1. SWAT Simulation

SWAT simulation enables us to prepare the input of the model and produce the result. The different hydrological elements of the basin can be simulated by running numerical simulation. The surface SF was generated as a result of the PMB simulation, which was carried out using ArcSWAT version 10.19.

3.2.2. SWAT Calibration and Uncertainty Analysis Program (SWAT-CUP)

SWAT calibration, validation, sensitivity, and uncertainty analysis were performed using SWAT-CUP. It is an open database program. It simplifies calibration and validation by reducing the time required for each process. It is assisted with graph-generating functionality. The SWAT-CUP is the bridge between SUFI-2 and SWAT.

3.2.3. Sensitivity Analysis

The change in model behaviour in relation to given input parameters is described as sensitivity analysis (SA). For model calibration and validation, seven parameters with initial ranges were chosen for SF simulation using the literature and the SWAT user guide. After considering the potential parameter ranges, SWAT-CUP utilised the sequential uncertainty fitting version two (SUFI-2) algorithm to obtain the most desirable parameters within 95 percent uncertainty ranges. The SA of the parameters can be achieved with the help of local sensitivity (LSA) or global sensitivity (GSA). In the present study, GSA was adopted to identify the most sensitive parameter for calibration and validation. The SA of the parameters was determined based on statistical measurements such as p values and t tests. The parameter’s sensitivity is directly proportional to the value of the t test and inversely proportional to the p value [34].

3.2.4. Model Calibration

Model calibration helps in determining which parameterized model is suitable for a certain set of discrete variables and lowers the uncertainty of forecasts. The parameters are adjusted within their acceptable ranges either by substitution, addition, or multiplication of the original values. The parameters were changed till the best simulation associated with observation was found [35].

3.2.5. Model Validation

Validation is the process of comparing calibrated parameters to objective datasets without making any changes to the values throughout the calibration process. The method of proposing that the model delivers suitably correct estimates at a specific location is known as “model validation” [36].

3.3. Model Performance Criteria

The results obtained from the calibration and validation processes are evaluated using different statistical matrices. In the present study, we evaluated the model’s performance using a number of statistical matrices such as R², NSE, PBIAS, RSR, P factor, and R factor, respectively [37,38]. Table 4 shows the range of performance ratings based on statistical metrics.

• Coefficient of determination (R²): (R²) calculates the correlation value between measured and estimated values by comparing the combined scattering of the measured and estimated series to the single scattering. Its value ranges between zero and one. The correlation between the observed and estimated SF can be understood with the value of (R²). Low correlation is described with a value close to zero, whereas a high correlation is indicated by a value to close to one.

$${\mathrm{R}}^2=\frac{{[{\sum }_{i=1}^n\left(Qi-\breve{Q}\right)·\left(Si-\breve{S}\right)]}^2}{\sqrt{{\sum }_{i=1}^n{\left(Qi-\breve{Q}\right)}^2·{\sum }_{i=1}^n{\left(Si-\breve{S}\right)}^2}}$$

(2)

• Nash–Sutcliffe efficiency (NSE): It is among the most extensively used hydrological model statistical measures. Its value varies from infinity to one, with one representing a perfect model. As the value approaches zero, the performance of the model degrades.

$$\mathrm{NSE}=1-\left[\frac{{\sum }_{i-=1}^n{(Qi-Si)}^2}{{\sum }_{i=1}^n{(Qi-\check{Q})}^2}\right]$$

(3)

• Percent bias (PBIAS): PBIAS shows the mean tendency of simulated outcome to be smaller or larger than the observed data. Its optimal value is zero. A positive PBIAS value means the model is an underprediction of the results, whereas a negative value suggests overprediction.

$$\mathrm{PBIAS}=\frac{{\sum }_{i=1}^n\left(Qi-Si\right)}{{\sum }_{i=1}^{n}Qi} \times 100$$

(4)

• Ratio of root-mean-square error to measured standard deviation (RSR): RSR is selected as a complimentary statistical metric to RMSE. The optimal value of RSR is zero. However, the higher the RSR, the lower the performance of the model.

$$\mathrm{RSR}=\frac{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{Qi}-\mathrm{Si})}^2/\mathrm{n}}}{\sqrt{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{Qi}-\breve{\mathrm{Q}})}^2/\mathrm{n}}}$$

(5)

where Qi = observed discharge, Q ̆ = Mean observed discharge,

Si = Simulated discharge, S ̆ = Mean simulated discharge.

The Soil and Water Assessment Tool (SWAT) uses various factors to simulate the hydrological processes that occur in a watershed, including the surface runoff generation factor (p) and the routing factor (r). Here are some common methods used to estimate these factors:

• Estimating p-factor: Curve Number (CN) Method: This method is based on the relationship between the antecedent moisture condition and the amount of rainfall that becomes runoff. The CN value is calculated using soil type, land use, and hydrologic soil group. Once the CN value is obtained, it is used to estimate the surface runoff using the SCS (Soil Conservation Service) equation.
• Estimating r-factor: Muskingum-Cunge Method: This method is based on the principle of flood routing and uses the Muskingum-Cunge equation to estimate the travel time and attenuation of flood waves. The equation requires input parameters such as channel length, channel slope, and channel roughness.

Table 4. Recommended performance of statistical metrics [39].

Performance	RSR	NSE	PBIAS
Very Good	0 ≤ RSR ≤ 0.5	0.75 < NSE ≤ 1	PBIAS < ± 10
Good	0.5 ≤ RSR ≤ 0.6	0.65 < NSE ≤ 0.75	± 10 ≤ PBIAS < ± 15
Satisfactory	0.6 ≤ RSR ≤ 0.7	0.50 < NSE ≤ 0.65	± 15 ≤ PBIAS < ± 25
Unsatisfactory	RSR > 0.7	NSE ≤ 0.50	PBIAS ≥ ± 25

Table 4. Recommended performance of statistical metrics [39].

Performance	RSR	NSE	PBIAS
Very Good	0 ≤ RSR ≤ 0.5	0.75 < NSE ≤ 1	PBIAS < ± 10
Good	0.5 ≤ RSR ≤ 0.6	0.65 < NSE ≤ 0.75	± 10 ≤ PBIAS < ± 15
Satisfactory	0.6 ≤ RSR ≤ 0.7	0.50 < NSE ≤ 0.65	± 15 ≤ PBIAS < ± 25
Unsatisfactory	RSR > 0.7	NSE ≤ 0.50	PBIAS ≥ ± 25

4. Results and Discussion

In this research, we are understanding the climate change and LULC variation effect on streamflow in current and future scenarios under CMIP6 climate scenarios for the Pravara Mula Basin, India. Therefore, 29 years of data from 1990 to 2018 were included as input data in the SWAT model with DEM, LULC, and slope. The main aim of this study is to find out the LULC impact on stream flow for sustainable development and planning purposes in the basin area. Currently, many climate challenges are impacting the LULC, agriculture, and water resource management. This model can be helpful to understand future information using the SWAT model. The details of the results are discussed in sections below.

4.1. LULC Accuracy Assessment Using Kappa Analysis

Kappa analysis is a distinct analysis of variance used in accuracy evaluations. Kappa analysis produces a Kappa Coefficient, which is a way of measuring agreement or accuracy. A Kappa coefficient of one indicates almost perfect agreement, while a zero indicates poor agreement.

Table 5. Rating criteria as per kappa analysis.

Sr. No	Kappa Coefficient	Agreement
1	<0	Poor
2	0–0.2	Slight
3	0.21–0.4	Fair
4	0.41–0.6	Moderate
5	0.61–0.8	Substantial
6	0.81–1.0	Almost Perfect

illustrates the kappa coefficients along with their strength of agreement [40].

Accuracy Assessment of LULC Map (2010) and LULC Map (2018)

Table 6. Accuracy assessment using kappa analysis for LULC map (2010) and LULC (2018), respectively.

LULC (2010)	Barren Land	Water Bodies	Built-Up Area	Vegetation	Agricultural Land	Total (User)
Barren Land	18	0	1	1	0	20
Water bodies	1	18	0	0	1	20
Built-Up Area	6	0	10	4	0	20
Vegetation	2	0	0	16	2	20
Agricultural Land	0	0	0	4	16	20
Total (producer)	27	18	11	25	19	100
LULU (2018)	Barren land	Water bodies	Built-Up Area	Vegetation	Agricultural Land	Total (user)
Barren Land	19	0	0	1	0	20
Water bodies	1	19	0	0	0	20
Built-Up Area	5	0	11	4	0	20
Vegetation	2	0	0	18	0	20
Agricultural Land	0	0	0	4	16	20
Total (producer)	27	19	11	27	16	100

reveals the accuracy assessment performed using kappa analysis for the years 2010 and 2018. In kappa analysis, the overall accuracy of the LULC map and kappa coefficient (T) are determined.

The year of 2010 LULC classes are Barren Land, Water Bodies, Built-Up Area, Vegetation, and Agricultural Land. The table is divided into two sections: the user section and the producer section. In the user section, the values represent the actual land cover, while in the producer section, the values represent the classified land cover. The total for each column is 20, representing the percentage of the land that falls into that category. For ex-ample, in the user section, 18% of the land was classified as Barren Land, while in the producer section, 27% of the land was classified as Barren Land. Kappa analysis is a statistical method used to assess the agreement between the actual and classified land cover. The value of kappa ranges from −1 to 1, with 1 indicating perfect agreement and 0 indicating random agreement.

$$\begin{gathered} \mathrm{Overall} \mathrm{accuracy}=\frac{Total number of correctly classified pixels}{Total number of reference pixels} \times 100 \\ =80/100\times 100=80\% \end{gathered}$$

(6)

$$\begin{gathered} \mathrm{Kappa} \mathrm{coefficient} \left(\mathrm{T}\right)=\frac{\left(TS\ast TCS\right)-\sum \left(column total\times row total\right)}{T{S}^2-\sum \left(column total\times row total\right)} \times 100 \\ =\frac{100\times 80-\left(27\times 20+18\times 20+11\times 20+25\times 20+19\times 20\right) }{100\times 100-\left(27\times 20+18\times 20+11\times 20+25\times 20+19\times 20\right) } \times 100 \\ =75\% \mathrm{or} 0.75 \end{gathered}$$

(7)

The Kappa Coefficient (T) obtained above using Equation (7) is 0.75, which is of substantial agreement as per the rating shown in Table 5. As a result, the LULC prepared for the year 2010 in this study is correct and can be used for further analysis.

In the LULC classes 2018, the values represent the actual land cover, while in the producer section, the values represent the classified land cover. The total for each column is 20, representing the percentage of the land that falls into that category. For example, in the user section, 19% of the land was classified as Barren Land, while in the producer section, 27% of the land was classified as Barren Land. Table 6 shows that in 2018, there was no Barren Land in the user section, while in the producer section, 19% of the land was classified as Barren Land. Additionally, the percentage of Built-Up Area increased from 10% in 2010 to 11% in 2018.

$$\begin{gathered} \mathrm{Overall} \mathrm{accuracy}=\frac{Total number of correctly classified pixels}{Total number of reference pixels} \times 100 \\ =83/100\times 100=83\% \end{gathered}$$

$$\begin{gathered} \mathrm{Kappa} \mathrm{coefficient} \left(\mathrm{T}\right)=\frac{\left(TS\times TCS\right)-\sum \left(column total\times row total\right)}{T{S}^2-\sum \left(column total\times row total\right)} \times 100 \\ =\frac{100\times 83-\left(27\times 20+19\times 20+11\times 20+27\times 20+16\times 20\right) }{100\times 100-\left(27\times 20+19\times 20+11\times 20+27\times 20+16\times 20\right) } \\ =78.75\% \mathrm{or} 0.7875 \end{gathered}$$

The Kappa Coefficient (T) obtained above using Equation (7) is 0.7875, which is of substantial agreement as per the rating shown in Table 5. As a result, the LULC prepared for the year 2018 in this study is correct and can be used for further analysis.

4.2. Sensitivity Analysis Using SUFI-2

To determine the optimum value for all the parameters, SUFI-2 algorithm was employed. In comparison with other methods, it requires a smaller number of simulations. For the calibration, it generates “n” numbers of combinations of parameters by using the Latin hypercube method. The process is runed until an optimal value is found. It describes the uncertainty with the help of the p and r factors [41]. The p factor is 95PPU, which is determined by 2.5% and 97.5% of an output variance. Whereas r factor is the ratio of mean thickness and standard deviation. Therefore, the SUFI-2 algorithm works on maximising these two factors to get the optimal parameter.

There are numerous stream parameters, but it is not necessary that each one of them be sensitive. As a result, prior to actually using it for calibration and validation, it was essential to recognize the utmost sensitive parameters. In this study, GSA is adopted for SA, which is performed by plotting dotty plots. If plots show a trend, it suggests high sensitivity, and if there is scattering, it means low sensitivity. SA is better understood with a t test, and the sensitivity of the parameter is directly proportional to the t test value [40]. To determine the specific ranges of the parameters, 1500 combinations of the parameters were calculated for each iteration. With the help of six iterations, we were able to determine the best range and optimal parameter. The validation process was completed by adopting these best ranges. There were in total seven parameters from the SWAT processes that were accounted for and put into simulations. These seven parameters are Curve number (R__CN2.mgt), Groundwater “revap” coefficient (V__GW_REVAP.gw), Soil evaporation compensation factor (V__ESCO.hru), Manning′s n value for the main channel (V__CH_N2.rte), Surface runoff lag time (R__SURLAG.bsn), Maximum canopy storage (R__CANMX.hru), and soil available water content (R__SOL_AWC(..).sol). Tables S3 and S4 show the fitted, minimum, and maximum values of the parameters. The parameter having the highest t stat and the lowest p value is the most sensitive. Hence, from

Table 7. Sensitivity Parameters (For the period 2003–2010) and (For the period 2013–2018).

Period (2003–2010)				Period (2013–2018)
Sr. No	Parameter	t Stat	p-Value	Sr. No	Parameter	t Stat	p-Value
1	R__CN2.mgt	27.81852	0.000000	1	R__CN2.mgt	17.81852	0.000000
2	V__GW_REVAP.gw	−1.11035	0.267024	2	V__GW_REVAP.gw	0.595071	0.551886
3	V__ESCO.hru	8.851582	0.000000	3	V__ESCO.hru	3.944408	0.000083
4	V__CH_N2.rte	−2.73454	0.006320	4	V__CH_N2.rte	0.416615	0.677019
5	R__SURLAG.bsn	0.437570	0.661760	5	R__SURLAG.bsn	0.168647	0.866096
6	R__CANMX.hru	0.249366	0.803111	6	R__CANMX.hru	0.615628	0.538233
7	R__SOL_AWC(..).sol	0.217374	0.827946	7	R__SOL_AWC(..).sol	0.052711	0.957968

, it can be seen that the parameters that are the most sensitive are CN2 and ESCO.

4.3. Model Calibration

Monthly model outputs were prepared for the period from 2000–2018 and were divided into two decades, i.e., 2000–2010 and 2011–2018. The model calibration was completed for both decades. For the first decade, 2000–2002 was set up as the warm-up period, whereas the calibration period was set as 2003–2007. Similarly, for the second decade, 2011–2012 was set as a warm-up period and 2013–2016 as the calibration period. The statistical metrics such as R², NSE, PBIAS, p-factor and r-factor were found for both of the calibrated periods. The obtained values of all the statistical metrics for the calibration period of 2003–2007 are (R² = 0.98), (NSE = 0.98), (PBIAS = 4.3), (RSR = 0.13), (p-factor = 0.72), and (r-factor = 0.19). Similarly, for the calibration period of 2013–2016, the values are (R² = 0.98), (NSE = 0.98), (PBIAS = 1.1), (RSR = 0.13), (p-factor = 0.71), and (r-factor = 0.44). From the results, it can be assessed that the model performed very well for both of the calibration periods. This result agrees with similar work conducted using SWAT [42,43,44,45]. Figure 4a,b illustrates the monthly calibrated SF.

Table 8. Performance analysis of the model during the calibration and validation period.

Sr. No	Parameter	Calibration (2003–2007)	Validation (2008–2010)	Calibration (2013–2016)	Validation (2017–2018)
1	R2	0.98	0.98	0.98	0.81
2	NSE	0.98	0.98	0.98	0.79
3	PBIAS	4.3	4.1	1.1	16
4	RSR	0.13	0.15	0.13	0.46
5	p-factor	0.72	0.64	0.71	0.5
6	r-factor	0.19	0.21	0.44	0.16

explains the performance of the calibrated model.

4.4. Model Validation

The validation follows the same steps as the calibration. The monthly step model was prepared for the period 2000–2018 and was divided into two decades, i.e., 2000–2010 and 2011–2018. The model validation was completed for both decades. For the first decade, the validation period was set at 2008–2010. Similarly, for the second decade, 2017–2018 was the validating period. The obtained values of all the statistical metrics for the validation period of 2008–2010 are (R² = 0.98), (NSE = 0.98), (PBIAS = 4.1), (RSR = 0.15), (p-factor = 0.64), and (r-factor = 0.21). Similarly, for the validation period of 2013–2016, the values are (R² = 0.81), (NSE = 0.79), (PBIAS = 16), (RSR = 0.46), (p-factor = 0.5), and (r-factor = 0.16). From the results, it can be concluded that the model performed very good for both of the validation periods. This result agrees with similar work conducted using SWAT [42,43,44,45]. Figure 5a,b illustrate the monthly validated streamflow. Table 8 explains the performance of the validated model.

The R² and NSE values of the model were consistently high across all time periods, indicating a good fit of the model to the observed data. However, the PBIAS values varied significantly across different phases, with the lowest value of 1.1 observed in the Calibration (2013–2016) period and the highest value of 16 in the Validation (2017–2018) period, indicating a decrease in model performance. The RSR values remained low across all periods, indicating that the model’s performance was consistent in both calibration and validation phases. However, the p-factor values decreased from 0.72 in the Calibration (2003–2007) period to 0.5 in the Validation (2017–2018) period, indicating reduced reliability of the model predictions. The r-factor values increased from 0.19 in the Calibration (2003–2007) period to 0.44 in the Calibration (2013–2016) period, indicating increased uncertainty associated with the model parameters. However, the r-factor value decreased to 0.16 in the Validation (2017–2018) period, indicating a reduction in the model’s uncertainty. The results of the model performed well in the calibration and validation phases, as indicated by high R² and NSE values and low RSR values. However, there were significant differences in the PBIAS values across different periods, indicating that the model’s performance varied over time. Additionally, the reliability of the model predictions decreased over time, as indicated by decreasing p-factor values and increasing r-factor values.

Figure 6a,b compares monthly observed and simulated SF hydrographs along with precipitation bar graphs, to help with the efficient understanding of the simulated results. These results clearly show that the simulated streamflow responded efficiently to the precipitation events. These figures also show that there is a similar trend between the monthly streamflow hydrographs that were observed and those that were simulated, although the model was unable to capture the entire shape of the streamflow hydrographs over the validation period (2017–2018). Hydrographs underpredict their peaks and recession limbs. The inadequacy in the matching of the peak flows could be associated with the limitations of rain gauges used to capture rainfall.

4.5. Projected Changes in Precipitation and Temperature for the Future Period (2020–2100)

The effect of climate change on precipitation and temperature were analysed using three GCMs (ACCESS-CM2, BCC-CSM2-MR, and CanESM5) from CMIP6 under three different scenarios (ssp245, ssp370, and ssp585) for nf (2020–2040), mf (2041–2070), and ff (2071–2100) in comparison to the baseline period (1990–2018). Figure 7 shows the typical pattern of monthly average precipitation projections for three GCMs under three scenarios for three future time horizons, respectively. Under all of the scenarios, it is expected that the projected precipitation for the three GCMs will increase. Precipitation is projected to increase more than in other months during the monsoon season (June to September). For ACCESS-CM2, average annual precipitation is likely to increase by 9.23%, 19.30%, and 25.80% under ssp245, ssp370, and sp585, respectively. Similar increases in average annual precipitation projection are seen for BCC-CSM2-MR (12.09%, 17.61%, and 20.12%) and CanESM5 (11.85%, 23.29%, and 30.83%) under ssp245, ssp370, and ssp585, respectively.

Table 9. % Increase in average annual precipitation, average annual temperature, and average annual streamflow for three GCMs under three different scenarios for nf (2020–2040), mf (2041–2070), and ff (2071–2100) respectively.

Average Annual Precipitation
GCM	Near Future			Mid Future			Far Future
	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585
ACCESS-CM2	8.19%	19.24%	21.52%	14.57%	22.27%	27.00%	14.95%	24.36%	30.00%
BCC-CSM2-MR	8.04%	18.71%	22.03%	12.09%	20.73%	26.57%	15.47%	23.41%	27.37%
CanESM5	9.80%	20.52%	22.32%	12.79%	24.37%	27.79%	16.67%	27.01%	33.30%
Average Annual Temperature
GCM	Near Future			Mid Future			Far Future
	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585
ACCESS-CM2	0.57%	5.32%	7.48%	0.79%	5.14%	11.41%	1.28%	6.56%	14.86%
BCC-CSM2-MR	1.08%	3.53%	5.71%	1.59%	4.26%	7.09%	2.84%	6.02%	9.56%
CanESM5	0.51%	2.35%	4.53%	0.70%	3.32%	6.97%	0.82%	4.89%	9.39%
Average Annual Streamflow
GCM	Near Future			Mid Future			Far Future
	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585	ssp245	ssp370	ssp585
ACCESS-CM2	34.29%	57.38%	32.73%	40.02%	33.90%	70.30%	51.28%	46.30%	71.92%
BCC-CSM2-MR	64.01%	33.70%	41.58%	67.13%	32.05%	30.41%	71.32%	68.01%	79.49%
CanESM5	69.88%	71.13%	77.62%	82.18%	72.05%	72.76%	64.62%	60.71%	71.34%
Ensemble Mean Model	56.04%	59.65%	68.15%	63.17%	45.39%	57.84%	72.21%	64.13%	80.74%

reveals % increase in average annual precipitation for three GCMs under three scenarios for three future time horizons, respectively.

Figure 8 shows the typical pattern of monthly average temperature projections for three GCMs under three scenarios for three future time horizons, respectively. As can be seen, under all of the scenarios for the three different time horizons, it is expected that the projected temperature for the three GCMs will increase. During the summer (March to May), temperatures are projected to rise higher than in other months. It is evident that the temperature will have a considerable impact, particularly in months with lower precipitation. For ACCESS-CM2, the average annual temperature is likely to increase by 4.46%, 5.78%, and 7.57% under ssp245, ssp370, and sp585, respectively. Similar increases in average annual temperature projection are seen for BCC-CSM2-MR (3.44%, 4.31%, and 6.14%) and CanESM5 (2.33%, 3.66%, and 5.03%) under ssp245, ssp370, and sp585, respectively. Table 9 reveals % increase in average annual temperature for three GCMs under three scenarios for three future time horizons, respectively. The findings unequivocally reveal that the temperature under ssp585 is comparatively higher than it was under ssp245 and ssp370, highlighting that ssp585 represents higher vulnerability to changing climate.

In order to attain a thorough understanding of how precipitation and temperature are projected to change over time in the PMB, the results of this study were evaluated by comparing them with a handful of other studies within the same and similar basins. The study’s findings were found to be consistent with those of [28], who predicted increases in rainfall and maximum temperature for future periods under three scenarios (A1B, A2, and B2) based on two models, CGCM3 and HadCM3. Similar findings were also identified in a number of other studies [32,46,47].

4.6. Projected Changes in Streamflow for the Future Period (2020–2100)

The calibrated SWAT model was employed to predict the probable impact of CC on the SF in three future time periods: nf (2020–2040), mf (2041–2070), and ff (2071–2100). Three GCMs (ACCESS-CM2, BCC-CSM2-MR, and CanESM5) from CMIP6 under three different scenarios (ssp245, ssp370, and ssp585) were employed.

4.6.1. Projected Changes in Streamflow under ACCESS-CM2

For predicting the future SF, we have kept the LULC of 2018 constant and adopted future precipitation and temperature from ACCESS-CM2 under the scenarios ssp245, ssp370, and ssp585. The simulated SF for the baseline period (2000–2018) is considered a reference. Figure 9 illustrates the change in monthly average streamflow (m³/sec) for three different time horizons in three scenarios under ACCESS-CM2. The annual average SF is expected to increase under all of the three scenarios ssp45, ssp370, and ssp585 by 50.84%, 44.75%, and 53.72%, respectively. For the ssp245 scenario, the runoff is expected to increase by 34.29% for nf (2020–2040), 40.02% for mf (2041–2070), and 51.28% for ff (2071–2100). In a similar manner, an expected increase in SF is observed for the ssp370 and ssp585 scenarios under different time horizons. For the ssp370 scenario, the SF is expected to increase by 57.38% for nf (2020–2040), 33.90% for mf (2041–2070), and 46.34% for ff (2071–2100). Furthermore, for ssp585, the SF is seen to show a probable increase of 32.73% for nf (2020–2040), 70.30% for mf (2041–2070), and 71.92% for ff (2071–2100). Table 9 reveals the % increase in average annual streamflow for three GCMs under three different scenarios for three different time horizons. This result is in agreement with similar work conducted on similar basins [48,49,50,51].

4.6.2. Projected Changes in Streamflow under BCC-CSM2-MR

The future projection of SF was calculated using precipitation and temperature from BCC-CSM2-MR under the scenarios ssp245, ssp370, and ssp585 and keeping the LULC of 2018 constant. Figure 9 illustrates changes in monthly average streamflow (m³/sec) for different future time horizons for three scenarios under BCC-CSM2-MR. The SF is expected to increase under all three scenarios, i.e., ssp245, ssp370, and ssp585 by 67.88%, 45.60%, and 77.31%, respectively. For the ssp245 scenario, the SF is expected to increase by 64.01% for nf (2020–2040), 67.13% for mf (2041–2070), and 71.32% for ff (2071–2100). In a similar manner, an expected increase in SF is observed for the ssp370 and ssp585 scenarios under different time horizons. For the ssp370 scenario, the SF is expected to increase by 33.70% for nf (2020–2040), 32.05% for mf (2041–2070), and 68.01% for ff (2071–2100). Furthermore, for ssp585, the SF is seen to show a probable increase of 41.58% for nf (2020–2040), 30.41% for mf (2041–2070), and 79.49% for ff (2071–2100).

4.6.3. Projected Changes in Streamflow under CanESM5

For predicting the future SF, we have kept the LULC of 2018 constant and adopted future precipitation and temperature from CanESM5 under the scenarios ssp245, ssp370, and ssp585. Figure 9 illustrates the change in monthly average streamflow (m³/sec) for different future time horizons for three scenarios under CanESM5. The SF is expected to increase under all three scenarios, i.e., ssp245, ssp370, and ssp585 by 48.51%, 67.96%, and 83.12%, respectively. For the ssp245 scenario, the SF is expected to increase by 69.88% for nf (2020–2040), 82.18% for mf (2041–2070), and 64.62% for ff (2071–2100). In a similar manner, an expected increase in SF is observed for the ssp370 and ssp585 scenarios under different time horizons. For the ssp370 scenario, the SF is expected to increase by 71.13% for nf (2020–2040), 72.05% for mf (2041–2070), and 60.71% for ff (2071–2100). Furthermore, for ssp585, the SF is seen to show a probable increase of 77.62% for nf (2020–2040), 72.76% for mf (2041–2070), and 71.34% for ff (2071–2100).

The results showed that the monthly average streamflow during future periods increased more under the CanESM5 model than it did under the ACCESS-CM2 and BCC-CSM2-MR models. This increase was due to the results of the precipitation projection, which showed a more incremental change under the CanESM5 model than under the ACCESS-CM2 and BCC-CSM2-MR models. Additionally, throughout the ff period (2071–2100), this increase was greater. This significant increase throughout the ff period could be attributed to the usage of a constant LULC (2018).

4.7. Mean Ensembled Model for Future Scenarios under Near Future (2020–2040), Mid Future (2041–2070), and Far Future (2071–2100) Time Horizons

In the present study, the mean ensemble model is developed by taking the mean of all three scenarios for the three different GCMs. With the help of such a model, we can determine the net effect of each scenario on the SF for the future projection in three different horizons: nf (2020–2040), mf (2041–2070), and ff (2071–2100), respectively.

Figure 10 illustrates the trend in the mean ensemble monthly average streamflow in (m³/sec) for three scenarios: nf (2020–2040), mf (2041–2070), and ff (2071–2100), respectively. In this model, we have tried to take the mean of all three GCMs under all the scenarios for the nf, mf, and ff periods. By doing this, an overall net result was obtained. According to [52], the ssp245 scenario faces the challenge of moderate mitigation and adaptation, the ssp370 scenario faces the challenge of high mitigation and adaptation, and the ssp585 scenario faces the challenge of high mitigation and low adaptation. For all the cases, increment in SF is projected for the nf period with reference to the baseline period. The increase in SF for the future projection period using different mean ensemble models for the nf period under three different scenarios such as ssp245, ssp370, and ssp585 are 56.04%, 59.65%, and 68.15%, respectively. The mean increase in SF is 64.19% for the nf period (2020–2040). For the mf period, the result of the mean ensemble model of three different scenarios reveals the increase in SF for the future projection period under the different scenarios ssp245, ssp370, and ssp585 are 63.17%, 45.39%, and 57.84%, respectively. Furthermore, the mean increase in SF is 47.33% for the mf (2040–2070). Similar increases in SF for the future projection period using different mean ensemble models for ff period under three different scenarios such as ssp245, ssp370, and ssp585 are 72.21%, 64.13%, and 80.74%, respectively. The mean increase in SF is 70.59% for the ff (2071–2100). Thus, with the help of this mean ensemble model, a combined increase in SF is obtained under multiple GCMs and various scenarios. Table 9 reveals the increase in average annual streamflow for the mean ensemble model under three different scenarios for nf (2020–2040), mf (2041–2070), and ff (2071–2100), respectively.

5. Discussion

We have examined the combined impact of climate change and LULC on the streamflow in the PMB. For this, the climatic data was divided into two periods: 2000–2010 and 2011–2018. In order to evaluate the impact of changing climate and LULC on the SF of the PMB, we selected baseline LULC (2018) for that time period because the LULC class varies gradually in the region. Between 2011 and 2018, we observed an increase in the catchment’s annual average SF runoff. We observed a direct correlation between changes in SF and changes in the catchment’s rainfall pattern. According to [53], by 2018 the climate affects SF more than LULC change. For example, land use such as the catchment’s built-up area (165.81 km²) and water body (32.42 km²) from 2010 had increased to (554.64 km²) and (66.56 km²), respectively.

There are a range of uncertainties when coupling climate models (RCM/GCM) with SWAT models to study the effects of climate change on SF, such as the choice of a GCM and plausible scenarios [54,55,56]. Despite these shortcomings, every attempt was taken in this study to minimize the uncertainties involved with hydrological simulation and climate prediction, in order to comprehend the likely impact of changing climate on SF. However, it is highly recommended to use multiple GCMs and ensemble models inclusive of all probable scenarios in order to examine a wide variety of climatic changes. Thus, in this study to predict the impact of changing climate on the SF for the future period, three GCMs (ACCESS-CM2, BCC-CSM2-MR, and CanESM5) from CMIP6 were employed under three different scenarios (ssp245, ssp370, and ssp585). Additionally, the findings of this study agree with those of previous related studies [57,58].

The study area is expected to have an overall increase in projected temperature and precipitation. Refs. [52,59,60] all revealed findings that were similar. The projected air temperature under ssp585 is higher compared to ssp245 and ssp370, highlighting that ssp585 represents higher vulnerability to changing climate. As a result, the hydrological cycle becomes active and results in significant rainfall [52]. Compared to other months, the monsoon season is predicted to have higher projected precipitation. The findings of [61,62] seem to support this conclusion. The annual average SF is expected to increase for all three GCMs under all three scenarios ssp45, ssp370, and ssp585. Table 9 reveals the details of % increase in SF under each scenario for all three future time horizons. Additionally, the findings of this study are consistent with other related studies [49,63,64].

According to a study by [63] in the Upper Godavari River, future precipitation and SF will increase. This is similar to what we found. However, the % increase in the annual average SF for future periods varies significantly. This is because different climatic models and reference periods were used, which could cause discrepancy. For the reference period (1985–2001), Ref. [63] considered two models under three scenarios. Hence, in order to reduce the uncertainty involved in using multiple GCMs and scenarios, in this study we have developed a mean ensemble model by taking the mean of each scenario of all GCMs under different time horizons, i.e., nf (2020–2040), mf (2041–2070), and ff (2071–2100). The results of the mean ensembled model as shown in Table 9 reveal a greater increase under the CanESM5 model than the ACCESS-CM2 and BCC-CSM2-MR models. This is due to the future precipitation projections of the CanESM5 model showing a greater increase than the ACCESS-CM2 and BCC-CSM2-MR models.

The findings of this study are centred on three GCMs under three scenarios (ssp245, ssp370, and ssp585). For a better understanding of the uncertainty involved, Ref. [64] suggest that more GCMs and emission scenarios should be included in hydrological studies. In order to simplify the study’s findings, it is suggested that additional GCMs and scenarios be included. Moreover, the SWAT model was calibrated and validated for a single site due to a lack of available data. Although the model proved effective in simulating the monthly SF, it is advisable to recalibrate the model on availability of data for multiple sites throughout the basin. Furthermore, LULC (2018) had been used to project future SF. Due to the PMB’s increasing population, changes in LULC are going to be significant in the future; however, this study did not take these changes into account.

We observed that, in comparison to the baseline period (2000–2018), the SF in the PMB increased all through the monsoon season (June–September) of the future period. Therefore, it is essential to preserve the surface water in the PMB in order to maintain the long-term sustainability of water management and meet the projected LULC and climate change. According to the many future climatic scenarios for the PMB, this study may be helpful to policymakers in developing informed laws and adaptive actions to reduce risk related to climate change.

6. Conclusions

In this study, hydrological modelling using SWAT was conducted to assess the impact of climate change on SF for the Parvara Mula Basin (PMB). Sensitivity analysis using the SUFI-2 algorithm in SWAT-CUP revealed that CN2 and ESCO are the most sensitive parameters. The SWAT model, coupled with SUFI-2, was established for the PMB, which was calibrated and validated using the observed SF data from the Newasa gauging site, which forms the basin’s outlet. The statistical matrices R², NSE, PBIAS, RSR, p-factor, and r-factor were employed to evaluate the performance of the model. Based on the results, the SWAT model performed very well in modelling the SF for the PMB. This model can be further extended to work for similar basins throughout the country.

The probable impact of climate change on SF for the future period (2020–2100) was analysed in three different time horizons: nf (2020–2040), mf (2041–2070), and ff (2071–2100). During this period, the LULC (2018) was kept constant, and precipitation and temperature projections were obtained from three GCMs which were ACCESS-CM2, BCC-CSM2-MR, and CanESM5 under three different scenarios which were ssp245, ssp370, and ssp585. Precipitation and temperature projections show increasing trends for all three climate scenarios (ssp245, sp370, and ssp585). Precipitation is projected to increase more than in other months during the monsoon season (June to September). The temperature rise is most significant under ssp585, confirming that ssp585 is the most vulnerable to climate change. All the models reveal a significant increase in SF for the coming years. According to ACCESS-CM2, the SF is expected to increase under all three scenarios by 44.75% to 53.72%. Similarly, the SF is expected to increase by 45.80% to 77.31% and 48.51% to 83.12%, according to BCC-CSM2-MR and CanESM5, respectively. The mean ensembled model determines the net effect of scenarios on the SF for all three future projections. The results obtained from the mean ensembled model also reveal an increase in the SF for all three future projections, i.e., nf (2020–2040), mf (2041–2070), and ff (2071–2100) to be 64.19%, 47.33%, and 70.59%, respectively. Moreover, the CanESM5 model produces more realistic results than the ACCESS-CM2 and BCC-CSM2-MR models, so the SF evaluated with this model can be utilised for water management strategies and planning efforts in the future. This study’s findings can be useful in developing water management practises and preventing the pessimistic effects of climate change in the PMB.