Appraisal of Climate Change and Its Impact on Water Resources of Pakistan: A Case Study of Mangla Watershed

: Water resources future, due to climate change. These results indicate that the study of climate change’s impact on the water resources under a suitable downscaling technique is imperative for proper planning and management of the water resources. the development of a hydrological model for the transboundary Mangla with the use of daily precipitation, temperature, and streamﬂow data (ii) future projection of precipitation and temperature from GCM data using three statistical downscaling techniques and their comparison, (iii) estimation of future climate change impacts on the transboundary Mangla streamﬂow. for of e ﬀ ects


Introduction
Most of the countries in South Asia are observing water stress due to global atmospheric changes. The rising population in urban areas, agriculture, and mismanagement of water resources and climate change has included Pakistan in the countries that are worst affected by climate change [1]. According to a report by the International Monetary Fund (IMF), Pakistan positioned third in the world among the countries prone to severe water scarcity [2]. The International Panel for Climate Change (IPCC) reported a 0.72 • C increase in the total temperature from 1951 to 2012. An expected rise of 1 • C to 3 • C till 2050 and 2 • C to 5 • C is likely to occur until 2100, based on the different greenhouse gases emission scenario. The severe temperature changes affect the land cover and ultimately change the streamflow patterns [3]. Rivers are providing more than fifty percent of the global water requirement [4]. However, river flows are associated with long-term fluctuations in rainfall and temperature, especially in areas where snowmelt is the principal component of the total runoff [5,6].
General circulation models (GCMs) are widely in use to study future climate change drifts. These models simulate climate variations based on possible projected greenhouse gas emission rates. The spatial resolution of almost all existing GCM models is around 150-300 km, and the spatial resolution of every single GCM alters when comparing with other GCMs [7]. To accurately apprehend the impact of climate change on water resources, the bias correction of the projected results from different climatic models is often performed [8,9]. Bias correction is beneficial, especially for hydrological modeling studies where streamflow directly connects with precipitation [10][11][12]. GCMs present large-scale forecasts for several climatic variables [3], but numerous climatic variables are not determined efficiently through the coarser resolution. To overcome this issue, different downscaling techniques are often applied to downscale the projected GCMs data to a fine resolution, but these techniques also provide systematic deviations [13][14][15]. Many statistical downscaling techniques have been introduced to eliminate systematic deviations [16,17]. Specific downscaling procedures establish a connection between regional-scale climatic components with the local scale climatic components. By studying the transfer of moisture at a regional scale and air blowing rate, the temperature and precipitation at a local scale can be predicted [18][19][20][21]. In statistical downscaling techniques, the spatial resolution of the GCMs is not considered, so the calculation of bias correction coefficient must be done effectively by using long period observed and historical climatic data [22,23]. Multiple studies deal with the drifts of hydro-climatic components of various basins in the upper Indus basin, using several statistical downscaling techniques, and a continuous rise in temperature was projected [24][25][26][27]. Nevertheless, the temperature is progressing at varying rates in various sub-basins. The history explains numerous rainfall drifts in various sub-basins present in the upper Indus basin [28].
To understand the hydrologic processes and interaction of water balance components under climate change scenarios in Mangla Watershed, the Soil and Water Assessment Tool (SWAT) model is used in this study. SWAT is proficient in modeling a single catchment or a system of hydrologically connected sub-catchments. The GIS-based interface model, ArcSWAT, defines the river network and the point of catchment outflow, and the distribution of sub-catchments and hydrological response units (HRU) [12,[29][30][31][32][33][34][35]. The calculation of the time of concentration by SWAT is done by adding overland flow time (time is taken for flow from the remotest point of the sub-basin to reach the water channel) and channel flow time (time taken from upstream channel to the watershed outlet) [36]. The projected climate is significantly affected by the selection of GCMs [35]. In this study, five global circulation models: The Australian Community Climate and Earth-System Simulator version 1 (ACCESS 1.0), Community Climate System Model (CCSM4), Hadley Centre Global Environment Model version 2 (HadGEM2), Max Planck Institute for Meteorology, Earth System Model Low Radiation Emission (MPI-ESM-LR), and Model for Interdisciplinary Research on Climate, Earth System Model (MirocESM), were selected based on their spatial resolution, and two different emission scenarios, RCP 4.5, and RCP 8.5 were chosen to reproduce future streamflow by applying the hydrological model SWAT (ArcSWAT-2012).
Several studies have defined precipitation and temperature variations as dominating factors that may affect water quantity and quality [37][38][39] and reported how climate change affects river flows [21,27,40]. Few studies have reported climate change effects on the average annual streamflow in the Mangla Watershed [32,41], however, the literature on the impact of future climate change on the annual and seasonal base flows and high flows in the watershed is missing. Moreover, investigation of the impact of future climate change on streamflow by considering the elevation band in watersheds is also missing. Previously, only biased corrections of the downscaled GCMs data had been performed without studying the different downscaling techniques. The current study has also used multiple downscaling techniques for projecting future water resources in the study area. The main objectives of this study are (i) the development of a hydrological model for the transboundary Mangla Dam watershed with the use of daily precipitation, temperature, and streamflow data (ii) future projection of precipitation and temperature from GCM data using three statistical downscaling techniques and their comparison, (iii) estimation of future climate change impacts on the transboundary Mangla Dam basin's streamflow. This study could be helpful to the water resource managers for the better administration of water resources, by modifying the management plans to mitigate the effects of climate change in the future.

Study Area
The transboundary Mangla Watershed ranges from 73 • 55' to 75 • 35' east of longitude and 33 • 25' to 34 • 40' north of latitude and is situated in the northeastern part of Pakistan and western part of the Himalaya. The geographical location of the watershed and considered weather station along with the land use and soil types in the watershed are presented in Figure 1. The Jhelum River is a main tributary of the Indus River and the Mangla Dam is built on this River. The maximum elevation in the Mangla watershed is about 5840 m, whereas the minimum elevation is 182 m above the mean sea level (a.m.s.l.). The total area of Mangla watershed is about 33,490 km 2 . Immense icecaps are present in the watershed, which provides large water volume to the Jhelum River by melting the ice, which eventually contributes to the Mangla Dam reservoir. The Mangla Dam was constructed in 1967 with an accommodation limit of 6.5 billion cubic meters of water. Almost 56% of the rear of the Mangla Watershed belongs to the Jhelum River, while 12% belongs to Poonch River, as shown in Figure 1. The total inflows into the Mangla dam are 1699.01 m 3 /s, whereas outflows are 566.34 m 3 /s. The hydro-power generation capacity of the Mangla dam is about 1310 MW. The mean population density in the watershed varies from 350 to 1000 people per square kilometer in the mountainous and cultivated regions, respectively. The average Tmax of Mangla Watershed ranges between 10.6-31.6 • C, and the average minimum temperature ranges from 4-25 • C. Soil data were downloaded from the website [42] land use data, and DEM data were downloaded from [43,44].

Climate Data
Daily observed data of climate variables such as Tmax, Tmin, and precipitation were obtained from the statistics files of the Pakistan Meteorological Department (PMD) and Water and Power Development Authority (WAPDA). Daily rainfall data for years 1981-2010 were provided by Pakistan Meteorology Department to quantify the outflow of the basin. Figure 2 shows the average annual temperature trends from 1981 to 2010 at Mangla Watershed. The maximum temperature increased at a rate of 0.373 °C per decade (Figure 2a), and the rate of rising of minimum temperature was observed as 0.139 °C per decade (Figure 2b).

Climate Data
Daily observed data of climate variables such as T max , T min , and precipitation were obtained from the statistics files of the Pakistan Meteorological Department (PMD) and Water and Power Development Authority (WAPDA). Daily rainfall data for years 1981-2010 were provided by Pakistan Meteorology Department to quantify the outflow of the basin. Figure 2 shows the average annual temperature trends from 1981 to 2010 at Mangla Watershed. The maximum temperature increased at a rate of 0.373 • C per decade (Figure 2a), and the rate of rising of minimum temperature was observed as 0.139 • C per decade ( Figure 2b).

Climate Data
Daily observed data of climate variables such as Tmax, Tmin, and precipitation were obtained from the statistics files of the Pakistan Meteorological Department (PMD) and Water and Power Development Authority (WAPDA). Daily rainfall data for years 1981-2010 were provided by Pakistan Meteorology Department to quantify the outflow of the basin. Figure 2 shows the average annual temperature trends from 1981 to 2010 at Mangla Watershed. The maximum temperature increased at a rate of 0.373 °C per decade (Figure 2a), and the rate of rising of minimum temperature was observed as 0.139 °C per decade (Figure 2b).

Global Climatic Models (GCMs) Data
In this study, five global circulation models; ACCESS, CCSM4, HadGEM2, ESMLR, and MirocESM were selected based on their spatial resolution [45]. Two representative concentration pathways (RCPs), Atmosphere 2020, 11, 1071 5 of 32 RCP 4.5, and RCP 8.5 were considered for these five GCMs. The data were downloaded from [46]. The details about these GCMs can be found from the user guide of CMIP5 [7].

Landuse and Soil Data
The soil data were downloaded from the Food and Agriculture Organization (FAO) website [42]. The downloaded harmonized soil 30 arc-second raster database covers the entire globe and displays the composition and characterization of different soil parameters, such as soil texture, ion exchange capacity, water holding capacity, soil depth, soil pH, soil temperature, lime and gypsum contents, granulometry, etc. The information about land use land cover (LULC) was downloaded from the United States Geological Survey (USGS) website [44]. The land-use data is a raster dataset, with a fine spatial resolution of 30 m by 30 m and covers a detailed 27 land cover and vegetation classification. The land use data provide detailed information about several types of agricultural lands, forests, vegetation, and grasslands present in any specific part of the globe. They also provide details about the snow-covered area, water bodies, artificial or urban areas, and barren lands. The hydrological model SWAT was parameterized using remotely sensed satellite datasets on soil, land cover, and topography. Based on the USGS soil classification system, the Mangla watershed have more than 45 percent of agricultural lands. Table 1 revealed that the most prominent land-use existing in the watershed is croplands, either artificially or irrigated covering, almost 45.38 percent of the area. The second widely spread class is grasslands, or green lands covered by needle leaves plants and covers 18.85 percent of the watershed area. The permanent snow-covered region covers approximately 3.91 percent of the Mangla watershed. Almost 1300 km 2 of the watershed is under permanent ice-caps and snow reserves. Details of soil data are presented in Table 2. The most prominent soil type in the watershed is loamy soil, which covered up to 71.04% of Mangla Watershed. The second classification represents the clayey soil, covers 24% of the total area, and other soils or water reservoirs surround the remaining 5% of the entire area. The GleyicSolonchak soil under group C holds 47.08% area of the watershed. Soil group C denotes those soils having low infiltration speeds when they remain fully saturated and uniform, with textures from medium-fine to moderately thick. These soils allow a reasonable water transportation rate. The principal soil collection comprises the Gleyic Solonchak group, which includes almost 47.087% of the watershed. The other groups hold Calcaric Phaeozems that cover about 23.429% of the area. The Mollic Planosols, Haplic Chernozems, Haplic Solonetz, Calcic Chernozems, Gelic Regosols, Luvic Chernozems, Lithic Leptosols, and Dystric Cambisols occupy 20.694%, 1.889%, 1.616%, 1.491%, 1.163%, 1.162%, 0.702%, 0.418%, and 0.342% of the total area of the watershed. Soil characteristics such as soil texture, soil bulk density, soil composition, soil electric conductivity, and soil possible water retention were collected from the FAO website. Wetland forests (WETF) are tree type crops [47]. The SWAT model has a built-in function to handle WETF by using several crop parameters, such as the ratio of erosion from land cropped to the clean tilled fallow continuous (USLE C factor), temperature responses, leaf area development, energy biomass conversion, light interception, stomatal conductance, canopy height, and root depth, plant nutrient content, and harvest. Detailed information on these parameters is given in [47].

Statistical Downscaling
The large-scale gridded GCMs data were downscaled to predict the climatic conditions at a local scale. Statistical downscaling can be performed by setting mathematical relations between large-scale GCMs climatic variables and the local climatic variable. Statistical downscaling is much easy to interpret than dynamical downscaling. However, statistical downscaling depends massively on historical climatic data and gauged data. In this study, three different types of downscaling techniques were used.

1.
Change factor (additive for temperature and multiplicative for precipitation) 2.
Linear scaling (additive for temperature and multiplicative for precipitation) 3.
Distribution mapping (additive for temperature and multiplicative for precipitation) In change factor downscaling, the bias correction for the temperature (maximum and minimum) was done by adding the difference of monthly changes among the future and base years GCM in Atmosphere 2020, 11, 1071 7 of 32 the observed temperature data. The bias correction coefficient was calculated for precipitation by dividing the monthly mean of GCM future precipitation data by GCM baseline precipitation data. The calculated bias correction coefficient was multiplied by observed precipitation data to obtain the modified/adjusted/bias-corrected precipitation [48]. The bias correction coefficient for precipitation in linear scaling downscaling was calculated by dividing long term monthly mean of gauged precipitation data with that of mean monthly GCMs, then this bias correction coefficient was multiplied with daily GCMs raw data to get corrected GCMs data [49]. Similarly, for the time-series data of temperature, an additive quotient was calculated by subtracting the long-term monthly mean of the GCMs simulated data from the long term monthly mean of observed data. The calculated bias correction coefficient was added in the daily GCMs raw data to obtain corrected GCMs data [49]. Distribution mapping (DM) is a commonly used bias-correction technique employed in various climate change studies [12,33,34,50]. Its basic working principle is to reproduce a transfer function based on the monthly mean value of GCMs data to the gauged data. The Gamma transfer function and Gaussian transfer function for precipitation and temperature data need to be reproduced, respectively, for correcting the bias from the downscaled future GCMs data [49]. Table 3 describes the five selected global climatic models used in this study. The GCMs were selected based on their spatial resolution. The projected climatic parameters vary extensively for different GCMs [35] and multiple GCMs are recommended for future projections, as single GCMs may cause uncertainty in the projected results. The data from five GCMs were downscaled by using three downscaling techniques, change factor downscaling, linear scaling downscaling, and distribution mapping downscaling. To investigate the suitability of the downscaling techniques, the downscaled data from all three downscaling techniques were compared with observed data using different statistical parameters, such as standard deviation and correlation coefficient. Figure 3 explains that there is a strong correlation between all three downscaling techniques and gauge data, i.e., above 0.56. In the case of the distribution mapping downscaling technique, the results were better than linear scaling and change factor downscaling techniques; the correlation of distribution mapping downscaling with observed data is above 90% for the GCM ESMLR in case of minimum temperature. Based on the correlation coefficient, standard deviation and root mean square deviation, as shown in Figure 3, it was concluded that the GCMs, ESMLR, and MirocESM has a maximum value of correlation coefficient (CC) as 0.83 and 0.80 respectively and minimum value of root mean square deviation (RMSD) 4.41 and 4.61 respectively for Tmax. ESMLR and MirocESM also have a maximum value of CC as 0.91 and 0.89 and a minimum value of RMSD as 3.14 and 3.42 respectively for Tmin. After the analysis of precipitation data, the CC of ESMLR and MirocESM also resulted in maximum reading as 0.81 and 0.74 under the distribution mapping downscaling technique. The GCMs, MirocESM, and ESMLR were found to be more accurate than other GCMs, and distribution mapping downscaling was found to be a more suitable downscaling technique for the transboundary of Mangla watershed to study climatic parameters such as maximum and minimum temperature and precipitation. The linear scaling downscaling technique is found to be less accurate than change factor downscaling [51]. The findings of this research consistent with previous studies such as [51] predicted that linear scaling is less suitable than change factor, whereas [52] studied that the distribution mapping performed better than other downscaling techniques. CMHyd software was used for linear scaling and distribution mapping, in which data were used directly in NetCDF files format [53]. Climate Model Data for Hydrologic Modeling (CMhyd) software was downloaded from the website of SWAT [54]. Using observed data of 30 years and historical data for 30 years, future data were downscaled. In CMhyd software, the overlapping of observed and historical data periods should be as long as 20 to 30 years [53]. If the overlapping period is less than 10 years, then results will not be considered appropriate. More details of the CMhyd software can be found in the user manual of CMhyd [53]. Change factor downscaling is showing some ambiguities and the correlation factor is less than 0.75 for model "ACCESS". On the other hand, there is a strong correlation for model "ESMLR" i.e., more than 0.8 in all three downscaling techniques.

SWAT Model Description and Setup
The SWAT model [55], developed by the United States Department of Agriculture, is broadly being used for the assessment of the footprint of climate change in water resources studies. The main objectives of SWAT model development are to simulate the quality and quantity of surface and groundwater and to predict the environmental impact of land use, land management practices, and climate change [56]. Previously, numerous researchers utilized the SWAT model for hydrological modeling studies in the various basin and sub-basin scales [29,31,45,[57][58][59][60][61]. The digital elevation model (DEM) is used for watershed delineation to investigate several properties of the basin, such as soil, slope, elevation, length of flow, streams created the longest path, etc. SWAT requires specific information on water abstraction that includes climate, land use, and land cover practices, topography, and river basin management to simulate hydrological processes. Similar properties within sub-basins can be combined into HRUs. The HRUs are spatially separated, and their replication results are added directly to the sub-basins for the entire basin path. In addition to the frequently used HRU discretization, the SWAT can also be discretized based on the groups or grids [56]. To perfectly reproduce the transfer of deposits, pesticides, and nutrients, the hydrological sequence, as prognosticated with the model, should define the historical trends of the basin. The detail about the model setup can be found in the user manual of ArcSWAT [60]. Figure 1 and Table 4 provide the details of slope classes in the Mangla Watershed. Three slope classes were selected. The first-class ranges from 0 to 20%, the second class ranges from 21 to 40%, the third class ranges from 41 to 60%, fourth ranges from 61-80%, and fifth covers the slopes greater than 80%. First-class shows the minimum slope, while the fifth class slope represents the hilly areas having higher slopes or mountainous areas. The overall slope within the watershed is higher, and due to greater slopes, the accurate and precise specification of elevation is mandatory to effectively study the flow regimes within a partially or fully snow-covered watershed. In this study, to the precise specification of the elevation, the concept of elevation bands was employed.

Elevation Bands
The precise specification of alterations in elevation plays a crucial role in the hydrological modeling studies. Elevation bands play a significant role in the SWAT model for streamflow generation within a fully or partially snow-covered region, as the snow usually melts first at the lower elevation and later at a higher elevation. Similarly, for precipitation, it will fall as fresh snow at a higher altitude, and as rainfall on lower altitudes [62]. As a result of the interdependence of snow-fall and rainfall on elevation, the SWAT model simulations by considering the elevation bands will generate accurate results. The upper limit of the elevation bands is 10 in SWAT 2012, so in this study, ten elevation bands were incorporated to get the possible accuracy in results. The inter-comparison between the results obtained by consideration of elevation bands and without elevation bands was studied by [63] and interpreted that the results in fully or partially snow-covered watershed will not be accurate without considering elevation bands.

Model Accuracy Criteria
The performance of a model can be estimated by comparing the observed streamflow and model-simulated streamflow. There are several parameters to check the performance of the model such as coefficient of determination, the efficiency of Nash-Sutcliffe, and percentage bias. By looking at the observed and calibrated hydrographs, a researcher can recognize the adaptability of the model if the model is too fast or is underestimated. However, to quantitatively evaluate the performance of the model, numerical quantities related to the efficiency of the model are required.
Coefficient of determination (R 2 ): The coefficient of determination R 2 estimates the percentage of distinction in the observed data that is introduced in the results of the simulated model.
Efficiency coefficient Nash-Sutcliffe: The efficiency coefficient used by Nash-Sutcliffe to evaluate the future predictive strength of almost all hydrological models. The main limitation of the Nash-Sutcliffe efficiency includes the fact that the variations between the observed outflow and the simulated outflow values are determined by the squared values. This leads to an overestimation of the performance of the model through outflow peaks and an underestimation of the basic output flows. To check the model performance, the following parameters are considered: If the NSE value for calibration and validation is above 0.65, then the results are very good. However, if the NSE values fall below 0.5, then the results must not be considered satisfactory for the strong footing of the model [64][65][66][67]. If the PBIAS result of the model is less than 0.1, then the results are very good and increase in the values to 0.15. The results remain in an acceptable range, but when the values of PBIAS cross 0.25, then the result falls into an unsatisfactory range; calibration must be repeated by using different parameters along with a different parametric range [68,69]. If the Pearson correlation coefficient value is between 0.81 and 1, then the model performance will be very strong. If the performance confidence index is between 0.75 and 0.85, then the model performance will be very strong [67].

Model Calibration and Validation
SWAT model calibration was done using SWAT-CUP SUFI2 Algorithm for a period of 15 years, starting from January 1981 to December 1995. The objective function R 2 was used during the calibration period. Three years were set-up as a warm-up period. The same parameters along with their values and the same number of iterations are further used to validate the model. The validation was done for a period of 15 years, from January 1995 to December 2010. Several parameters such as p-factor, r-factor, R 2 , NSE, PBIAS were calculated to check the performance of the SWAT model. P-factor represents the percentage of gauged data wrapped by the simulated 95% uncertainty band. R-factor represents the thickness of the 95% uncertainty band. The values of the p-factor should be greater than 0.70, whereas the values of the r-factor should be closer to zero. The range of p-factor is 0 to 1, whereas the range of r-factor lies between 0 and infinity [70].

Model Calibration without Considering Elevation Bands
The results of calibration and validation without considering elevation bands in Mangla watershed show a poor model performance, as represented in Figure 4. SWAT model calibration and validation performance evaluation indices, such as p-factor, r-factor, R 2 , NSE, and PBIAS, did not result in an acceptable range [67]. Figure 4 depicts the poor correlation between observed and simulated flows, as these results were simulated without considering elevation bands. The results of calibration without considering elevation bands for p-factor, r-factor, R 2 , NSE, and PBIAS were 0.43, 1.46, 0.68, 0.46, and −32. The results of the model performance indices, such as p-factor, r-factor, R 2 , NSE, and PBIAS, show the poor performance of the model in the Mangla Watershed without the use of elevation bands. So, in this research we recommend the hydrological modeling researchers to consider precise elevation bands to obtain accurate simulated results.  The results of the model performance indices, such as p-factor, r-factor, R 2 , NSE, and PBIAS, show the poor performance of the model in the Mangla Watershed without the use of elevation bands. So, in this research we recommend the hydrological modeling researchers to consider precise elevation bands to obtain accurate simulated results.

Calibration and Validation Considering Elevation Bands
Calibration and validation results revealed a strong model performance that can be potentially utilized to investigate the impacts of climate change on stream outflows ( Figure 5). A strong relationship between model-simulated and observed flows was exhibited, which shows a strong footing of the SWAT model for future flow projections. The values for p-factor, r-factor, R 2 , NSE, PBIAS are 0.77, 0.99, 0.80, 0.78, and 1.1, respectively. These results show the strong footing of the model for future prediction. The R 2 value is 0.80, which shows a strong relationship between the simulated and observed flow. Pbias should be between −20 to +20, and for perfect model footing, it should be close to zero. The pbias value for calibration is closer to the ideal, which is 1.1, which shows a great model performance for future projections. The NSE value in SWAT-CUP results is 0.78, which also represents the accuracy of model calibration.

Calibration and Validation Considering Elevation Bands
Calibration and validation results revealed a strong model performance that can be potentially utilized to investigate the impacts of climate change on stream outflows ( Figure 5). A strong relationship between model-simulated and observed flows was exhibited, which shows a strong footing of the SWAT model for future flow projections. The values for p-factor, r-factor, R 2 , NSE, PBIAS are 0.77, 0.99, 0.80, 0.78, and 1.1, respectively. These results show the strong footing of the model for future prediction. The R 2 value is 0.80, which shows a strong relationship between the simulated and observed flow. Pbias should be between −20 to +20, and for perfect model footing, it should be close to zero. The pbias value for calibration is closer to the ideal, which is 1.1, which shows a great model performance for future projections. The NSE value in SWAT-CUP results is 0.78, which also represents the accuracy of model calibration.
The results of R 2 , NSE, pbias, p-factor, and r-factor, for calibration and validation, are presented in Table 5. These results depict that the SWAT model performance in the Mangla watershed is very strong [67]. The initial values of some parameters, such as curve number 2 for soil conservation services (CN2), the alpha factor for base flow in bank storage (days) (ALPHA_BF), delay in groundwater in days (GW_DELAY), the minimum depth of water in the shallow aquifer essential for backflow (mm) GWQMN, and groundwater revap coefficient (GW_REVAP), were taken from the literature [30]. The parameters used during the model calibration are shown in Table 6. The northern part of the Mangla Watershed is a partially snow-covered region [32] so some snow cover parameters, such as the temperature of snowfall (SFTM), the Base temperature of snowmelt (SMTMP), the maximum rate of snowmelt over a year (SMFMX), the minimum rate of snowmelt over a year (SMFMN), the minimum amount of snow water resembles 100% of snow cover (SNOCOVMX), and the volume of snow that corresponds to 50% of snow cover (SNO50COV), were also utilized. The minimum initial value given to SFTM is −5 and the maximum value given is 5. The parameter SNOCOVMX was initially set between 0 and 400, whereas the minimum and maximum values designated to SNO50COV were 0.1 and 0.6. More details about these parameters can be found in the user manual of SWAT-CUP_2012 [70]. The results of R 2 , NSE, pbias, p-factor, and r-factor, for calibration and validation, are presented in Table 5. These results depict that the SWAT model performance in the Mangla watershed is very strong [67]. The initial values of some parameters, such as curve number 2 for soil conservation services (CN2), the alpha factor for base flow in bank storage (days) (ALPHA_BF), delay in groundwater in days (GW_DELAY), the minimum depth of water in the shallow aquifer essential for backflow (mm) GWQMN, and groundwater revap coefficient (GW_REVAP), were taken from the literature [30]. The parameters used during the model calibration are shown in Table 6. The northern part of the Mangla Watershed is a partially snow-covered region [32] so some snow cover parameters, such as the temperature of snowfall (SFTM), the Base temperature of snowmelt (SMTMP), the maximum rate of snowmelt over a year (SMFMX), the minimum rate of snowmelt over a year (SMFMN), the minimum amount of snow water resembles 100% of snow cover (SNOCOVMX), and the volume of snow that corresponds to 50% of snow cover (SNO50COV), were also utilized. The minimum initial value given to SFTM is −5 and the maximum value given is 5. The parameter SNOCOVMX was initially set between 0 and 400, whereas the minimum and maximum values designated to SNO50COV were 0.1

Performance Evaluation of The SWAT Model in Terms of Low and High Flows
The simulated low streamflow and high streamflow by the SWAT model were compared with observed low stream flows and observed high streamflow, respectively. The results indicated a strong correlation coefficient of 0.94 and 0.98 for low and high streamflow, respectively. The values of the coefficient of determination (R 2 ) also indicated that the SWAT model performed well in simulating low and high streamflow. The R 2 values were 0.88 and 0.96 for low streamflow and high streamflow respectively. Figure 6 shows the comparison between observed and simulated low streamflow and observed and simulated high streamflow. The values of the coefficient of determination (R 2 ) also represent that the observed and simulated low flows and observed and simulated high flows have a strong relationship. These indices depict that the SWAT model is performing well in simulating low flows and high flows.  Figure 6 shows the comparison between observed and simulated low

Seasonal Change in Temperature and Precipitation
The variation in mean seasonal rainfall is limited in spring and winter while shifting in rainfall is more during autumn and summer seasons for five GCMs. The rainfall during the moon-soon season is predicted to be more intense, while there will be a large increase in temperature; up to 8.4 °C.

Seasonal Change in Temperature and Precipitation
The variation in mean seasonal rainfall is limited in spring and winter while shifting in rainfall is more during autumn and summer seasons for five GCMs. The rainfall during the moon-soon season is predicted to be more intense, while there will be a large increase in temperature; up to 8.4 • C. From Figure 7, it is concluded that in the 20th century the temperature is relatively less than in the 21st century; meanwhile, the temperature is increasing at a relatively steeper rate as compared with the previous period. In the year 1990, the average temperature of the whole watershed was 24.5 • C, whereas the expected average annual temperature for the 2030 scenario is 26.2 • C. It can be observed that there is almost a 1.7 • C increase in temperature in four decades. However, it can also be observed that in 2070, the average temperature will reach 29.54 • C, as there will be a rise of almost 3.34 • C in the next four decades. For 2100, the temperature will touch 32.08 • C, from which it can be observed that the temperature will rise to 2.54 • C during the period from 2070 to 2100. Considering RCP (4.5), it can be observed that temperature will tend to rise at a higher rate during the middle of the century, as compared to the end of the century.  Figure 7 revealed the average temperature from 1990 to 2100 within the transboundary Mangla Watershed. The best adopted downscaling technique for the watershed, i.e., distribution mapping downscaling, was used to predict these results. Figure 7 represents that the maximum rise in temperature till 2100 will be 8.9 °C. It is a large increase in temperature which will definitely increase evapotranspiration and also cause glaciers to melt and will severely affect freshwater resources in the future.  Watershed. The best adopted downscaling technique for the watershed, i.e., distribution mapping downscaling, was used to predict these results. Figure 7 represents that the maximum rise in temperature till 2100 will be 8.9 • C. It is a large increase in temperature which will definitely increase evapotranspiration and also cause glaciers to melt and will severely affect freshwater resources in the future. Figure 8 shows the projected change in precipitation (mm) over the Mangla watershed from the year 1990 to 2100. In 1990, there was less precipitation in most of the watershed, and it increases continuously in the future from 2021 to 2100. The northern part of the Mangla Watershed located at a higher altitude registered more rainfall as compared to the southern part of the watershed. In the late century, there will be more precipitation due to an increase in temperature. The increase in temperature ultimately causes evapotranspiration and precipitation to surge.
Atmosphere 2020, 11, x; doi: FOR PEER REVIEW 16 of 33 a higher altitude registered more rainfall as compared to the southern part of the watershed. In the late century, there will be more precipitation due to an increase in temperature. The increase in temperature ultimately causes evapotranspiration and precipitation to surge.    Figure 9a shows the average annual precipitation of ensembled GCM models [71] from 1981 to 2010 under three downscaling techniques for the RCP 4.5 radiative forcing emission scenario. The average annual precipitation under change factor downscaling ranges from 3.2 mm to 5.5 mm. The average annual precipitation from ensembled model output ranges from 2.78 mm to 4.89 mm under distribution mapping downscaling, and 2.54 mm to 9.4 mm under linear scaling. Figure 9b shows the GCMs ensembled average annual precipitation under RCP 8.5. The change factor downscaling ranges between 3.14 mm and 5.47 mm, distribution mapping ranges from 2.44 mm to 4.65 mm, and linear scaling ranges between 2.87 mm and 9.33 mm.
Atmosphere 2020, 11, x; doi: FOR PEER REVIEW 17 of 33 distribution mapping downscaling, and 2.54 mm to 9.4 mm under linear scaling. Figure 9b shows the GCMs ensembled average annual precipitation under RCP 8.5. The change factor downscaling ranges between 3.14 mm and 5.47 mm, distribution mapping ranges from 2.44 mm to 4.65 mm, and linear scaling ranges between 2.87 mm and 9.33 mm.

Future Flow Projection
The streamflow generation by applying CMIP5 GCMs data in SWAT is quite uncertain. The increase in streamflow on a larger scale has been observed in all five GCMs. The percentage increase in flow following RCP 4.5 can be observed in Table 7. In the summer season, the flow will tend to decrease, while in winter, it tends to rise. Both flood and drought scenarios can be observed in the prospect in the watershed. Using the global climatic model, ACCESS, from 2021 to 2039, the annual flow will tend to rise by 60.19%. The stream flows are predicted to rise by 87.62% during the years 2040-2069 and 97.17% from 2070-2099. During the summer season (June, July, and August), these flows tend to decrease considerably for 2021-2039 scenarios, as the flows will be 121.74% more than base period flows, whereas, from 2040-2069, scenarios showed a 62.35% decrease in streamflow as compared with the base period. In the last decade, the flows in summer will reach to 52.02%. Winter season (including December, January, and February) showed a remarkable increase in streamflow, 35.57% in 2021-2030, 154.25% in the mid-century, and 200.05% in the late century. The annual and seasonal flow exceedance percentage was described in Table 7, regarding five GCMs under RCP 4.5 during three periods; 2021 to 2039, 2040 to 2069, and 2070 to 2099. The annual flow will tend to increase up to 155.2%; the flow during winter will increase up to 186.29%. In summer, it will increase up to 121% and in autumn it will increase up to 175.97%. There is a maximum increase during the spring season which may lead to floods, and during the winter season, the streamflow tends to shrink, which may cause droughts. Flows will upsurge at an abrupt rate in the prospect. For instance, it is revealed that, through the year 2021 to 2040, because of ACCESS, the annual average flows will rise to 86.36%, while for the year 2030-2070, it will rise at a higher grade of 88.52%. The flow will

Future Flow Projection
The streamflow generation by applying CMIP5 GCMs data in SWAT is quite uncertain. The increase in streamflow on a larger scale has been observed in all five GCMs. The percentage increase in flow following RCP 4.5 can be observed in Table 7. In the summer season, the flow will tend to decrease, while in winter, it tends to rise. Both flood and drought scenarios can be observed in the prospect in the watershed. Using the global climatic model, ACCESS, from 2021 to 2039, the annual flow will tend to rise by 60.19%. The stream flows are predicted to rise by 87.62% during the years 2040-2069 and 97.17% from 2070-2099. During the summer season (June, July, and August), these flows tend to decrease considerably for 2021-2039 scenarios, as the flows will be 121.74% more than base period flows, whereas, from 2040-2069, scenarios showed a 62.35% decrease in streamflow as compared with the base period. In the last decade, the flows in summer will reach to 52.02%. Winter season (including December, January, and February) showed a remarkable increase in streamflow, 35.57% in 2021-2030, 154.25% in the mid-century, and 200.05% in the late century. The annual and seasonal flow exceedance percentage was described in Table 7, regarding five GCMs under RCP 4.5 during three periods; 2021 to 2039, 2040 to 2069, and 2070 to 2099. The annual flow will tend to increase up to 155.2%; the flow during winter will increase up to 186.29%. In summer, it will increase up to 121% and in autumn it will increase up to 175.97%. There is a maximum increase during the spring season which may lead to floods, and during the winter season, the streamflow tends to shrink, which may cause droughts. Flows will upsurge at an abrupt rate in the prospect. For instance, it is revealed that, through the year 2021 to 2040, because of ACCESS, the annual average flows will rise to 86.36%, while for the year 2030-2070, it will rise at a higher grade of 88.52%. The flow will upsurge by 110.55% in the course of the year 2070 to 2100. In the summer season, a decrease was observed, whereas, in spring and winter, flows are showing a remarkable increase. By studying the results of MIROCESM, it can be observed that, on an annual basis, streamflow will increase by 121.56% from 2021-2039, as compared with the base period 1981 to 2010. In the year 2040-2070, the streamflow will increase more and reach the peak value of 137.31%, and at the end of the century, it will touch 142.87%. In the summer season, stream flows are showing a decline, as in the first 30 years, the average streamflow is projected to be 107% more than the base streamflow; in the next period, it will decline to 99.13%, and in the year 2070-2099, it will decline by more than 97.95%. Figure 11a shows the percentage increase in stream flows as compared with the base flow for the RCP 4.5 scenario by downscaling using the change factor technique. The streamflow under ACCESS increment ranges from −67.9% to 112.2%, ESMLR ranges from −12.2% to 111.3%, and CCSM4 ranges from −20.9% to 96.3%; HadGEM2 ranges from −54.6% to 56.1% and MirocESM ranges from −49.6% to 67.5%. In Figure 11b, an annual percentage increase in stream flows, as compared with the base flow by using the change factor downscaling technique under RCP 8.5, is presented. The extreme values for streamflow increase under ACCESS ranges from 9.03% to 243.2%, ESMLR arrays from 28.09% to 165.33%, CCSM4 ranges from 31.3% to 165.3%, HadGEM2 ranges from 28.09% to 166.8% and MirocESM series from 48.7% to 239.26%. Figure 11c shows an annual percentage increase in stream flows as compared with baseflow under the linear scaling downscaling technique for RCP 4.5.
The extreme values for streamflow under ACCESS increment ranges from −54.3% to 145.3%, ESMLR ranges from −87.9% to 345.3%, CCSM4 ranges from −45.5% to 112.7%, HADGEM2 ranges from −27.9% to 53.9%, and MirocESM ranges from −51.7% to 220.3%. Figure 11d shows an annual percentage increase in stream flows as compared with the base flow by using linear scaling downscaling under RCP 8.5. The extreme values for streamflow under ACCESS increment range from −47.6% to 113.2%, ESMLR ranges from −72.0% to 172.3%, CCSM4 ranges from −33.5% to 103.3%, and HadGEM2 ranges from −54.2% to 99.4%. Figure 11e shows an annual percentage increase in stream flows as compared with the base flow by using distribution mapping downscaling under RCP 4.5. The extreme values for streamflow under ACCESS increment range from −63.7% to 139.2%, ESMLR ranges from 44.8% to 140.73%, CCSM4 ranges from −52.7% to 114.3%, HadGEM2 ranges from −46.3% to 168.9% and MirocESM ranges from −44.8% to 140.7%. Figure 11f shows an annual percentage increase in stream flows as compared with the base flow, by using distribution mapping downscaling under RCP 8.5. The extreme values for streamflow under ACCESS increment ranges from −61.5% to 112.7%, ESMLR ranges from −53.6% to 112.7%, and CCSM4 ranges from −54.3% to 96.6%, HadGEM2 ranges from 32.49% to 164.425% and MirocESM ranges from −45.7% to 56.9%. It is observed that the possibility of the appearance of low and high flow could be more leading in the prospect within the Mangla watershed following both RCP scenarios. However, average flows during the summer season were predicted to decrease by all five GCMs. This reduction in summer outflow is usually possible due to the reduction in average yearly precipitation. The difference in high outflow and outflow duration trajectories explains that the repetition of water abundance will further increase their quantities into the prospect that will generate plenty of administration obstacles in the watershed. Not only economic losses but also losses of life will be caused by flooding in the future.      Figure 14 shows the average annual flow (January to December); winter (October to March), Summer (August to September); Winter (December, January, and February); Spring (March, April, and May); Summer (June, July, and August) and Autumn (September, October, and November) during the year 2021 to 2100. Figure 14 revealed that there will be more streamflow during the summer season, whereas there will be very little streamflow during the autumn season. The streamflow during the summer season is relatively more due to the melting of glaciers but as we proceed to the next period 2040-2069 the streamflow during the summer season show a remarkable declining trend. A large volume of water in streams is observed in the spring season and it will tend to rise in the future. By considering distribution mapping, the average annual flows are 1887 m 3 /s flows during the autumn season, 1435 m 3 /s in the summer season 2263 m 3 /s, and in the spring season, stream flows will be 2280 m 3 /s. There are a few uncertainties in the generated flows because the climatic data of temperature and precipitation from several GCMs vary due to their different resolutions. All of the five GCM are predicting huge stream flows during the spring season, while there will be very little streamflow during the autumn season. The multi-model ensemble means (EM) method was used to ensemble the five GCMs time series data under three types of downscaling techniques and two different emission scenarios. Figure 13 revealed that in EM, under change factor downscaling, GCMs are projecting high flows, ranges between 735 m 3 /s to 2169 m 3 /s under RCP 4.5 and 662 m 3 /s to 2149 m 3 /s under the RCP 8.5 scenario.  Figure 14 shows the average annual flow (January to December); winter (October to March), Summer (August to September); Winter (December, January, and February); Spring (March, April, and May); Summer (June, July, and August) and Autumn (September, October, and November) during the year 2021 to 2100. Figure 14 revealed that there will be more streamflow during the summer season, whereas there will be very little streamflow during the autumn season. The streamflow during the summer season is relatively more due to the melting of glaciers but as we proceed to the next period 2040-2069 the streamflow during the summer season show a remarkable declining trend. A large volume of water in streams is observed in the spring season and it will tend to rise in the future. By considering distribution mapping, the average annual flows are 1887 m 3 /s flows during the autumn season, 1435 m 3 /s in the summer season 2263 m 3 /s, and in the spring season, stream flows will be 2280 m 3 /s. There are a few uncertainties in the generated flows because the climatic data of temperature and precipitation from several GCMs vary due to their different resolutions. All of the five GCM are predicting huge stream flows during the spring season, while there will be very little streamflow during the autumn season. The box-plot represents the value of Q0 as lower extreme, Q25% as lower quartile, Q50% as median, Q75% as upper quartile, and Q100% as the upper extreme of average annual flows. Figure  15  The box-plot represents the value of Q0 as lower extreme, Q25% as lower quartile, Q50% as median, Q75% as upper quartile, and Q100% as the upper extreme of average annual flows. Figure 15 Table A1 represents the peak, mean, and low streamflow from the year 2021 to 2100 under the RCP-4.5 emission scenario. The peak flows occur commonly in the summer season (during June, July, and August) due to the melting of snow from the northern part of the watershed, and low flows usually occur during the winter season due to less snow melting. HadGEM2 depicts the maximum

Discussion
In this study, future climate change was assessed under RCP (8.5) and RCP (4.5) emission scenarios. Change of temperature and precipitation and their probable effects on water resources at Mangla Watershed was studied on a broad spectrum. Five GCMs, along with their ensembles, were selected, and three types of downscaling techniques were compared to estimate the change in precipitation and temperature. After a comparison of these downscaling techniques with base data, it was concluded that the distribution mapping is the best downscaling technique to study climatic parameters within the watershed, however, linear scaling downscaling was found to be a less effective downscaling technique, as the results are in close agreement with previous studies [51]. Downscaled data were used to generate stream flows using a semi-distributed hydrological model SWAT. In this research, evapotranspiration was neglected due to data scarcity and unreliability. Observed data were gathered from seventeen metrological stations, five of them in the Indian region. The scarcity of metrological stations within Mangla Watershed can cause the low performance of the calibration. The land use classification and soil type were kept constant throughout the simulations. This type of assumption can affect the prediction of streamflow. Uncertainty exists as a well-known problem toward most of the hydrological modeling considerations, particularly on large scales. The bulk part of these uncertainties appears in the forecasted precipitation and temperature datasets, due to the distinctive resolutions of several GCMs. The streamflow generated by using these datasets also shows uncertainties. In the prospect, as interpreted by [31], the projection of climatic parameters using CMIP5 GCM is considerably uncertain. In the present research, five GCMs of distinct variation in resolutions are used to predict future climatic scenarios. Many uncertainties were found in the maximum and minimum temperature, as well as in rainfall after using three downscaling techniques applied on CMIP5 GCMs. The uncertainties occurred in maximum temperature using the RCP 4.5 scenario ranges between 2.58 • C and 5.30 • C, while in minimum temperature, it ranges between 0.8 • C and 2.0 • C. These uncertainties under the RCP 8.5 scenario are even more probable. The maximum temperature under the RCP 8.5 scenario ranges between 4.33 • C and 8.90 • C, whereas the minimum temperature ranges between 3.10 • C to 8.60 • C [21]. These temperature variations ultimately affect the amount of precipitation. The percentage increase in precipitation compared to the baseline period ranges between 45 to 128.4% under the RCP 4.5 scenarios and 92 to 180.3% in the RCP 8.5 scenario. As predicted by [72], it is hard to exactly predict the climatic variables, particularly rainfall, and consequently, uncertainties in predicted temperature and rainfall will lead to uncertainties in streamflow [31]. The results of five GCMs under both RCPs 4.5 and RCP 8.5 scenarios, illustrates that the change in temperature and precipitation frequencies will affect the water resources and streamflow within the transboundary, directly or indirectly. Besides climate change, other anthropogenic changes can also have a significant impact on streamflow, such as future management practices and land-use change significantly influencing water availability [61,73].
The projected GCMs show an increase in the precipitation amount in the study area compared to base years and are in close agreement with previous studies [21,24,74]. The rainfall is greatly affected by various transmission systems, such as monsoons [75]. The moisture transportation, inter-annual, and inter-decadal change of the monsoon are the principal agents that can induce an accession in the rainfall frequencies in the Indus basin region [45]. An increase in precipitation amount is also vulnerable to temperature change at a different inter-annual and inter-decadal time [29]. These factors are performing a major role in the increase of precipitation; further investigations by performing different climatic modeling techniques ought to be carried in the prospect to study these aspects in more detail. The projected increase in precipitation and temperature tends to influence the streamflow within the transboundary Mangla dam watershed and is in close agreement with the previous work carried out in the study area [21,24,45]. This increase in streamflow in the rivers will tend to cause floods [76], so several preventions, such as the construction of new dams, are strongly recommended. The results of this study are consistent with the results of precipitation and temperature in the past [45]. These projected climatic variables and stream flows could be valuable for policymakers and water resource administrators for the proper administration of water reserves. Some of the predicted streamflows showed a large increase that may cause floods in the prospect, to avoid these floods this study can be valuable for water resource planners for the management of water resources construction of new water storage structures, and modernizing the designs and storage capabilities of existing dams.

Conclusions
Pakistan stands among one of the countries with significant and more prominent water stress, and its existing water resources are extremely vulnerable to climate change. In this investigation, the possible consequences of climate change on water reserves in the Mangla Watershed have been evaluated under RCP 4.5 and RCP 8.5 scenarios using the outputs of five GCMs. The forecasts for climatic variables and their probable influences on water resources within the boundary of Mangla Watershed were examined up to the year 2099. Under CMIP5, five GCMs and their ensembles were used, considering RCP 4.5 and RCP 8.5 emission scenarios to forecast temperature, precipitation, and streamflow within the watershed, to attain substantial precipitation, maximum and minimum temperature data series following various climate change situations. Change factor, linear scaling, and distribution mapping downscaling techniques were exercised to downscale the future projections of climatic parameters. Furthermore, the mentioned three downscaling techniques were inter-compared and decided the best downscaling technique for the watershed. The hydrological model SWAT, after calibration, was applied to produce streamflow based on the values of the climatic parameters that were downscaled by the specified downscaling techniques. Subsequently, the streamflow generated by using these downscaled data in SWAT was analyzed, comprehensively, for all five models and their ensembles. The several noticeable results extracted from this research were compiled as subsequently:

•
The results regarding calibration and validation criteria, i.e., NSE, PBIAS, and R 2 , are satisfactory within the monthly period by incorporating elevation bands. The calibrated hydrological model SWAT precisely re-generates streamflow within the Mangla Watershed.

•
Maximum temperature and minimum temperature are projected to increase in the future from 2021 to 2099 for all five GCMs, under both RCP (4.5) and RCP (8.5) emission scenarios.

•
The projected precipitation is more uncertain and obscure. All five GCMs and their ensemble are predicting an increase in precipitation frequencies from June to September within the Mangla Watershed with a significant increase of (219%) in June.

•
The projected average annual flow will increase constantly for all five GCMs and their ensemble, especially in the Spring season with a 193% increase, but there will be a decrease in the streamflow during the autumn season. Moreover, there will be an excess of high flows and low flow events for the projected flows.
Hence, based on these outcomes, to predict climatic parameters in the future, a more certain GCMs ensemble scenario could be used. An increase in flow will cause high flood risk in the watershed. By applying good management practices and by constructing dams and reservoirs, Pakistan may be able to generate more hydropower and also produce more food for its residents. The principal outcome of the research is that the Mangla watershed will be expected to face more flooding in the future, and average flows are supposed to decrease. The maximum temperature, minimum temperature, and precipitation peak will also tend to increase in the future. The outcomes of this research are helpful for climate change researchers, development planners, and policymakers toward designing and accomplishing suitable water administration systems to mitigate the influences of climate change.