Impact Assessment of Future Climate Change on Streamflows Upstream of Khanpur Dam, Pakistan using Soil and Water Assessment Tool

Abstract

The study aims to evaluate the long-term changes in meteorological parameters and to quantify their impacts on water resources of the Haro River watershed located on the upstream side of Khanpur Dam in Pakistan. The climate data was obtained from the NASA Earth Exchange Global Daily Downscaled Projection (NEX-GDDP) for MIROC-ESM model under two Representative Concentration Pathway (RCP) scenarios. The model data was bias corrected and the performance of the bias correction was assessed statistically. Soil and Water Assessment Tool was used for the hydrological simulation of watershed followed by model calibration using Sequential Uncertainty Fitting version-2. The study is useful for devising strategies for future management of Khanpur Dam. The study indicated that in the future, at Murree station (P-1), the maximum temperature, minimum temperature and precipitation were anticipated to increase from 3.1 °C (RCP 4.5) to 4.0 °C (RCP 8.5), 3.2 °C (RCP 4.5) to 4.3 °C (RCP 8.5) and 8.6% to 13.5% respectively, in comparison to the baseline period. Similarly, at Islamabad station (P-2), the maximum temperature, minimum temperature and precipitation were projected to increase from 3.3 °C (RCP 4.5) to 4.1 °C (RCP 8.5), 3.3 °C (RCP 4.5) to 4.2 °C (RCP 8.5) and 14.0% to 21.2% respectively compared to baseline period. The streamflows at Haro River basin were expected to rise from 8.7 m³/s to 9.3 m³/s.

1. Introduction

The temperature of the Earth has increased significantly in recent decades with the decade of 2000s being the warmest [1]. The high emission of greenhouse gases (GHGs), particularly carbon dioxide that traps heat, poses a crucial threat to the environment [2]. It hinders the radiating efficiency of the surface of the Earth thereby affecting the overall temperature of the Earth [3]. The continued emission of these gases will not only increase the temperature but will also instigate long-term changes in the climate [4]. Global warming will disturb the global water cycle while altering the hydrological regimes that will affect management of fisheries, ecology, industrial and municipal water supply, human health and water energy balance [5]. The hydrological systems are under great threat due to changing climatic conditions in terms of variation in the availability of water as well as increased rate of precipitation [6]. Long-term changes in climate may alter the hydrological cycle to an extent that it may affect the intensity and frequency of floods as well as droughts. According to the United Nations report [7], climate change will increase the probability of natural disasters particularly related to water resources. These disasters pose negative social impacts and largely risk the economy of the region. Climate change impacts vary from region to region [8].

Global climate change processes have drastic effects on the climate of Pakistan in social, economic and environmental aspects [9]. Germanwatch, in its “Global Climate Risk Index”, has ranked Pakistan in third place after Haiti and Philippines, as a country most impacted by climate change effects [10]. Environmental dependency for fulfilling basic living requirements, increased rate of population growth, poverty and lesser resources to adjust to the negative effects of climate changes are also contributing to changes in the climate of the country [9]. Climate change is causing a significant rise in the average mean sea level with increase in global average temperature and variation in precipitation patterns regionally [11]. In Pakistan, during the past century, the average annual temperature has been observed to increase by 0.6 °C [12]. However, the rate of change varies throughout the country. In Northern Pakistan, the rate of increase in temperature is higher as compared to the Southern region [9]. The agrarian land of Pakistan receives most of the rainfall during the months of June till September. A twenty-five percent increase in precipitation observed in the previous century has been attributed to the impacts of climate change [12]. Despite numerous research activities, it is vital to quantify the effect of change in temperature, precipitation and atmospheric concentration of carbon dioxide on the water resources. At the watershed level, this information is useful for proper planning and decision-making [3]. In the past, Pakistan has already suffered from the extreme impact of climate change in the form of drought of 1998–2001 and flood of 2010 [9].

The river flows are affected by changes in climatic conditions [13]. Akhtar et al. [14] studied climate change impact on the water resources of Hindukush-Karakorum-Himalaya in Pakistan under three different glacial coverages. They predicted an increase in temperature, precipitation and discharge with increased risk of flooding. Azim et al. [15] also predicted an increase in discharge and sediment yield for Naran watershed in Pakistan. Kunhar River basin is likely to experience increased flooding as well as droughts; low flows are projected to decrease and high flows are expected to increase in the future [6]. Meanwhile, the flows in Chitral watershed in Pakistan are expected to decreases in future by up to 4.2% [13].

To assess impact of climate change on the water resources, climate model data is integrated into hydrologic models. The hydrological processes are simulated by linking the water yield and the climate changes. Based on the selected hydrological model and the climate scenario, the results of the study may greatly vary [16]. The Soil and Water Assessment Tool (SWAT) is a semi-distributed, continuous-time, physical model facilitating long-term hydrologic simulations. Although SWAT is not a three-dimensionally distributed model, it is comprised of spatially separated parts forming sub-basins and hydrologic response units (HRUs) [17]. The effects of changes in carbon dioxide gas, temperature, precipitation and humidity can be related to the plant growth and its effect on the evapotranspiration, snow and runoff can be evaluated [16]. Recently, SWAT is used in many studies for the simulation of watershed’s hydrological processes with focus on climate change impact assessment. SWAT was used by Abbaspour et al. [18] to study climate change impact on Iran’s water resources. Garee et al. [19] studied climate change impact on the streamflows of Hunza watershed in Pakistan. According to their study, the overall temperature is anticipated to rise form 1.39 °C to 6.58 °C with 31% increase in precipitation causing 5–10% increase in streamflows by the end of the century.

The current study aimed to quantify the probable effects of climate change on the water resources of Haro River watershed in Pakistan by using SWAT. Haro River, flowing partly through Punjab and Khyber Pakhtunkhwa provinces in Pakistan, is located on the upstream side of Khanpur Dam. Khanpur dam is a major source of supplying drinking water to the twin cities of Islamabad and Rawalpindi. It also fulfils the irrigation demand for Punjab and Khyber Pakhtunkhwa. It is important to study the effects of future climate projection on the streamflows of Haro River watershed, as it will affect the flow at Khanpur Dam. The variation in future streamflows will aid in devising future policies and management measures for Khanpur Dam. To achieve the study objectives, bias corrected climate projections for the baseline period (1976–2005) and future period (2006–2095) from NEX-GDDP model MIROC-ESM, under two emission scenarios i.e., RCP 4.5 and RCP 8.5, were integrated in a calibrated and validated SWAT model. The study area and data description are discussed in Section 2. The methodology adopted to achieve the study goals is presented in Section 3. Section 4 includes the results and discussion on the results. The conclusion, limitations and the recommendations are presented under Section 5.

2. Study Area and Data Description

2.1. Study Area

Haro River flows partly through the provinces of Punjab and Khyber Pakhtunkhwa. It originates from Changla Gali and Nathia Gali and is drained by Lora Haro and Stora Haro. The Lora Haro drains through Murree hills whereas Stora Haro drains from the hills of Nathia Gali. At Jabri, the Lora Haro merges into Stora Haro. Haro River flows through narrow and deep gorges and is dammed at Khanpur [20] (Figure 1). The Dam was built in 1983 with the primary purpose of supplying drinking water to Islamabad and Rawalpindi as well as to the industrial units at Taxila. Throughout the year, the springs and streams fulfil the demand of the domestic water supply. At the time of droughts, the stream discharge is significantly reduced and the seepage through the spring retains the perennial flow [21]. The major source of streamflows in the region is precipitation.

2.2. Data Description

2.2.1. Observed Climate Data

The daily data for maximum and minimum temperature and precipitation for the two nearby stations at Murree (P-1) and Islamabad (P-2) from 1987–2005, was obtained from Pakistan Meteorological Department (PMD). The wind speed, relative humidity and solar radiation data was obtained from National Centre for Environmental Prediction’s (NCEP) Climate Forecast System Reanalysis (CFSR) database available online under Global weather data for SWAT [22].

2.2.2. Climate Model Data

Climate projections, including precipitation and minimum and maximum temperature were obtained for the MIROC-ESM model under NEX-GDDP. The projections were divided into baseline period (1976–2005) and three future periods: early century or 2020s (2006–2035), mid-century or 2050s (2036–2065) and late century or 2080s (2066–2095).

2.2.3. Discharge Data

For Dhartian gauging station, the discharge data was acquired from the Surface Water Hydrology Project (SWHP) of the Water and Power Development Authority (WAPDA) from the year 1989 until 1991 and from 1996 until September 1998. The former period was used for the calibration of SWAT whereas the later period was used for the validation of SWAT.

2.2.4. Spatial Data

The spatial data, including digital elevation model (DEM), landuse data, and soil data, was incorporated in the SWAT. All of the spatial dataset was projected in similar coordinate system i.e., WGS_1984_UTM_Zone_43N.

The DEM of the study region was used by SWAT to develop the stream network and to delineate the watershed. It was acquired online from the database of CGIAR-Consortium for Spatial Information v4.1 [23]. It included DEM from Shuttle Radar-Topography Mission (SRTM) at the resolution of 90 m. Haro River lies mostly in the mountainous terrain with elevation ranging from 2786 m to 627 m above mean sea level.

Landuse data was obtained from the Water Resource Engineering Division of National Engineering Services Pakistan (NESPAK). The landuse was classified into six categories including forest, hills, grasses, open/barren land, agricultural land and water. The soil classification of the region was obtained from Harmonized World Soil Database (HWSD) version 1.2 available online [24]. The soil was classified as Eutric Cambisol with clay texture having soil name Be72-3c-3672. Based on the DEM, three slope classes were defined in SWAT model i.e., gentle slope, medium slope and steep slope.

3. Methodology

The methodology adopted to obtain the research objectives included the following steps:

• Selection of a climate model.
• Bias correction of baseline MIROC-ESM model data using linear scaling.
• Performance evaluation of bias correction method using mean, standard deviation, 90th percentile and 10th percentile after dividing the bias corrected baseline model data into calibration and validation periods.
• Bias correction of future model data using correction factors obtained in step 2.
• Comparison of baseline and future climate projections on monthly time step.
• Hydrological simulation of Haro River watershed in SWAT model using observed climate data obtained from PMD.
• Calibration, validation, sensitivity analysis and uncertainty analysis for the SWAT model using Sequential Uncertainty Fitting version-2 (SUFI-2) in SWAT Calibration and Uncertainty Program (SWAT-CUP).
• Hydrological simulation using calibrated SWAT model for baseline and future climate projections by incorporating the effect of change in the level of carbon dioxide under both emissions scenarios; RCP 4.5 and RCP 8.5.

3.1. Selection of a Climate Model

Five General Circulation Models (GCMs), (ACCESS 1.0, MPI-ESM-LR, CCSM4, CNRM-CM5 and MIROC-ESM) under NEX-GDDP were analyzed based upon the value of coefficient of determination (R²) of the climate model data against meteorological station data. A single model was selected for climate change impact assessment at Haro River watershed based upon highest R² value obtained for both stations, at quarterly scale.

3.2. Bias Correction of Climate Model Data

Linear scaling was used to correct the historical and future model data for any systematic bias. The model data was corrected based upon the observed meteorological data for precipitation, and minimum and maximum temperature at P-1 and P-2 from 1987–2005. In bias correction, the correction factors were computed based upon the difference between mean monthly observed values and model simulated historical data. The correction factors were then applied to the model simulated data [25]. Precipitation was corrected using multiplicative factors whereas additive factor was used to correct mean monthly values of minimum and maximum temperature [26]. The bias correction of precipitation and temperature time series was based upon Equations (1) and (2).

$$P_{sim\left(bc\right)}=P_{sim}·\left(P_{obs}/P_{his}\right)$$

(1)

$$T_{sim\left(bc\right)}=T_{sim}+\left(T_{obs}-T_{his}\right)$$

(2)

where; $P_{sim\left(bc\right)}$ and $T_{sim\left(bc\right)}$ are bias corrected monthly precipitation and temperature respectively. $P_{sim}$ and $T_{sim}$ are NEX-GDDP values for precipitation and temperature respectively. $P_{obs}$ and $T_{obs}$ are observed values for precipitation and temperature respectively. $P_{his}$ and $T_{his}$ are NEX-GDDP simulated historical values for precipitation and temperature respectively.

Steps for bias correction:

• The period of observation (1987–2005) was divided into two periods: calibration period for bias correction (1987–1996) and validation period for bias correction (1997–2005).
• Historical model data and observed data was used to compute the correction factors for precipitation and temperature for the calibration period.
• Performance evaluation for the linear scaling method for the calibration period was evaluated based upon comparison of mean, standard deviation, 90th percentile and 10th percentile between the observed data and the model data before and after bias correction.
• The bias corrected factors from the calibration period were applied to the validation period and the performance was evaluated for the validation period.
• The monthly correction factors were applied to the future model data from 2006–2095.

3.3. Comparison of Future Climate Projection

The climate change impact on mean monthly precipitation, minimum and maximum temperature was studied by evaluating the variation in future climate variables under two emissions scenarios: RCP 4.5 and RCP 8.5 with the historical/baseline period. The future period was divided into 30-year periods: early century (2006–2035), mid-century (2036–2065) and late century (2066–2095). RCPs provide the initiation point for climate change impact research [27]. RCP 8.5, the high-energy intensive scenario, represents highest range of predicted emissions in the order of 15–20 gigatonnes of carbon by end of the century [28]. The radiative forcing pathway of RCP 4.5 is an intermediate emission scenario, representing near-term climate projections [27]. In RCP 4.5, the radiative forcing level is stabilized at 4.5 W/m² after 2100, with equivalent carbon dioxide concentration (CO₂-eq) approximately equal to 650 ppm. In contrast, under RCP 8.5, the CO₂-eq was approximated to be more than 1370 ppm with radiative forcing level more than 8.5 W/m² in 2100 [28].

3.4. Hydrological Simulation Using SWAT

The hydrological modelling of Haro River watershed was carried out using ArcSWAT 2012. In SWAT, the water balance of the river basin is considered as the major driving force behind all processes occurring in the watershed. The land phase in the model deals with surface runoff, evapotranspiration, lateral flow, tributary channel, groundwater, ponds and return flow whereas the routing phase or the in-stream phase deals with the movement of water, sediments, nutrients, organic chemicals within and through the channel [29].

SWAT requires the climate dataset as well as the spatial data including DEM, landuse information, soil information and slope classification. All the spatial dataset needs to be projected in the same coordinate system. The climate data at station P-1 and P-2 from 1987–2005 was used to run the SWAT model with the outlet of the watershed defined at Dhartian gauging station. The DEM was used to delineate the watershed into nine subbasins whereas the spatial information of landuse, soil and slope was used by the model to create 96 lumped areas or HRUs. The model output was analyzed at a monthly time step. Surface runoff in SWAT was computed using SCS Curve Number method whereas the potential evapotranspiration was computed using the Penman-Monteith method. Variable storage method was used for the routing of water through main channel.

3.5. Calibration and Validation of SWAT

SWAT was calibrated using SUFI-2 optimization algorithm. The model was calibrated for discharge at Dhartian station from 1989–1991. The calibration was done in two iterations each having 500 simulations. The model calibration was conducted until satisfactory values for Nash Sutcliffe efficiency (NSE) and R² were obtained. The validation of the model was carried out from 1996 until 1998 using the parameter ranges obtained in the calibration period. A single iteration for the validation period was carried out by having a similar number of simulations as the last iteration of the calibration period.

SUFI-2 makes use of inverse modelling technique to take into account uncertainty in parameter and carry out calibration [30]. It provides parameter ranges that ensure best results of simulation by satisfying two criteria denoted by “P-factor” and “R-factor”. P-factor ensures that the 95 percent prediction uncertainty (95PPU) computed at the 2.5 percent and 97.5 percent cumulative distribution levels of the output variables encloses maximum observed data. The thickness of the 95 PPU band is denoted by the R-factor [31]. The strength of the calibration and validation is indicated by the P-factor and R-factor. The best calibrated ranges are finalized when a balance is reached in both the factors. P-factor value ranges from 0–1 but a value greater than 0.7 is regarded as acceptable for the calibration of discharge whereas an R-factor having values lesser than 1.5 is considered adequate [32].

NSE defines the fit between the observed and simulated data on 1:1 line [33]. The NSE value ranges between 1 and −∞ [34]. According to Moraisi et al. [35], NSE value lesser or equal to 0.5 reflects unsatisfactory performance of model whereas value greater than 0.75 shows very good model performance.

Coefficient of determination or R² evaluates the strength of a linear relation between observed and simulated data [33]. Its value varies from 0–1 where R² = 0 shows that no correlation exists between observed and simulated data. A perfect correlation is indicated by a value of 1 [34]. R² value greater than 0.5 is considered acceptable [35].

3.6. Incorporating Climate Change Impacts in SWAT

The climate change impact was considered by changing the level of CO₂ and the projected temperature and precipitation under each emission scenarios for three time periods: early century, mid-century and late century. The concentration of CO₂ level for each period is given in

Table 1. Projected concentration of CO₂ in ppm.

Period	RCP 4.5	RCP 8.5
Early century	424	436
Mid-century	487	541
Late century	523	804

.

The values for early and late century were obtained from a study by Li et al. [36], whereas the mid-century values were obtained from Meinshausen et al. [37].

4. Results and Discussion

4.1. Selection of Climate Model

The results of the R² for precipitation for five GCMs: ACCESS 1.0, MPI-ESM-LR, CCSM4, CNRM-CM5 and MIROC-ESM were 33.3%, 32.9%, 10.5%, 43.4% and 40.5% respectively at station P-1 and 52.5%, 0.2%, 7.5%, 31.5% and 54.7% respectively at station P-2. Hence, due to better performance of MIROC-ESM compared to the other models, it was chosen for climate change impact assessment at Haro River watershed.

4.2. Bias Correction

Figure 2 shows the variation in statistical indicators including mean, standard deviation, 90th percentile and 10th percentile for observed, non-bias corrected and bias corrected model data for precipitation, maximum and minimum temperature at stations P-1 and P-2. The correction factors found during the calibration phase (1987–1996) were applied to the validation period (1997–2005) and the performance during both phases was evaluated. At both stations, the bias corrected parameters during calibration period indicated more improvement as compared to the validation period. For precipitation, maximum and minimum temperatures, the non-bias corrected indicators showed underestimated values specifically for mean, standard deviation and 90th percentile at P-1 and P-2. The bias corrected mean precipitation at both stations became equal to the observed mean precipitation during the calibration period whereas non- bias corrected values showed a percentage difference of −44.3% and −25.2% at P-1 and P-2 respectively in comparison to the observed values.

The indicators for maximum temperature showed overestimated values for non-bias corrected model data at station P-1 for both calibration and validation periods, whereas for station P-2, the observed, non-bias corrected and bias corrected statistical parameter values demonstrated little variation. The standard deviation values at P-1 for non-bias corrected and bias corrected cases during calibration as well as validation periods were close to the observed values. At P-1, the mean values for the maximum temperature were improved with percentage difference of +27.2% (non-bias corrected) to −0.4% (bias corrected) during the calibration period and from +24.4% to −2.2% during the validation period with respect to the observed data. The 90th percentile at P-1was improved from a difference of +22.9% (non-bias corrected) to −1.5% (bias corrected) during the calibration period and with a difference of +18.2% (non-bias corrected) to −5.5% (bias corrected) with respect to observed data during the validation period. The 10th percentile was improved with a difference of +35% to +3.26% during calibration period and +34.4% to +2.5% during the validation period, for station P-1.

For minimum temperature, the non-bias corrected values for three indicators including mean, 90th percentile and 10th percentile were overestimated during calibration and validation period at both stations. At P-1, during the calibration period, the 90th percentile for non- bias corrected minimum temperature was almost equal to the observed value with a difference of +0.8% but bias correction increased the percentage difference to −6.0%. Overall, it was observed that the application of linear scaling improved the statistical indicators significantly at both stations. This indicated that the linear scaling technique used for the improvement of the selected NEX-GDDP model output is suitable to be used for the study region.

4.3. Comparison of Mean Monthly Future Meteorological Parameters at P-1 and P-2

The mean monthly amounts of precipitation, maximum and minimum temperature for the baseline and the future periods: early century (2020s), mid-century (2050s) and late century (2080s), under RCP 4.5 and RCP 8.5 for station P-1 are shown in Figure 3 below.

Mean monthly precipitation was expected to be the highest during the month of July followed by August for baseline and future periods. The maximum difference in precipitation, for both scenarios during the month of July was expected in the mid-century from 20.4 mm/day under RCP 4.5 to 25.0 mm/day under RCP 8.5. Under both scenarios, in July, no precipitation difference was expected during the early century but during the late century, the precipitation was expected to increase by 1mm/day under RCP 8.5 compared to that under RCP 4.5. Overall, the maximum increase in precipitation up to 90.6% was expected under RCP 8.5 in late century during the month of October compared to the baseline precipitation of 3.2 mm/day.

The maximum temperature was expected to be highest during the month of June. The late century temperatures were expected to be highest in future but the increment under RCP 8.5 was substantially higher than under RCP 4.5. The highest maximum temperature under RCP 8.5 was expected to be 34 °C approximately, 29.2% higher than the baseline temperature of 26.3 °C.

Under RCP 4.5, the highest minimum temperature was expected in the late century June with the increment of 26.0% compared to the baseline minimum temperature of 16.9 °C. Under RCP 8.5, the minimum temperature early in the century was expected to be lower than that under RCP 4.5 but the late century temperatures were anticipated to be higher. The highest minimum temperature was expected during the late century June with temperature expected to be 42.6% higher than the baseline temperature of 16.9 °C.

The mean monthly amounts of precipitation, maximum and minimum temperature for the baseline and the future periods under RCP 4.5 and RCP 8.5 for station P-2 are shown in Figure 4. The precipitation, at station P-2, was also expected to be highest during the month of July followed by August. Under RCP 4.5, the early century was anticipated to receive maximum precipitation of 20.2 mm/day, approximately 20.9% higher than the baseline precipitation of 16.7 mm/day whereas the late century and mid-century were expected to receive 18% and 14% increased precipitation as compared to baseline precipitation respectively. Under RCP 8.5, the early century precipitation was projected to be lower than that under RCP 4.5 but about 19% higher than the baseline precipitation. Mid-century was likely to receive most amount of precipitation, about 38% higher than the baseline precipitation followed by late century.

Similar to the case at station P-1, the maximum temperature at P-2 was also expected to be highest during the month of June. Under RCP 4.5, the maximum temperature late in the century was expected to reach as high as 43.6 °C, almost 4.7 °C higher than the baseline temperature of 38.9 °C. The early century temperatures under RCP 8.5 were anticipated to be lesser than that under RCP 4.5 whereas the mid-century temperatures under both scenarios were anticipated to be similar. Under RCP 8.5, the late century temperature was anticipated to reach up to 47 °C, approximately 8.1 °C higher than the baseline temperature.

The minimum temperature was expected to be highest during July at station P-2. In July, during the early century, under RCP 4.5, the minimum temperature was 1.1 °C higher than the baseline temperature of 23.6 °C but 0.6 °C higher than that under RCP 8.5. During mid-century, in July, the minimum temperature was anticipated to rise under both emissions scenarios but the temperature under RCP 8.5 was still expected to be 0.2 °C lower than that under RCP 4.5. In the late century, the minimum temperature under RCP 4.5 was expected to reach 27.6 °C whereas it was anticipated to be 30.1 °C under RCP 8.5; approximately 6.5 °C higher than the baseline minimum temperature.

Table 2. Comparison of climate variables at station P-1 and P-2 under RCP 4.5 and RCP 8.5.

Station ID	Climate Variable	Time Period	Values	Increase with Respect to Baseline
P-1	Precipitation	Baseline	5.6
		RCP 4.5	6.1	8.9
		RCP 8.5	6.4	14.2
P-2		Baseline	3.9
		RCP 4.5	4.5	15.4
		RCP 8.5	4.8	23.1
P-1	Maximum Temperature	Baseline	17.7
		RCP 4.5	20.8	3.1
		RCP 8.5	21.7	4.0
P-2		Baseline	28.6
		RCP 4.5	31.9	3.3
		RCP 8.5	32.7	4.1
P-1	Minimum Temperature	Baseline	9.4
		RCP 4.5	12.6	3.2
		RCP 8.5	13.7	4.3
P-2		Baseline	13.6
		RCP 4.5	16.9	3.3
		RCP 8.5	17.9	4.2

Note: Precipitation is given in mm/day, temperature in °C. Precipitation change is given as percentage (%).

summarizes average annual climate parameters for stations P-1 and P-2 for the baseline as well as under both emissions scenarios. In the future, at P-1, the precipitation was projected to rise with respect to baseline precipitation of 5.6 mm/day. At P-2, the precipitation was anticipated to increase more under RCP 8.5 than under RCP 4.5 with respect to baseline precipitation of 3.9 mm/day.

Overall, under RCP 4.5, the precipitation varied from 0.5 mm/day to 0.6 mm/day, at P-1 and P-2 respectively whereas under RCP 8.5, the precipitation varied from 0.8 mm/day to 0.9 mm/day, at P-1 and P-2 respectively. The minimum and maximum temperatures were expected to increase in the future under both scenarios with more rise expected under RCP 8.5 as compared to RCP 4.5. The temperature varied from 3.1 °C to 3.3 °C under RCP 4.5 and 4.0 °C to 4.3 °C under RCP 8.5.

In a study by Ikram et al. [38], the future increase in precipitation in Pakistan was predicted to be 2–3 mm/day under RCP 4.5 and 3–4 mm/day under RCP 8.5 whereas the temperature was projected to increase from 3–4 °C under RCP 4.5 and 3–8 °C under RCP 8.5. The results of the present study indicated the temperature extremes to be well within the projected limit by Ikram et al. [38], but the precipitation was underestimated. This could be due to the use of different GCM for the climate projections as use of different climate models for the study areas may suggest different results.

4.4. Water Balance in SWAT

The average annual water balance for the baseline period simulated in the SWAT hydrological model is represented in

Table 3. Average annual water balance for Haro River watershed (1976–2005).

Flow	Water Balance Components	Symbology	Computed Values (mm)	Total Flow(mm)
Inflow	Precipitation	PRECIP	1873.8	1873.8
Outflow	Water Yield	WYLD	1149.82
	Deep Aquifer Recharge	DA_RCHG	8.57
	Evapotranspiration	ET	713.8	1872.19
Losses				1.61

. The incoming flow in the watershed was almost equal to the outflows except a minor difference, which was due to the losses occurring in the watershed. Almost 61% of the incoming flow was discharged as water yield, which showed that the major portion of the precipitation was contributing to streamflow. In comparison, the water yield for the Simly Dam watershed in Pakistan was also estimated to be as high as 59% of the total inflows in the watershed [39]. Evapotranspiration accounted for almost 38% of the total precipitation in the watershed.

Figure 5 shows the variation in the water balance components expressed as percentages, computed with respect to precipitation amount in the baseline and three future periods under RCP 4.5 and RCP 8.5.

4.5. Calibration and Validation of SWAT Using SUFI-2 in SWAT-CUP

Table 4. Calibration Results for SWAT model.

Performance Indicator	P-Factor	R-Factor	NSE	R2
Calibration	0.83	0.88	0.80	0.82
Validation	0.79	0.84	0.77	0.77

summarizes the results of calibration and validation of the SWAT within SUFI-2. The calibration and validation results were based on monthly observed and simulated streamflows for the calibration period (1989–1991) and validation period (1996–September 1998). The model calibration results also fulfilled the criteria for P and R factors. The results of the SUFI-2 indicated that 83% of the observed data during calibration and 79% of the data during validation period was bracketed by the 95 PPU band. The criterion for the R-factor was also achieved. According to Abbaspour et al. [30], a P factor value of 0.80–1.00 indicates the least number of conceptual errors in the model as well as a higher quality of input data. NSE value of 0.80 and 0.77 and R² of 0.82 and 0.77 for calibration and validation respectively indicated good model performance.

Figure 6 shows the comparison of observed and simulated flows for Haro River watershed during the calibration and validation period. The simulated discharges during the calibration period were higher compared to the observed flows from July until September 1989 whereas the simulated flows were lower in April 1991 as well as from July until September 1991. Similarly, during the validation period, the simulated flows were much lower than the observed flows from April 1997 until June 1997 and in August 1997. This flow difference was due to the difference in observed precipitation and the climate model precipitation at both stations P-1 and P-2.

The initial and final ranges of the parameters considered for the calibration of discharge in SUFI-2 are given in

Table 5. Description of parameters with the initial and final parameter in SUFI-2.

Sr. No.	Parameters	Description	Initial Parameter Ranges		Final Parameter Ranges
			Minimum Value	Maximum Value	Minimum Value	Maximum Value
1	CN2	SCS Runoff Curve Number	−0.2	0.2	−0.2	0.02
2	ALPHA_BF	Base flow alpha Factor	0	1	0.32	0.98
3	GW_DELAY	Groundwater delay	30	450	30	278
4	GWQMN	Threshold depth of water in the shallow aquifer required for return flow to occur	0	2	0	1.1
5	GW_REVAP	Groundwater revap coefficient	0.02	0.2	0.02	0.19
6	ESCO	Soil evaporation compensation factor	0	1	0	0.58
7	CH_N2	Manning’s “n” value for the main channel	0	0.15	0.07	0.15
8	CH_K2	Effective hydraulic conductivity in main channel alluvium	0	130	0	66
9	SOL_AWC	Available water capacity of the soil layer	0	1	0.3	0.9
10	SOL_K	Saturated hydraulic conductivity	−1	1	−1	0.17
11	SOL_BD	Moist bulk density	−1	1	−0.3	1
12	SFTMP	Snowfall temperature	−5	5	−0.25	5
13	SMFMX	Maximum melt rate for snow during year	1.5	7.5	4.5	7.5
14	SMFMN	Minimum melt rate for snow during the year	1.5	7.5	3.9	7.5
15	SMTMP	Snow melt base temperature.	−5	5	−0.35	5
16	TIMP	Snow pack temperature lag factor	0	1	0.24	0.75

. The absolute values of the parameters identified in SUFI-2 were constrained with each iteration while computing the objective function. The final ranges of the parameters indicated the uncertainty ranges that gave the best value for the objective function i.e., NSE and R².

The results of the global sensitivity analysis in SUFI-2 are shown in Figure 7. Based upon least p-value and the largest absolute t-stat, GW_DELAY followed by CN2, SOL_K and SOL_BD were ranked as the most sensitive parameters. This indicated that parameters related to groundwater processes, surface runoff and soil processes have the highest impact on the results of the SWAT model for this particular watershed. A higher value of CN2 indicated its key role in the generation of streamflow for Haro River watershed as CN is used by SWAT to divide precipitation into overland flow or infiltration. For each HRU, surface runoff was generated within SWAT by using SCS-CN method [40].

4.6. Variation in Streamflows

The monthly variation in streamflows under RCP 4.5 and RCP 8.5 for the baseline and three future periods is shown in Figure 8. Overall, compared to the baseline discharge of 7.7 m³/s at Dhartian gauging station, in the future, the streamflows were expected to increase by 13.0% under RCP 4.5 and 20.8% under RCP 8.5. Under both scenarios, no shift in the peak discharges was anticipated in future. The trend in the streamflows was directly related to the trend in mean monthly precipitation at both stations P-1 and P-2. On the other hand, temperature did not show a major influence on the peak flows in the Haro River watershed in future. The maximum streamflows during the baseline and the future scenarios were expected during the month of July. Under RCP 4.5, the maximum flow was expected during the early century followed by the late century and mid-century whereas under RCP 8.5, the flow was expected to be maximum during the mid-century followed by the early and late century. This was found to be similar to the precipitation trend of Murree.

Figure 9 shows the box and whisker plot for the mean annual streamflows during the baseline and future periods under RCP 4.5 and RCP 8.5. Under RCP 4.5, during the late century the increase in average annual streamflow was highest among all future periods compared to the baseline flow. Under RCP 8.5, the maximum increment in streamflows was expected during early century with respect to the baseline period. The largest variation in streamflows was expected to be during the late century under RCP 4.5.

4.7. Comparison in Future Projections under RCP 4.5 and RCP 8.5

Overall, the results of climate projections under MIROC-ESM suggested an increase in precipitation as well as maximum and minimum temperatures at both stations, P-1 and P-2. The projected increase in future climate variables was higher under RCP 8.5 than under RCP 4.5, compared to the baseline period. The precipitation increase at both stations ranged from 4.5 mm/day to 6.1 mm/day under RCP 4.5 and 4.8 mm/day to 6.4 mm/day under RCP 8.5. Under RCP 4.5, the maximum temperature at both stations was projected to increase by 20.8 °C and 31.9 °C while the increment in minimum temperature was 12.6 °C and 16.9 °C at the end of 21st century. Similarly, under RCP 8.5, the maximum temperature at both stations was projected to increase by 21.7 °C and 32.7 °C while the minimum temperature was projected to increase by 13.7 °C and 17.9 °C. The streamflows in the Haro River watershed were also projected to increase in future with higher streamflows expected under RCP 8.5 (9.3 m³/s) compared to RCP 4.5 (8.7 m³/s). The higher range of climate and hydrologic changes under RCP 8.5 can be attributed to the highest radiative forcing level and concentration level of gases in this particular scenario among all scenarios [37].

5. Conclusions, Limitations and Recommendations

5.1. Conclusions

The study concludes that the linear scaling technique can be successfully applied to remove bias in the temperature and precipitation time series. The statistical indicators, after bias correction, showed considerable improvement compared to the observed time series. SWAT well simulated the hydrology of the data-scarce, high elevation region and was a good tool in evaluating the components of the water balance.

Precipitation was expected to increase in future under both emissions scenarios during early, mid and late century. Summer season was anticipated to receive maximum precipitation whereas winter season was expected to remain relatively dry. The minimum and maximum temperatures were also expected to be highest during early summer season. The increase in minimum temperature in future was more significant as compared to the increase in maximum temperature with respect to the baseline period.

The study indicated that in the future, the streamflows at Haro River watershed were expected to increase with maximum streamflows expected during the summer season. In the future, an increase in streamflows indicate that water scarcity at Khanpur Dam is less likely to be the major concern for policy makers. On the other hand, this increase may enhance erosion and sedimentation and may lead to a decrease in the storage capacity of the Khanpur Dam reservoir. The management strategies at Khanpur Dam may focus on adapting de-siltation plans as well as increasing water storage capacity of the reservoir. The stored water, during the high flow periods could be used to fulfill the demand during low flow period.

The changing streamflows at Haro River basin would affect the discharge at Khanpur Dam, which is a major source of fulfilling the domestic water supply needs of large population. The coordination among the involved institutions, for example, Water and Sanitation Authority (WASA), Capital Development Authority (CDA) and Rawalpindi Cantonment Board (RCB), in conceiving prudent schemes for effective utilization of water supply, is essential to fulfil the demands of all sub-sectors utilizing water in all the periods of the year.

Finally, the potential climate change impact on sedimentation should be studied and sedimentation plans should be proposed. The impact on the sediment yield can be analyzed in SWAT model using long-term suspended sediment data as well as discharge data in the historical period and use it to project future periods. This is an opportunity for continuing research for the Khanpur Dam. The integration of sedimentation plans, in future studies, may help to preserve the storage capacity of reservoir and encourage proper maintenance plans. An insight into devising watershed management plans focusing on delaying the siltation process in the reservoir may help to reduce the effect of climate change on the reservoir. The effect of changing landuse can also be integrated in the SWAT model to study its impact on discharge and sedimentation within the watershed.

5.2. Limitations and Recommendations

The discharge data was available only for seven years in total: 1989–1991 and 1996–1998 which made the calibration and validation challenging. However, the results of the calibrated parameter ranges when applied to an independent validation time series showed good model performance. Steps should be taken by authorities to set up gauging stations in the watershed, as Dhartian gauging station is the only discharging measuring station in Haro River, which is non-functional since 1998. The daily discharge data must be recorded in order to understand the hydrological processes occurring in the watershed in a better way.

The daily data for solar radiation, relative humidity and wind speed data was not available from local meteorological stations. Therefore, this data was extracted from the CFSR dataset for the similar time period at the stations corresponding to P-1 and P-2.

In the study, landuse was assumed constant. In future, an insight into climate change due to variation in landuse and land cover is recommended to be considered. Research in future may also focus on studying the effect of variation in climatic parameters on sedimentation.

Finally, there is inherent uncertainty in the use of one GCM to propagate climate change projections into projections of future flows. Future work should include a full ensemble of GCM model outputs to better capture the uncertainty of the future climate.