Statistical Downscaling and Hydrological Modeling-Based Runoff Simulation in Trans-Boundary Mangla Watershed Pakistan

Abstract

The economy of Pakistan relies on the agricultural sector which mainly depends on the irrigation water generating from the upper Indus river basin. Mangla watershed is a trans-boundary basin which shares borders of India and Pakistan, it comprises five major sub-basins, i.e., Jhelum, Poonch, Kanshi, Neelum and Kunhar. The runoff production of this basin is largely controlled by snowmelt in combination with the winter precipitation in the upper part of the basin and summer monsoon. The present study focusses on the application of a statistical downscaling method to generate future climatic scenarios of climatic trends (temperature and precipitation) in Mangla watershed. Statistical Downscaling Model (SDSM) was applied to downscale the Hadley Centre Coupled Model, version 3, Global Climate Model (HadCM3-GCM) predictions of the A2 and B2 emission scenarios. The Soil and Water Assessment Tool (SWAT) hydrological model was used for the future projected streamflows based on developing climate change scenarios by SDSM. The results revealed an increasing trend of annual maximum temperature (A2) at the rates of 0.4, 0.7 and 1.2 °C for the periods of 2020s, 2050s and 2080s, respectively. However, a consistent decreasing trend of temperature was observed at the high-altitude region. Similarly, the annual minimum temperature exhibited an increasing pattern at the rates of 0.3, 0.5 and 0.9 °C for the periods of 2020s, 2050s and 2080s, respectively. Furthermore, similar increases were observed for annual precipitation at the rates of 6%, 10%, and 19% during 2020, 2050 and 2080, respectively, for the whole watershed. Significant increasing precipitation trends in the future (2080) were observed in Kunhar, Neelum, Poonch and Kanshi sub-basins at the rates of 16%, 11%, 13% and 59%, respectively. Consequently, increased annual streamflow in the future at the rate of 15% was observed attributing to an increased temperature for snow melting in Mangla watershed. The similar increasing streamflow trend is consistent with the seasonal trends in terms of winter (16%), spring (19%) and summer (20%); however, autumn exhibited decreasing trend for all periods.

1. Introduction

Upper Indus River Basin (UIRB) serves as a main source of irrigation water for agricultural sector, which originates from the mountainous areas of Hindukush-Karakoram-Himalaya (HKH) range. Several studies highlighted the impacts of climate change over Upper Indus Basin (UIB) [1,2,3,4,5], for instance, the scientific assessment by the IPCC regarding predicted globally increased temperature between 1.5 and 4.5 °C in the next century with a doubling of the CO₂ concentration in the atmosphere [6]. The global climate change created a negative impact on the precipitation patterns, which ultimately declined the runoff on global and regional scale [7]. Such warming trends not only affect the global circulation patterns which are often simulated and predicted using the General Circulation Models (GCM) but also directly affect the local climatic conditions by altering the distribution and characteristics of precipitation and temperature. Hydrological response resulted from the climate change may have significant impact on existing water resources systems by changing the hydrological cycle [8]. Therefore, the hydrologic extremes like floods and droughts have the negative impacts on all watershed [9]. However, planning and organizing of the new water infrastructure is not developed according to the projected changes of precipitation, temperature and streamflows for an efficient water resource management.

Recently, several researchers have conducted the studies to analyze the impacts of climate impacts on natural resources through climate change and hydrological models and have found it as a complex process [10,11,12,13,14,15,16]. However, the output of GCMs cannot be used directly for local impact studies due to the coarse resolution (250 × 250 km). Many key challenges in the application of GCMs and hydrological models still exist [9]. Therefore, downscaling techniques have emerged as a means of deriving local-scale weather from regional-scale atmospheric predictor variables [17,18,19].

Many downscaling techniques are available in literature for the conversion of GCM outputs into daily meteorological variables. For instance, Mascaro et al. [12] stated that the efficacy of the downscaling techniques used to analyze the future climate variability depends upon the characteristics of the study area and data availability. Two main approaches, i.e., dynamic downscaling and statistical downscaling, are generally used to derive information from global climate scenarios at local or regional scale. Dynamic downscaling (DD) involves a nested regional climate model (RCM) and is suitable in the regions with complex topography and highly heterogeneous land covers [20]. However, statistical downscaling techniques develop a statistical link between the GCM’s and variability derived from the historical time series data on the local level [21,22]. Therefore, compared to other methods, statistical downscaling approach is relatively easy to use and provide station-scale climate information from the GCM scenarios [23].

Three types of statistical downscaling techniques were extensively used in previous studies, i.e., statistical regression methods, weather generators and weather organization systems [24,25]. Therefore, this method of downscaling is widely utilized in estimation of the anticipated effect under various types of study scenarios [26]. The main benefit of statistical downscaling is the low cost, facilitate the local level information under the scenarios of climate and its impacts on the water resources. This method ensembles the scenarios of climate which permit uncertainty. Among the large scale downscaling approach, Statistical Downscaling Model (SDSM) has already been used in many studies round the world because of its universal application in analyzing the climate change impacts [27,28,29,30,31,32,33].

Therefore, also, the performance of the statistical downscaling model is evaluated in simulating the daily temperature and precipitation time series data.

We generate runoff simulations in the Mangla watershed using the Soil and Water Assessment Tool (SWAT), for the period 2001 to 2010. The present study seeks to apply SWAT by employing available observed data collected from five major sub-basins, i.e., Jhelum, Poonch, Kanshi, Neelum and Kunhar as basic input to SWAT. Moreover, the present study aims to investigate the feasibility correlation of maximum and minimum temperature and precipitation at a local scale with large scale atmospheric variables using statistical downscaling method. Moreover, future climate change scenarios of temperature and precipitation on monthly, seasonal and annual basis were developed under outputs of HadCM3, GCM model emission scenarios (A2 and B2), which were utilized for analyzing the influence of climate parameters on a trans-boundary Mangla watershed. The core objective of the work is to understand the effect of temperature and precipitation on the fluctuation of runoff generated in Mangla watershed. The river runoff being highly sensitive to precipitation and temperature in the form of direct (precipitation) and delayed (snowmelt) runoff indicates direct impact of changing climate. We present the combination of a degree-day model SWAT and SDSM for the first time within Mangla watershed. It is expected that our findings will make an insightful contribution to directly inform the development of sustainable water use policies, in the context of a changing climate and improved water resources management.

2. Data and Methodology

2.1. Study Area

The Jhelum river basin is a trans-boundary watershed shared by Pakistan and India. The basin originates from Verinag-Spring situated on the southern slopes of the Himalaya Mountains ranges at the base of Pir-Panjal ranges in Jammu and Kashmir with a drainage area of 33,425 km² at Mangla dam. Around 55% area of this basin lies in Indian held Kashmir while the rest of the area is controlled by the government of Pakistan marked as Azad Kashmir. The total storage capacity of this Mangla dam is 5.8 MAF. The primary purpose of this dam is to irrigate the large agricultural lands of the Punjab and Sindh provinces. However, the hydropower is also generated as a by-product. There are five major sub-basins of Mangla watershed, i.e., Jhelum, Poonch, Kanshi, Neelum, and Kunhar drain into the Mangla reservoir (Figure 1). The tributaries of Jhelum, Neelum and Kunhar rivers join the mainstream at Domel, Muzaffarabad and Kohala Bridge, respectively (Figure 1). The Jhelum, Poonch and Kanshi rivers enter the Mangla reservoir in the Mirpur district. In the southeast, Mangla watershed shares borders between India and Pakistan. The runoff production of this basin mainly dictates by snowmelt and the winter precipitation in the upper part of the basin and monsoon precipitation affects the lower part of the basin.

2.2. Data and Analysis

The datasets used for runoff simulation to analyze the current and future climate change impacts on streamflow in Mangla watershed are shown in

Table 1. Type of data used in the present study and their source.

Data Type	Source	Spatial/Temporal Resolution	Description/Period of Record
Topography	USGS National Elevation Dataset	90 × 90 m	DEM (Elevation)
Land use data	European Space Agency (ESA) Global Land Cover http://ionia1.esrin.esa.int/	300 × 300 m	Classified land use such as forest, agriculture, crops, water, etc.
Soil data	FAO–International Soil Reference and Information Centre (ISRIC)	1 km	Classified soil and physical properties as sand, silt clay bulk density, etc.
Climate data	Pakistan Meteorological Department (PMD), Water and Power Development Authority (WAPDA)	Daily	Precipitation, Temperature, solar radiation, wind speed (1961–2010)
Hydrological data	Water and Power Development Authority (WAPDA)	Daily	Stream flows (1961–2010)
Climate Scenario Data (GCM-Had CM3)	Canadian Climate Scenarios Network (http://www.cics.uvic.ca/scenarios/sdsm/select.cgi)	2.5° latitude 3.75° longitude Daily	1961–2099

.

2.2.1. Climate Change Scenarios (Predictors)

Climate change scenarios were generated from the GCM outputs. These scenarios were developed for the climate (minimum, maximum temperature and precipitation) variables. The outcomes from the Hadley Centre Coupled Model, version 3 (HadCM3)—Global Climate Models (GCM) model for the A2 (high emission scenarios) and B2 (low emission scenarios) have also been utilized to estimate the future scenarios. Predictor variables provide daily information at large-scale obtained from HadCM3 outputs at 2.5° latitude 3.75° longitude gridded scale (http://ccds-dscc.ec.gc.ca/?page=pred-hadcm3).

2.2.2. Climate Station Data (Predictands)

Predictands are the observed climatic variables, i.e., temperature and precipitation. Daily maximum and minimum temperature data of thirteen (13) stations and fifteen (15) stations for precipitation data in the study area were used as predictands as given in

Table 2. Climatic stations in Mangla watershed.

Sr. No.	Station	Lat (degree)	Lon (degree)	Elevation (m)	Catchment	Max. Temp (°C)	Min. Temp (°C)	Precipitation (mm)
1	Murree	33.9	73.4	2206	Jhelum	17.3	8.9	1782
2	Bagh	34.0	73.8	1067	Jhelum	25.1	12.0	1482
3	Mangla	33.1	73.6	282	Jhelum	29.7	17.1	867
4	Gujar Khan	33.3	73.3	457	Kanshi	28.5	14.9	825
5	Kallar	33.4	73.4	518	Kanshi	28.5	16.4	973
6	Balakot	34.6	73.4	995.5	Kunhar	25.1	12.2	1623
7	Naran	34.9	73.7	2363	Kunhar	12.0	2.5	1793
8	Astore	35.2	74.5	2168	Neelum	15.7	4.1	488
9	Garidopatta	34.2	73.6	813.5	Neelum	26.0	12.5	1413
10	Muzaffarabad	34.4	73.5	702	Neelum	27.7	13.6	1484
11	Kotli	33.5	73.9	610	Poonch	28.4	15.3	1219
12	Palandri	33.7	73.7	1402	Poonch	22.9	12.1	1390
13	Rawalakot	34.0	74.0	1677	Poonch	20.6	9.2	1369
14	Khandar	33.5	74.1	1067	Poonch	-	-	1062
15	Sehr kakota	33.7	74.0	1125	Poonch	-	-	1342

. The data from these stations were collected from Surface water Hydrology Project (SWHP), WAPDA and Pakistan Meteorological Department (PMD) for the period 1961–2010.

2.2.3. Hydrological Data

Daily streamflow data for the period of 50 years (1961–2010) at 9 stations was collected from Water and Power Development Authority-Surface Water Hydrology Project (WAPDA-SWHP). The geographical locations of the stream gauge stations along with their respective sub-basins are presented in

Table 3. Stream gauges installed in Mangla watershed.

Sr. No.	Station	Latitude (degree)	Longitude (degree)	Elevation (m)	Basin Area (Km2)	River Name	Mean Annual (Cumec m3/s)
1	Naran	34.9	73.7	2400	1107	Kunhar	46
2	G-Habibullah	34.4	73.4	820	2433	Kunhar	102
3	Muzaffarabad	34.4	73.5	670	7412	Neelum	325
4	Chinari	34.2	73.8	1070	13,652	Jhelum	300
5	Domel	34.4	73.5	701	14,396	Jhelum	328
6	Kohala	34.1	73.5	560	24,464	Jhelum	776
7	Azad Pattan	33.7	73.6	485	25,967	Jhelum	1239
8	Palote	33.2	73.4	400	867	Kanshi	6
9	Kotli	33.5	73.9	530	3210	Poonch	128

. The mean annual streamflow at low-altitude stations was higher as compared to the high-altitude stations. The highest mean annual flow was observed in the months of August and June at low and high-altitude stations, respectively. At the low-altitude station, the source of stream flow is monsoon rainfalls during the months of June to August. Whereas, high-altitude stations affected by winter snow and glacier melt during March to May. Three sub-basins (Jhelum, Poonch, Kanshi) drain into Mangla reservoir with mean annual inflow of 967 m³ s⁻¹. (1961–2010).

2.2.4. Land Use and Soil Data

The land use maps were prepared by the processing of Global Land Cover data to identify the major land use classes in the watershed. The map shows 8 different types of land cover (Figure 2): 45% of the basin is covered by irrigated and rain-fed cropland, 20% by mosaic grassland, 6% by bare soil, and about 4.5% by Water, Snow and Ice. The soil Map was acquired from IRSIC, Global soil classes developed by the UNO agency Food and Agriculture Organization (FAO). Two major classes of soil were identified as the loamy non-calcareous (27%) and rock outcrop with mainly loamy partly gravelly soils (50%) in the watershed. The SWAT uses Agricultural Land-Close-grown (AGLC), Forest-Deciduous (FRSD), Forest-Evergreen (FRSE), Forest-Mixed (FRST), Pasture (PAST), Range-Grasses (RNGE), Urbanization (URBN) and Water as shown in Figure 2a.

2.2.5. Statistical Downscaling Model (SDSM)

This model was applied to estimate downscale the large-scale GCM data at local scale. SDSM was integrated with multiple regression and stochastic downscaling techniques [19,34], previously used in the many studies worldwide [35,36,37]. The model employs an ordinary least squared algorithm for the estimation of parameters of the regression equation. In downscaling the temperature data, a direct link among large-scale predictors and predictand at local level was assumed.

In this technique, most relevant predictor variables were selected by the approach of correlation matrix. The scatter plots were also used between predictors and prediction variables. SDSM model was set to consider the variance inflation, bias correction and transformation functions for screening [38,39] and the predictor variables with acceptable correlation value and probability were selected. Identified predictor variables for the selected stations and predictand variables are presented in

Table 4. Selected predictor variables with their respective predictand.

Station	Precipitation	Maximum Temperature	Minimum Temperature
Kotli	shum, p5_z, p8_f, p850, rhum	p500, temp	p500, temp, shum
Rawalakot	temp, p850, mslps, p5_f, r850	p500, temp	p500, temp, shum
Plandri	shum, p8_f, p850, p5_f, p_z, p5_u	p500, temp	p500, temp, shum
Khandar	temp, shum, p5_z, rhum, mslp, p_v, p_zh, p5_v	-	-
Sehr Kokata	temp, shum, rhum, p5_f, p5_u	-	-
Murree	p5-z, p8-f, r850, rhum, shum	p500, temp	p500, temp, shum
Bagh	p5-f, p5-u, p5th, shum	p500, temp	p500, temp, shum
Gujar Khan	mslp, p5-f, p5th, 9850, shum, temp	p500, temp	p500, temp, shum
Kallar	mslp, p5-u, p5th, r500, shum	p500, temp	p500, temp, shum
Mangla	mslp, p5th, p850 shum	p500, temp	p500, temp, shum
Balakot	p5-f, p5th, p8-f, r850, rhum, shum	p500, temp, shum	p500, temp
Muzaffarabad	mslp, p5-v, p850, r850, shum	p500, temp	p500, temp, shum
Naran	p-v, p500, p5th, shum, temp	p500, temp, shum	p500, temp, shum
Astore	p-f, p-u, p-v, p8-f, rhum, temp	p500, temp, shum	p500, temp, shum
Garidopatta	mslp, p-f, p-u, p-v, p8-f, rhum, temp	p500, temp, shum	p500, temp, shum

Mslp = mean sea level pressure; rhum= near surface relative humidity; shum = surface specific humidity; p5_z = 500 hpa voritcity; p8_f = 850hpa airflow strength; r850 = relative humidity at 850 hpa; r500 = relative humidity at 500 hpa; p5-f = 850 hpa airflow strength; p5-v= 500 hp meriodinal velocity; p5th = 500 hpa wind direction; p-f = surface air flow strength; p-u = surface zonal velocity; p-v = surface meridian velocity; temp= mean temperature at 2 m.

. This calibration and validation of SDSM relied on the outcomes from the second step. The selection of the predictor parameters utilized the base of the National Centre for Environmental protection (NCEP) data for the selected gridded box. After that, links between the specific and selected predictor variables are estimated. From this, multiple linear regression equation is determined.

In this investigation, the calibration was applied for the base period of 30 years (1961–1990) for most of the stations. The monthly based model used to investigate the temporal variation on the monthly basis. From 1991–2001, ten years of simulation were taken from the validation process. Three types of climate data set were utilized through the H3A2 and H3B2 scenarios for validation of SDSM model. Eight evaluation objective functions were applied in the simulation series as average (μ), coefficient of variation (C_v), skewness (C_s), R square (R²), Root Mean Square Error (RMSE), Nash-Sutcliffe Coefficient (Nr), Percent bias (PBIAS) and Mean Absolutely Error (MAE) to qualify the simulation results of temperature and precipitation series in each station as shown in

Table 5. Comparison of determination coefficient (R²), Root Mean Square Error (RMSE), Nash-Sutcliffe Coefficient (Nr), Percent bias (PBIAS) and Mean Absolutely Error (MAE) between observed and simulated results of maximum temperature, minimum temperature and precipitation for each station in the validation period (1991–2000).

Variable	R2		RMSE (mm or °C/Month)		Nr		PBIAS (%)		MAE (mm or °C/Month)
	Range	Mean	Range	Mean	Range	Mean	Range	Mean	Range	Mean
Precipitation
NCEP	0.17–0.6	0.36	46–172	93	−0.1–0.6	0.25	−11.1–29	−7.8	32–121	66
H3A2	0.17–0.4	0.27	47–171	105	−0.8–0.3	0.25	−10.8–26	−4.7	33–125	77
H3B2	0.11–0.5	0.28	46–167	103	−0.6–0.3	0.15	−10.4–14	−10	32–120	72
Max. Temperature
NCEP	0.88–0.97	0.93	1.5–3.4	2.2	0.7–1	0.92	−11–3.5	−2.4	1.2–2.8	1.7
H3A2	0.87–0.96	0.92	1.7–3.9	2.3	0.7–1	0.93	12–3.9	−1.4	1.3–3.2	1.7
H3B2	0.87–0.95	0.91	1.7–3.7	2.4	0.7–1	0.89	11–4.2	−0.7	1.3–3	1.8
Min. Temperature
NCEP	0.90–.98	0.95	0.9–2.6	0.92	0.86–0.98	0.92	40–14	−2.1	0.7–2.1	1.4
H3A2	0.92–0.97	0.94	1.2–2.8	1.9	0.86–0.97	0.93	18–8	−0.1	1–2.2	1.5
H3B2	0.92–0.96	0.94	1.4–2.9	2	0.87–0.96	0.92	15.3–9.1	1	1.1–2.3	1.5

.

It was found that R² values are higher than 0.27 for the precipitation time series on daily basis. These values are high as compared to the previous studies [12,24,40]. Such subject of downscaling is very difficult for the daily precipitation due to stochastic processes. It could be seen that precipitation was underestimated (PBIAS was up to −10) for all three sets of atmospheric data. Results indicated that R square values between simulated and observed temperature for the investigated weather stations increase above 90% while for downscaled precipitation are lower than 40%.

2.3. Hydrological Model

A semi-distributed (SWAT) model was used to simulate the stream flows of Mangla basin. This system application is integrated with the geographic information system (GIS), which consists of the various spatial data sets included soil type, land use pattern and topographic features. This model is efficient and uses the easily available datasets, enabling users to investigate the long run impacts. Furthermore, it is the physically based, long term series and efficient hydrological model [41]. This model is categorized as semi-distributed model because it divides sub-watersheds into further small units called hydrological response units (HRUs). The hydrologic response unit (HRUs) consists of homogeneous land cover, soil characteristics and management practices. Further subdivisions of sub-basins into HRUs enable to account for the impact of different land cover types, management practices and soil characteristics on the hydrologic cycle of any basin. The model generates the daily datasets of the climate variables using the arithmetical techniques which are derivative from the average monthly datasets. The basic input data for SWAT model run include the digital elevation model (DEM), soil map, soil properties, land use cover map, climatological variables (precipitation and temperature) and hydrological (streamflow) data.

2.4. Model Evaluation

For each station, the reliability of the model is assessed using different model performance criteria. The criteria used include the determination coefficient (R²), Nash-Sutcliffe Coefficient (N_r) and Percent bias (PBIAS). Analysis is common, and reference is made to applications in [15]. Each of the criteria is briefly described below:

• Nash-Sutcliffe Coefficient (N_r)

  $$N_r=1-\frac{{\displaystyle {\sum }_{i=1}^N{\left(S_i-O_i\right)}^2}}{{\displaystyle {\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2}}$$

  (1)
• Percent bias (PBIAS)

  $$PBIAS=\frac{{\displaystyle {\sum }_{i=1}^N\left(S_i-O_i\right)}}{\overline{O}}$$

  (2)

  where N is the number of observations; S_i and O_i are the simulated and corresponding observed values at time i; Ō is the average observation.

2.5. Time Scales

Three time-series were selected, i.e., 2020s (2011–2040), 2050s (2041–2070) and 2080s (2071–2099) to present the simulated climatic projections. The future climate change scenarios were also developed at monthly, seasonal and annual scales. Different seasons of Mangla watershed are defined as follows: Winter (December–February), spring (March–May), summer (June–August) and autumn (September–November). Annual mean is the average of January to December monthly means.

3. Results

3.1. Model Simulation

SWAT was simulated by using the precipitation and temperature data for the period of 1 January 1996 to 31 December 2010. During calibration procedure, some other parameters also need to be adjusted to achieve GOF, which are given in

Table 6. Model calibration and verification statistics for daily and monthly stream flow comparison.

Objective Function	Rivers
	Jhelum			Poonch			Kanshi
	Cal-P1	Cal-P2	Val	Cal-P1	Cal-P2	Val	Cal-P1	Cal-P2	Val
Daily Simulation
R2	0.836	0.819	0.752	0.883	0.810	0.710	0.804	0.826	0.697
NSE	0.83	0.79	0.62	0.837	0.77	0.68	0.85	0.81	0.74
PBIAS	3.8	4.6	7.3	9.1	−2.1	5.8	7.2	4.5	−8.2
Monthly Simulation
R2	0.893	0.871	0.826	0.929	0.941	0.902	0.963	0.893	0.894
NSE	0.69	0.823	0.729	0.857	0.871	0.93	0.95	0.89	0.786
PBIAS	−0.93	5.07	−5.6	10.12	0.05	7.46	5.54	3.64	5.13

Cal-P1 = Calibration Period 1st (1996–2010), Cal-P2 = Calibration Period 2nd (2006–2010). Val = Validation for the period (2001–2005).

. The observed and simulated values were plotted against each other in order to determine the goodness-of-fit (GOF) criterion of coefficient of determination (R²) for runoff. The (R²) values of rivers Jhelum, Poonch and Kanshi for daily time span 2001–2005 were 0.836, 0.715 and 0.804, and for 2006–2010 were 0.819, 0.789 and 0.826 as given in Table 6. Figure 3 shows that almost all the parameter values produce similar model efficiency. The main purpose of these plots was to demonstrate the parameter sensitivity as well as distribution of sampling points. The daily observed and simulated hydrographs during calibration periods at Azad Pattan, Kotli and Palote stations for Jhelum, Poonch and Kanshi rivers, respectively are shown in Figure 4.

The model was validated by comparing simulated and observed runoff hydrographs at the Azad Pattan, Kotli and Palote stations for Jhelum, Poonch and Kanshi rivers, respectively for five years as shown in Figure 5. The validation was carried out for the years 2001 to 2005. It was observed that the model simulated streamflows were consistent with the observed streamflows. The values of parameters for precipitation lapse rate (PLAPS) and Surface runoff lag time (SURLAG) at a given range of the model. Several simulation trials were run to calibrate the model successfully during calibration procedure, some other parameters also need to be adjusted to achieve GOF, which are given in

Table 7. List of calibrated parameters and their optimized value.

Parameter	Description	Range	Optimized Value of Parameters for Rivers
			Jhelum			Poonch			Kanshi
			Cal-P1	Cal-P2	Average	Cal-P1	Cal-P2	Average	Cal-P1	Cal-P2	Average
TLAPS	Temperature lapse rate (°C/km)	0−10	−3.1	−4.6	−3.85	−5.7	−6	−5.85	−7	−6.3	−6.65
PLAPS	Precipitation lapse rate (mm H2O/km)	0,100	13	17	15	9	11	10	8	9	8.5
SFTMP	Snowfall temperature (°C)	−5 + 5	3	3	3	2	1	1.5	−	−	−
SMTMP	Snow melt base temperature (°C)	−5 + 5	2.1	1.7	1.9	3.1	3.5	3.3	−	−	−
SNOEB	Initial snow water content in elevation band (mm)	0,300	200	190	195	150	150	150	−	−	−
TIMP	Snow pack temperature lag factor	0.01,1	1	1	1	1	1	1	−	−	−
SMFMN	Melt factor for snow (mm H2O/°C-day)	0,10	3	4	3.5	4.5	4.5	4.5	−	−	−
SMFMX	Melt factor for snow (mm H2O/°C-day)	0,10	4.2	4	4.1	5	5.3	5.15	−	−	−
SURLAG	Surface runoff lag time (hrs)	0,96	31.3	28.5	29.9	17.2	15.8	16.5	9	7.8	8.4
CN2	Runoff Curve Number	54–85	67	70	68.5	69	70	69.5	65	65	65
SOL K	Saturated hydraulic conductivity	0.0 to 2000						0			0
ALPHA BF	Base flow alpha factor	0.00 to 1.00	0.0143	0.014	0.01415	0.012	0.012	0.012	0.014	0.014	0.014
SOL AWC	Available water capacity	0.0 to 1.00						0			0
ESCO	Soil evaporation compensation factor	0.7–1.0	0.95	0.8	0.875	0.85	0.85	0.85	0.7	0.8	0.75
EPCO	Plant uptake compensation factor	0.4–0.9	0.7	0.7	0.7	−	−	−	0.7	0.7	0.7

.

3.2. Generation of Future Climate Change Scenarios

In this study, the duration from 1961–1990 was considered as a base period. The similar time period was used in various studies at a global and regional scale. The future scenario was arranged into three intervals as 2020s (2011–2040), 2050s (2041–2070) and 2080s (2071–2099). From the obtained data from H3A2 and H3B2. The climate variability was assessed using the above-mentioned base period for the entire Mangla catchment.

3.2.1. Downscaling Daily Temperature under Future Emission Scenarios

It was observed that under scenario A2, the annual maximum temperature would increase for the duration of the 2020s, 2050s and 2080s in the (A2) emission scenarios as shown in Figure 5. The results proposed an increasing temperature pattern at the rates of 0.4, 0.7 and 1.2 °C during three future time intervals, respectively. Some of the high-altitude sub-basins exhibited decreasing temperature; however, an increasing trend was observed at low elevated basins. Seasonal Tmax would be increased in all seasons by different magnitude in the whole basin. All sub-basins have the decreasing trend of winter Tmax under both scenarios as shown in Figure 6. The decreasing trends in Summer Tmax were also observed for the Kunhar, Neelum and Poonch basins. Highest positive changes were seen for the autumn Tmax for all sub-basins under both scenarios.

Figure 7 shows the absolute changes of Tmin annually and seasonally for the periods of the 2020s, 2050s and 2080s in the A2 and B2 emission scenarios for Mangla basin and its sub-basin. The changes in Tmin (scenario A2) in the Mangla sub-basin were observed as 3%, 4% and 8%, however, for Kunhar basin, such changes were expected as 7%, 10% and 16% during the 2020s, 2050s and 2080s, respectively. The Neelum sub-basin will exhibit the changes as 5%, 4% and 10%, respectively during these three future intervals. As far as low altitude eastern sub-basins Poonch and Kanshi are concerned, changes would be expected as −2%, 3%, 5% and 3%, 5% and 7%, respectively for the same time intervals. Tmin during winter and spring seasons showed decreasing tendency in the 2020s and 2080s under both future scenarios; however, increasing trends in 2050s were observed. Spring season exhibited the highest warming trend for all three-time intervals.

3.2.2. Downscaling Daily Precipitation under Future Emission Scenarios

The changes in precipitation based on H3A2 and H3B2 datasets were compared to the base period. The results revealed that the precipitation would be increased annually during the periods of 2020, 2050 and 2080 with the values of 6%, 10%, and 19% in the Mangla basin under the H3A2 scenario, respectively (Figure 8). Similarly, the precipitation (A2) would be increased annually in Kunhar, Neelum, Poonch and Kanshi at the rates of 16%, 11%, 13% and 59%, respectively, for the period of 2080s. For B2 scenario, the precipitation changes in Mangla and Kunhar sub-basins were observed as 6, 11%, 16% and 4%, 9%, 15% during three future intervals; however, Neelum sub-basin exhibited precipitation changes as 3%, 5% and 10%. Similarly, the low altitude sub-basins Poonch and Kanshi showed −2%, 1% and 2% and 17%, 39% and 51%, respectively. The changes in precipitation seasonally in Mangla basin and its sub-basin in the scenarios of A2 and B2 presented a complex behavior in different seasons. We found decreasing trends of precipitation in the winter and spring, while summer and autumn precipitation would be increased with different magnitude under both scenarios in Mangla basin. We noted the highest increase during the summer up to 52% and 47% for A2 and B2, respectively, in the 2080s.

3.3. Impact of Climate Change on Water Availability in the Future

The accessibility of water highly depends on the quantity of precipitation and evapotranspiration over Mangla catchment. Therefore, variations in the precipitation and temperature directly affect the soil moisture patterns. A2 and B2 scenarios were utilized to analyze the effect of climate on the natural resources management, particularly water. SWAT model was calibrated for a period of 10 years (1 January 1996 to 31 December 2000) and 1 January 2006 to 31 December 2010). After calibration, daily precipitation, minimum and maximum temperature data of future climate change scenarios was used in SWAT model for three periods, 2011–2040, 2041–2070 and 2071–2099. SWAT was also simulated with the downscaled climate variables of temperature and precipitation datasets. However, wind speed, solar radiation and relative humidity along with the changes in the land use pattern were considered constant throughout the future intervals to assess the future water availability. These simulation results are shown in Figure 9 and Figure 10.

Annul flows for the whole Mangle reservoir would be increased by 15% and 9.2% in the investigated scenarios A2 and B2, respectively, during the period of current and until the 2080s. The observed changes in seasonal streamflow also indicated the same type of pattern and magnitude that vary with respect to location as shown in the Figure 10. From the detailed observations, the future scenario has an effect on the hydrological system of the Upper Indus Basins, which increased the 15% streamflow annually. A similar increasing trend is observed for all seasons except the autumn season. In the 2080s, streamflows would decrease in the autumn season up to 17%. The monthly streamflows may have changes from −44% to 44% for high emission scenarios A2. It was pointed out that 44% is a very high value. Similar results have been reported by different researchers such as −6.1% to 89.9% depending on the RCM in Lake Winnipeg watershed (LWW), Canada [42]; the Western region of Saudi Arabia has an increase of average daily runoff of 353%, and each of the remaining regions has more than a 100% increase [43]; annual projected streamflows would be decreasing up to 61.3%, 49.8%, 56.6% and 44.6% in the Rock River (RRW), Illinois River (IRW), Kaskaskia River (KRW) and Wabash River (WRW) watersheds in the Midwestern United States, respectively, by the end of this century [44]. The temperature increases of 1 to 4 °C have resulted in an increase of about 16.6%, 33.08%, 49.66% and 68.57%, respectively, in current mean annual total flow volume of Hunza basin [45].

4. Discussion

The present study focusses on the application of a degree-day model and statistical downscaling method to generate future climatic scenarios of climatic trends (temperature and precipitation) in Mangla watershed. Statistical Downscaling Model (SDSM) was applied to downscale the HadCM3-GCM predictions of the A2 and B2 emission scenarios. The SWAT hydrological model was used for the future projected stream flows based on developing climate change scenarios by SDSM. The hydrological model shows an early spring snow melt effect, which is attributed to the increasing mean and peak flow in March and April and a corresponding decrease in August. All simulations also showed an increase in the winter low flows, especially in the months of November and December, while some of the simulations show a decrease in low flow in the summer and early autumn. However, the models simulated based on the SDSM downscaled data resulted in an increase in mean annual flow of about 8%, 7% and 41% for Rivers Jhelum, Poonch and Kanshi during the current decade and the 2080s under the scenarios of HadCM3-A2, respectively. The annual streamflow would be 5%, 1% and 36% for Rivers Jhelum, Poonch and Kanshi during the current decade and the 2080s in HadCM3-B2, respectively.

Some studies [46] suggested growing annual Snow-covered area (SCA) trend in the monsoon influenced sub-basins of Mangla basin such as Jehlum basin, attributed to decreased liquid precipitation with increased solid precipitation. Previous studies concluded a reduced SCA during winter and autumn but increasing summer SCA within Hunza [47], and Gilgit basins [46], during 2001–2012, which is still expanding up to 2018 [48]. The growing summer SCA is associated with winter warming and summer cooling effects [49]. Some studies associated the expansion of SCA with increasing precipitation trends [48] and significantly reduced temperature [1], at high altitude stations in Gilgit River basin. Our results partially validated the expansion of SCA attributed to negative temperature trends at high altitude observed in the present study. These findings are also consistent with significant decreases in maximum temperature within the Gilgit basin at Gilgit, Gupis, Ushkore and Yasin [1]. The decreasing temperature in summer decreased the melting rate for new snow as well as for permanent snow and glaciers. The recent temperature and precipitation trends in the Gilgit basin were comprehensively discussed by [1,2]. They found that both annual and seasonal temperatures had significantly decreased at Gilgit, Gupis, Ushkore, and Yasin over the period 1991–2013.

We found that SWAT performs excellently in terms of Nr, R² and PBIAS values. Hydrological model efficiency, particularly degree-day model’s efficiency, can be improved significantly by using field-based TLAPS, PLAPS, snow melt base temperature and snow melt factor [48]. Our findings are inconsistent with the results of [48,50], as they suggested less model efficiency of degree-day model (SRM) over the Hunza and Gilgit River Basins during (July–August). SWAT simulated peak flows excellently during the whole ablation period including peak ablation (July–August). We also observed that SWAT seems more sensitive towards PLAPS and TLAPS as compared to other parameters during calibration. Specifically, the TLAPS is the most subtle parameter; previous studies suggest −0.68 °C/m for the Karakoram Range. However, our updated analysis proposed average TLAPS ranging from −0.35 to −0.55 °C/m for the Mangla watershed.

Some previous studies claimed that summer flow increases at the rates of 33% and 16% increase in Hunza and Shyok rivers, respectively when temperature raised to 1 °C [49]. These authors also proposed that climate variability behaved differently on streamflows for different basins based on selected methods, hypothesis and basin characteristics. Similarly, some studies suggested increasing flows since 1986–2010 and 2001−2012 attributed to increased temperature based on similar hydrological modelling approach (degree-day models) UIB and neighboring basins [48,50,51]. The most recent cryospheric modelled-based simulation study [48] proposed that the flows were not increasing according to the temperature increase as there would be less snowmelt and glacier melt due to intensified warming, so overall water availability in terms of cryospheric extent will be decreased during 2070–2100. Our findings are partially consistent with the above-mentioned results as increased future flows are consistent with the increased temperature; however, increased future flows were also observed in all future time intervals. They associated reduced flows to reduced glacier area in the future, but Mangla basin being a snow-dominant region would be able to provide substantial melt water due to seasonal snowmelt.

Climate change impacts on the hydrologic response of the UIB has also been observed, using station data and model-based studies [1,2,46,47,48]. However, a number of limitations have been highlighted, associated with precipitation underestimation and bias calibration methods, which have curtailed the development of policy to inform more sustainable water resources management in the region [48]. Studies have also highlighted the lack of spatiotemporal data coverage and quality of datasets employed, which impact the variability in simulations and uncertainty in estimating snowmelt contribution to runoff. Previous hydrological modelling studies within UIB [46,49,50,51], assessing the impact of climate variability on projected future flows, found significant biases in runoff contribution, which they attributed to the selection of large catchment areas, poor quality/sparse precipitation data and choice of DDF and TLR [52]. More recent studies have sought to address these shortcomings by employing more specific DDF and TLR applied at sub-basin scales [48].

5. Conclusions

The present study concluded that the statistical downscaling approach combined with degree-day hydrological model (SWAT) can be successfully used at the sub-basin scale for trans-boundary Mangla watershed. The SDSM method was used for the current and future simulation using SWAT. Our downscaled climatological data provided an excellent model performance during all time intervals. The specific conclusions from this study are the following.

The annual maximum temperature would increase for the 2020s, 2050s and 2080s under the A2 emission scenarios. We observed prominent changes in terms of temperature and precipitation throughout the 21st century. Warming trends were observed most noticeably at lower regions, while higher altitude regions exhibited less warming intensity. Similarly, lower regions will become drier in the future; however, higher regions will exhibit wetting trends (positive precipitation trends). We found that annual maximum temperature will increase from −1.5 to 0.83 °C decade⁻¹, while annual minimum temperature is expected to increase from −0.74 to 0.52 °C decade⁻¹. The annual precipitation will decrease 110 to 50 mm decade⁻¹; however, annual streamflow will increase from −66 mm to 17 mm decade⁻¹.

The increased future summer flows could be helpful for hydropower and irrigation. However, water security remains questionable, depending mainly upon the seasonal snow accumulation and subsequent temperature for melting in Mangla watershed. Moreover, Pakistan’s available water storage capacity will be unable to meet such increased flows in the future. The excessive water will be a disaster rather than a beneficial source by damaging standing crops, resulting in a severe drought in terms of food and water scarcity. Water management facilities in Pakistan are still unable to store surplus water during the peak-water/flooding season to serve during drought/low-flow periods.