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

: 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 streamﬂow (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 2 , 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 streamﬂow (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 streamﬂow 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 streamﬂow evaluated with this model can be used for the PMB’s future water management strategies. Thus, this study’s ﬁndings may be helpful in developing water management strategies and preventing the pessimistic effect of CC in the PMB.


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 2 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, 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.

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 2 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.

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 shows the resolution and source for each data set that has been used for model setup in SWAT.  [30]. For the PMB, bias-corrected daily precipitation and temperature data were acquired from (https://zenodo.org/record/3874046#.YOQWg0kzZPa).

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 shows the resolution and source for each data set that has been used for model setup in SWAT.  [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).

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. 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.

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: where: SW t = Final Soil water content in (mm). SW o = Initial Soil water content in (mm). t = Time in days. Rday = Precipitation in (mm). QSurf = Surface runoff in (mm). ETi = Evapotranspiration in (mm). WSeep = Percolation in (mm). 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.

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. 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: where: SWt = Final Soil water content in (mm). SWo = Initial Soil water content in (mm). t = Time in days. Rday = Precipitation in (mm). QSurf = Surface runoff in (mm).

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.

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.

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.

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).

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.

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. 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. The details about the scenarios used in this study are explained in Table S2. 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 Tables 3 and 4. Since it is hard to predict the precise Water 2023, 15, 1753 9 of 27 future LULC, therefore for the present study to analyse the climate change impact, LULC was assumed to be unaffected for the future.

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.

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.

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].

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].

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].

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 2 , 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 2 ): (R 2 ) 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 2 ). Low correlation is described with a value close to zero, whereas a high correlation is indicated by a value to close to one.
• 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.
• 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.
• 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.
where Qi = observed discharge, value of (R 2 ). Low correlation is described with a value close to zero, whereas a high correlation is indicated by a value to close to one.
 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.
 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.
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.
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.  value of (R 2 ). Low correlation is described with a value close to zero, whereas a high correlation is indicated by a value to close to one.
 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.
 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.
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.
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.

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 = 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.

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.

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 illustrates the kappa coefficients along with their strength of agreement [40].  Table 6 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.
Overall accuracy = Total number o f correctly classi f ied pixels Total number o f re f erence pixels × 100 = 80/100 × 100 = 80% Kappa coefficient (T) = (TS * TCS)−∑(column total×row total) 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  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.

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, it can be seen that the parameters that are the most sensitive are CN2 and ESCO.   [42][43][44][45]. Figure 4a,b illustrates the monthly calibrated SF. Table 8 explains the performance of the calibrated model.

Model Validation
The validation follows the same steps as the calibration. The monthly step model was  [42][43][44][45]. Figure 5a,b illustrate the monthly validated streamflow. Table 8 explains the performance of the validated model.  The results of the model performed well in the calibration and validation phases, as indicated by high R 2 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.

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).

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 reveals % increase in average annual precipitation for three GCMs under three scenarios for three future time horizons, 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.  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]. 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].

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 3 /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]. 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 3 /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%, 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 3 /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).

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 3 /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). 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 3 /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.

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

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 2 ) and water body (32.42 km 2 ) from 2010 had increased to (554.64 km 2 ) and (66.56 km 2 ), 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)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(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.

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 2 , 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.

Supplementary Materials:
The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/w15091753/s1, Table S1. FAO soil code details for the PMB, Table S2. Details about different SSPs adopted in the study, Table S3. Minimum, Maximum, and Fitted Value of the parameters (for the period 2003-2010), Table S4 Funding: Research and funding were provided under the project "Determining the potential of watercourses for the production of electric energy from micro, mini, and pico hydropower plants", from the University North, Croatia.

Institutional Review Board Statement:
This research is original and does not include any findings from other sources, with the exception of cited work. As a result, ethical standards are followed throughout the work.

Data Availability Statement:
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
This certifies that there are no academic or commercial conflicts of interest.