*Article* **Hydrological Response of the Kunhar River Basin in Pakistan to Climate Change and Anthropogenic Impacts on Runoff Characteristics**

**Muhammad Saifullah 1,\*, Muhammad Adnan 2, Muhammad Zaman 3, Andrzej Wał ˛ega 4, Shiyin Liu 2,\*, Muhammad Imran Khan 5, Alexandre S. Gagnon <sup>6</sup> and Sher Muhammad <sup>7</sup>**


**Abstract:** Pakistan is amongst the most water-stressed countries in the world, with changes in the frequency of extreme events, notably droughts, under climate change expected to further increase water scarcity. This study examines the impacts of climate change and anthropogenic activities on the runoff of the Kunhar River Basin (KRB) in Pakistan. The Mann Kendall (MK) test detected statistically significant increasing trends in both precipitation and evapotranspiration during the period 1971–2010 over the basin, but with the lack of a statistically significant trend in runoff over the same time-period. Then, a change-point analysis identified changes in the temporal behavior of the annual runoff time series in 1996. Hence, the time series was divided into two time periods, i.e., prior to and after that change: 1971–1996 and 1997–2010, respectively. For the time-period prior to the change point, the analysis revealed a statistically significant increasing trend in precipitation, which is also reflected in the runoff time series, and a decreasing trend in evapotranspiration, albeit lacking statistical significance, was observed. After 1996, however, increasing trends in precipitation and runoff were detected, but the former lacked statistical significance, while no trend in evapotranspiration was noted. Through a hydrological modelling approach reconstructing the natural runoff of the KRB, a 16.1 m3/s (or 15.3%) reduction in the mean flow in the KRB was simulated for the period 1997– 2010 in comparison to the period 1971–1996. The trend analyses and modeling study suggest the importance of anthropogenic activities on the variability of runoff over KRB since 1996. The changes in streamflow caused by irrigation, urbanization, and recreational activities, in addition to climate change, have influenced the regional water resources, and there is consequently an urgent need to adapt existing practices for the water requirements of the domestic, agricultural and energy sector to continue being met in the future.

**Keywords:** climate change; Kunhar River Basin; streamflow; trend analysis; Soil and Water Assessment Tool (SWAT); anthropogenic impacts

**Citation:** Saifullah, M.; Adnan, M.; Zaman, M.; Wał ˛ega, A.; Liu, S.; Khan, M.I.; Gagnon, A.S.; Muhammad, S. Hydrological Response of the Kunhar River Basin in Pakistan to Climate Change and Anthropogenic Impacts on Runoff Characteristics. *Water* **2021**, *13*, 3163. https://doi.org/10.3390/ w13223163

Academic Editor: Aizhong Ye

Received: 22 September 2021 Accepted: 2 November 2021 Published: 9 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Human development has been increasing at a substantial rate since the Industrial Revolution [1], at the expense of anthropogenic activities causing land use changes, an increase in emissions of greenhouse gases (GHGs), and therefore climate change [2]. Both changes in climate and anthropogenic activities can affect hydrological processes of a river catchment [3,4], thus impacting on water resources, hydropower production [5,6], and crop yield [7]. According to the Intergovernmental Panel on Climate Change (IPCC) [8], GHG emission have caused a 0.85 ◦C increase in global mean temperature with associated changes in precipitation over the period from 1880 to 2012, thus impacting on the hydrological characteristics of the catchment [9], notably the volume of runoff [10] and peaks value [11], in addition to other runoff characteristics [12].

The identification of trends in precipitation and evapotranspiration time series (or more commonly temperature as a proxy for evapotranspiration) can help us understand the complex temporal variability of streamflow. Canchala et al. [10], for example, examined the variability of river flow of two Colombian rivers, which they corelated with various indices of atmospheric teleconnections as well as precipitation. Similarly, a trend analysis was performed on the runoff of the Athabasca River Basin in western Canada [11]. The author found that the decreasing trend in streamflow in recent decades was coherent with the temperature and precipitation trends, and that the trends in hydrological variables that the catchment have recently experienced are projected to continue under climate change.

The generation of runoff and its characteristics are not only affected by climatic changes, but also by anthropogenic activities [13]. For this reason, land-use/land-cover (LULC) types and changes in the latter are important to monitor, as the hydrological response of forested land, urbanized land [14], and cultivated land [15,16] to a precipitation event does vary. Applying a catchment-based approach with regard to managing the hydrology of a river basin, including water resources availability, for the benefits of all users is thus necessary [17].

There is a growing body of research on simulating the hydrological response and estimating changes in the characteristics of streamflow to the impacts of climate change and anthropogenic activities. The impacts of anthropogenic activities on streamflow were examined by [10,11,18,19], for instance, while other studies have focused on assessing the impacts of climate change on the hydrological response [20–24]. However, there is a limited number of studies that investigated the hydrological response to the impacts of both climate change and anthropogenic activities in river basin [2,9,25–30], those that have been performed to date have used various methods, including hydrological modeling [31,32], statistical techniques [33,34], empirical methods [35,36] and paired catchment techniques [37,38] and paired years methods [39,40], with each method having its limitations. For instance, in the case of hydrological modeling, the models require long-term and detailed data for their calibration, which are not always available. Moreover, the calibration of the model with limited input data can causes uncertainties and discrepancies in the model outputs. Empirical methods for their part, provide physical interpretation and have fewer data requirements, hence they are usually preferred in developing countries where there are fewer resources to monitor the hydrological conditions.

Many catchments have their source in mountainous regions and hence these regions play an important role in the water balance of areas located further downstream. The majority of inhabitants of the Kunhar River Basin (KRB), located in the western Himalaya of Pakistan, for instance, rely on agriculture for their livelihoods, with the water required for the latter originating from upstream mountain areas, and used for irrigation, power generation, recreation, and municipal use [41]. The region also has great potential for tourism and for further development in hydropower generation and agriculture, which may put additional pressure on water supply in the future.

The KRB has experienced changes in LULC types in recent decades. Woods burning to meet the domestic energy needs and a growing population have caused deforestation in the uphill areas of the basin [42], with forests being replaced by grasses and shrub [43]. Moreover, slope agriculture is becoming more and more common, causing further changes to the vegetation type over the basin [44]. A study conducted by Saifullah et.al. [45] determined the threshold levels and a climate-sensitivity model for the hydrological regime of the KRB. The current study is directly associated with the water-related ecosystem that is designed in accordance with the Sustainable Development Goals of the United Nations (SDGs 6.3), to integrate climate change measures and policies (SDGs 13.2), and anthropogenic impacts (SDGs 15.1.1).

The objectives of this study are: (1) To examine the temporal variability and trends in hydro-climatic variables over the KRB (2) to quantify of the relative contribution of climate change and anthropogenic activities on the variability of runoff in the KRB.

#### **2. Materials and Methods**

#### *2.1. Study Area*

This study was conducted in the KRB, a high-altitude catchment in Northern Pakistan (Figure 1). The Kunhar River is 171 km long; it originates in lake Lulusar in the Kaghan Valley of Khyber Pakhtunkhwa passing through the town of Jalkhand, Bata Kundi, Naran, Kaghan, Kwai, Balakot, and Garhi Habibullah and exiting into the Jhelum River at Rara [22]. The Kunhar River is an important source of water for the Mangla reservoir, which contributes nearly 11% of its water [46]. The drainage area of KRB is approximately 2600 km2, it is mountainous [44], with elevation ranging from 672 to 5192 m above sea level [47].

**Figure 1.** Geographical location of the KRB in Pakistan.

Mangla is the second-largest reservoir in Pakistan and its storage is used to irrigate nearly six million hectares of the country's agricultural land, in addition to producing nearly 1000 MW of hydroelectricity [41,48]. The region experiences mild summers and cold winters. Average annual maximum temperature at Naran, Balakot, and Muzaffarabad is 12.3, 24.9, and 28.4 ◦C, respectively, whereas the average annual minimum temperature is 3.2, 12.4, and 13.5 ◦C, respectively [22]. Annual rainfall at Muzaffarabad and Balakot is on average, 1351 and 1531 mm, respectively [43]. The basin also receives a substantial amount of precipitation in the form of snowfall in the winter (December-March) [44].

The vegetation of the KRB is divers, consisting of coniferous and broadleaf forests, grasses and shrubs, in addition to bare soil and snow and glaciers, and water bodies (Figure 2, Table 1). Three dominant land cover types dominate the basin: grasses/shrubs, bare soil/rocks, and dense coniferous forests, with crops covering only 4.35% of the basin (Table 1).

**Figure 2.** Land-cover classes of the KRB.



#### *2.2. Datasets and Pre-Processing*

A Digital Elevation Model (DEM), with a 30 m resolution, was downloaded from the website of the United States Geological Survey (USGS) [49] and used to delineate the boundaries of the KRB. The land-cover data over the basin, for their part, were downloaded from the website of the International Center for Integrated Mountain Development (ICIMOD) [50], while the Harmonized World Soil Database v 1.2 was acquired from the website of Food and Agriculture Organization (FAO) [51]. There are two gauging stations and three meteorological stations in the KRB (Table 2). The Water and Power Development Authority (WAPDA) provided daily streamflow data from the two gauging stations, while daily temperature and precipitation data were obtained from the Pakistan Meteorological Department (PMD). Wind speed, solar radiation, and relative humidity were extracted from the ERA5 reanalysis dataset at an hourly time scale [52]. The analyses were performed at the monthly and annual time scales, but some methods required daily data, e.g., the hydrological model, as described below.

**Table 2.** List of hydro-meteorological data of the KRB.


Note: WAPDA refer to the Water and Power Development Authority and PMD to Pakistan Meteorological department.

#### *2.3. Methods*

2.3.1. Mann Kendall Trend Test

Mann-Kendall (MK) trend test [53,54], calculated using Equations (1) and (2), is a commonly used non-parametric test to detect trends in a climatological and hydrological time series.

$$S = \sum\_{i=1}^{n-1} \sum\_{j=i+1}^{n} s \mathfrak{g} n \left(\mathbf{x}\_{j} - \mathbf{x}\_{i}\right) \tag{1}$$

$$\text{sgn}(\mathbf{x}\_{j} - \mathbf{x}\_{i}) = \begin{cases} +1, & \mathbf{x}\_{j} > \mathbf{x}\_{i} \\ \mathbf{0}, & \mathbf{x}\_{j} = \mathbf{x}\_{i} \\ -1, & \mathbf{x}\_{j} < \mathbf{x}\_{i} \end{cases} \tag{2}$$

where *xi* and *xj* denote the data values at times *i* and *j*, respectively, *n* specifies the length of the data set. A positive value of *S* indicates an increasing trend, while a negative value refers to a decreasing trend. The variance of *S*, *Var(S*), is calculated using Equation (3), assuming a normally distributed time series with *n* > 10 and.

$$Var(S) = \frac{n(n-1)(2n+5) - \sum\_{i=1}^{n} t\_i i(i-1)(2i+5)}{18} \tag{3}$$

where *ti* denotes the number of data ties. The test statistic, *Z*, is then calculated using Equation (4).

$$z = \begin{cases} \frac{S - 1}{\sqrt{var(S)}} & S > 0\\ 0 & S = 0\\ \frac{S - 1}{\sqrt{var(S)}} & S < 0 \end{cases} \tag{4}$$

In a two-tailed test, the standard Z value is compared with the standard normal distribution table at a 5% significance level (α). The null hypothesis (*Ho*) is rejected if |Z| > |Z1-α/2|, meaning that the trend is statistically significant, otherwise *Ho* is accepted, i.e., there is no presence of a statistically significant trend in the time series at the 95% confidence level.

#### 2.3.2. Innovative Trend Analysis

The Innovative Trend Analysis (ITA) was proposed by Sen [55]. The technique consists of the dividing a time series into two subsets. Hence, the observed hydrological time series was divided into two sub-series after having re-arranged the time series in ascending order. Based on the Cartesian coordinate system, the first sub-series (*Xi*) is drawn on X-axis while the second (*Xj*) is on Y-axis as shown in Figure 3. Generally, the hydrological time series is described as trendless if it is located on the 1:1 (45◦) straight line whereas it is a decreasing trend in time series if the data values are found below the triangular area of the straight line (45◦), and an increasing trend if it is above it.

**Figure 3.** Interpretation of the ITA technique.

#### *2.4. Change Point Analysis*

A change-point analysis, as described in [11], was conducted on the time series with the purpose of determining the presence of abrupt changes in time series. Test determines the number of change points and estimate their time occurrence. This study used the Combinations of the Sum Boxes (CUSUM) plot to determine the change point in a time series. This procedure is described in [56] procedure. The accumulated sum of data points i.e., X1, X2,... ,X24 is calculated and CUSUM graphs are produced. Moreover, in general, the average value displays periods that are greater than the average. Mostly, the values are below the average and the downward slope segment is displayed by the sudden changes in CUSUM direction, which points out an average or sudden change. The period of the CUSUM column is straightforward when the average or sudden change does not occur.

#### *2.5. Double Mass Curve Analysis*

The consistency of of hydrological time series is checked by comparing data for a single station using double mass curve analysis [57]. It can be used to correct the unstable precipitation data. It is basically a fraction of the accumulated image arithmetic figures in a variable different from the accumulated figures of another variable or the same variable occasionally. The relative mass ratios of these variables change concerning the interrelations between the variables. This may be due to changes in the data collection procedure or might be the physical changes that mark the relationship. It must not be thought that all discrepancies shown by a double mass curve were inconsistent due to changes in data collection methods or errors in data collection.

#### *2.6. Flow Duration Curve*

Flow Duration Curve (FDC) [58] is a graphical representation of the overall historical variability in river flow. FDC is a corresponding distribution function of daily flows. The probability that a certain value will be exceeded in a predefined future period is called the probability of exceedance. It is used to predict extreme events such as floods, hurricanes, and earthquakes. The exceedance probability *P* can be calculated using Equation (5).

$$P = \left[\frac{M}{(n+1)}\right] \times 100\tag{5}$$

where *M* specifies the ranked position on the listing (dimensionless) and *n* denotes the number of events during the data record. The exceedance predicts the probability of a given streamflow [11]. Further detail about FDC can be found in a study conducted by [11].

#### *2.7. Eco-Hydrological Framework*

This framework focuses on studying the interactions between water and ecological systems. The annual hydrological budget for a catchment can be expressed by Equation (6):

$$P = ET + Q + D + \Delta S \tag{6}$$

where *P* stands for precipitation, *ET* for the evapotranspiration, *Q* for streamflow, *D* for deep groundwater losses, and Δ*S* for the changes in storage. Deep groundwater losses and changes in water storage can reasonably be assumed to be zero over a long time period. The available water and energy in an agricultural watershed can be assessed through the excess water (*PEx*) and excess energy (*EEx*), which can be estimated using the following equations:

$$P\_{Ex} = \frac{(P\_{\parallel} - ET)}{P} \tag{7}$$

$$E\_{Ex} = \frac{(PET - ET)}{PET} \tag{8}$$

where *PET* denotes potential evapotranspiration. The value of *PEx* and *EEx* can range from 0 to 1. *PET* denotes potential evapotranspiration. The eco-hydrological analysis represented by Figure 4 is an example of the conceptual model that is applicable for understanding the soil water conservations measures and climatic variability impacts on watershed hydrology [59].

#### *2.8. Climate Elasticity Model*

Climate change impacts can also be assessed using the climate elasticity model [35,60,61]. Schaake [60] presented the concept of climate elasticity to assess the sensitivity of streamflow to climate changes. It can be defined as the relative change in streamflow divided by the relative change in a climate variable, precipitation for instance. Schaake [60] defined the precipitation elasticity, *ε <sup>p</sup>*, through Equation (9):

$$
\varepsilon\_p(P, Q) = \frac{\frac{dQ}{Q}}{\frac{dp}{P}} = \frac{dQ}{dP}\frac{P}{Q} \tag{9}
$$

where *P* and *Q* represent precipitation and streamflow, respectively, and *ε <sup>p</sup>* denotes the precipitation elasticity. The precipitation elasticity of streamflow is a random variable that depends on *P* and *Q*. A non-parametric estimator of precipitation elasticity was defined by [35]. Another similar study conducted by [62] defined an estimator of precipitation elasticity of streamflow by Equation (10) given below:

$$\frac{\Delta Q\_i}{\overline{Q}} = \varepsilon\_p. \frac{\Delta P\_i}{\overline{P}} \tag{10}$$

where Δ*Qi* and Δ*Pi* represents a change in annual streamflow and precipitation in comparison to the long-term average of *Q* and *P*, respectively. The current study used historical data of precipitation and streamflow from 1971 to 1996 for the calculating the precipitation elasticity. It also considers the impacts of evapotranspiration in the climate elasticity model by adding evapotranspiration to Equation (10), resulting in Equation (11):

$$\frac{\Delta \mathbf{Q}\_i}{\overline{\mathbf{Q}}} = \varepsilon\_p. \,\,\frac{\Delta P\_i}{\overline{P}} + \varepsilon\_{\varepsilon t}. \,\,\frac{\Delta ET}{\overline{ET}}\tag{11}$$

where *ε<sup>p</sup>* and *εet* designate precipitation and evapotranspiration elasticity of streamflow, respectively, in multiple regression systems whereas Δ*ET* represents a change in annual mean evapotranspiration in comparison to the long-term mean evapotranspiration (*ET*).

**Figure 4.** Tomer Schilling framework.

#### *2.9. Statistical Model*

The linear regression between averaged annual precipitation (*Pref*) and annual runoff (*Qref*) can be represented by Equation (12):

$$Q\_{ref} = \mathbf{a}P\_{ref} + \mathbf{b}ET + \mathbf{c} \tag{12}$$

where 'a and b' represents the regression equation constant and 'c' denotes the regression intercept. The coefficients of the equation are determined by the least square method. Climate variability is due to external factors such as a, change in LULC and climate variability as represented by Equation (13):

$$
\Delta Q\_{total} = \Delta Q\_p + \Delta Q\_L \tag{13}
$$

where Δ*Qtotal* denotes a change in observed mean annual streamflow which is the result of a change in streamflow due to climate change (Δ*Qp*) and anthropogenic activities (Δ*QL*). Equation (14) specifies the difference between streamflow of the reference period and change period and the output is a change in average annual streamflow Δ*Q*.

$$
\Delta Q = Q\_2 - Q\_1 \tag{14}
$$

where Δ*Q* denotes the change in average annual streamflow, *Q*<sup>1</sup> denotes the streamflow of the reference period whereas *Q*<sup>2</sup> denotes the streamflow during the change period.

#### *2.10. Hydrological Modeling*

The Soil and Water Assessment Tool (SWAT) is a semi-distributed river basin scale model that can be used to compute the impacts of land-use changes by using climatic and streamflow data over a long period [32,63]. The SWAT model is based on the physical characteristics of the basin. In the Arc SWAT model, a watershed is divided into multiple subbasins and then each subbasin is further sub-divided into Hydrological Response Units (HRUs). The SWAT model simulates streamflow based on the following water balance Equation (15) [64]:

$$SW\_t = SW\_0 + \sum\_{i=1}^{t} (R\_{day} - Q\_{surf} - E\_a - W\_{sep} - Q\_{\mathcal{g}w}) \tag{15}$$

where *SWt* denotes the final soil water content (mm); *SW*<sup>0</sup> denotes the initial soil water content on the day *i*; '*t*' designates the time (days); *Rday* specifies the amount of precipitation on the day *i* (mm); *Qsurf* represents the amount of precipitation on the day *i* (mm); *Ea* denotes the amount of evapotranspiration on the day *i* (mm); *Wseep* specifies the amount of water entering the vadose zone from the soil profile on the day *i* (mm) and *Qgw* represents the amount of return flow on the day *i* (mm).

The simulated runoff in the SWAT model can be achieved at the HRUs level by using the Soil Conservation Service (SCS) curve number method [65]. This method mainly depends on the soil, land cover, and hydrological characteristics. The basic equation of the SCS curve number method:

$$Q\_{surf} = \frac{\left(R\_{day} - 0.2S\right)^2}{\left(R\_{day} + 0.8S\right)}\tag{16}$$

where *Qsurf* denotes the daily surface runoff (mm), *Rday* specifies the daily rainfall (mm) and *S* represents the retention parameter (mm). The parameter '*S*' changes due to change in soil, land-use, water conservation measures, slope, and also with time due to variations in the soil water contents of the watershed. The retention parameter can be defined by Equation (17):

$$S = 25.4 \left( \frac{1000}{\text{CN}} - 10 \right) \tag{17}$$

where *S* denotes the retention parameter which gives the value of the drainable volume of soil water per unit area of the saturated thickness (mm/day); *CN* represents the curve number. A more detailed explanation of the SWAT model can be found [65].

The current study divided the KRB into 15 subbasins and 223 HRUs were created based on the areas with homogeneous soil types, vegetation, and slope (Figure 5). The number of subbasins divisions depends on the stream networks of the watershed. Generally, the runoff values are not affected by the number of subbasins [66]. The land use, soil, and slope classification of the KRB in the Arc SWAT model are displayed in Figure 5. The major land-use type in the KRB is pasture whereas the major soil type is lithosols. The slope map shows that a major part of the basin lies at a slope greater than 60% (Figure 5).

#### 2.10.1. Model Setup

The SWAT model requires the preparation of the forcing data such as soil type, land use, DEM, and climatic data. The key steps in the application of the SWAT model are displayed in Figure 6 and include; (a) watershed delineation and subbasin feature derivation; (b) HRU definition; (c) climatic inputs and weather generator; (d) simulation of the SWAT model; (e) sensitivity and uncertainty analysis, and (f) calibration and validation of the SWAT model.

The surface runoff volume was computed using the SCS curve number method [67]. The Penman-Monteith method was used for the calculation of evapotranspiration [68,69]. The SWAT model was run at the monthly time scale. The SWAT calibration and uncertainty program (SWAT-CUP) [70] was used for calibration, validation, sensitivity, and uncertainty analysis using the sequential uncertainty fitting algorithm (SUFI-2 algorithm). The SWAT model was calibrated for the period (1972–1981) and validated (1983–1996) against the observed flows at Gari Habibullah streamflow gauging station in SWAT-CUP.

The sensitivity analysis was performed in SWAT-CUP to determine the parameters that streamflow is most sensitive to. A total of 18 parameters of the SWAT model that relate to surface runoff, groundwater, and snowmelt were selected for model calibration in Table 3. The parameters related to snowmelt such as SMFMX, SMFMN, SFTMP and SMTMP, TLAPS, and others such as CN2, PLAPS, and ALPHA\_BF were found to be more sensitive in the KRB as the basin is mainly snow-fed.

#### 2.10.2. Model Performance and Evaluation Criteria

The performance of the SWAT model was assessed using the Nash-Sutcliffe model efficiency (NSE) and coefficient of determination (*R*2):

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} \left(Q\_0 - Q\_s\right)^2}{\sum\_{i=1}^{n} \left(Q\_0 - \overline{Q}\_o\right)^2} \tag{18}$$

$$\mathcal{R}^2 = \left[ \frac{\sum\_{i=1}^n \left( Q\_o - \overline{Q}\_o \right) \left( Q\_s - \overline{Q}\_s \right)}{\sqrt{\sum\_{i=1}^n \left( Q\_o - \overline{Q}\_o \right)^2 \sum\_{i=1}^n \left( Q\_s - \overline{Q}\_s \right)^2}} \right]^2 \tag{19}$$

where *i* denotes the time step and *n* specifies the total number of simulated time steps. *Qo* and *Qs* represent the observed and simulated streamflow values, respectively. NSE specifies how well the graph between observed and simulated data fits the 1:1 line. The

value of NSE is considered satisfactory when NSE is higher than 0.5 [71]. The greater value of NSE, the better is the model accuracy. The model's simulations are considered perfectly fit if NSE = 1 [72].

Whereas, *R*<sup>2</sup> measures the strength of the linear relationship between the simulated and observed streamflow time series. *R*<sup>2</sup> value can range between 0 and 1. The model is considered ideal if the *R*<sup>2</sup> value is equal to 1, and satisfactory if *R*<sup>2</sup> > 0.6 [71].

**Figure 6.** Flowsheet diagram of the SWAT model.

**Table 3.** List of parameters used for calibrating of the SWAT model.



**Table 3.** *Cont*.

Note: where 'v' designates that the parameter value is replaced by a given value, whereas 'r' specifies that the parameter value is multiplied by (1 + a given value).

#### **3. Results and Discussion**

#### *3.1. Trends in Hydro-Climatic Variables*

The statistical indices explaining the annual variations in observed precipitation, evapotranspiration, and runoff (from 1971 to 2010) of the KRB are listed in Table 4. Table 4 shows the combined annual trend analysis of all the variables. It was observed that the coefficient of variance (Cv) of precipitation was higher than that of runoff and evapotranspiration whereas the Cv of evapotranspiration was found the lowest i.e., 0.12 as compared to precipitation (0.22) and runoff (0.21). Moreover, the ratio of precipitation to runoff was found smaller than the ratio of precipitation to evapotranspiration. Overall, the mean value of precipitation was found to be 14% higher than runoff whereas runoff was found to be 70% higher than evapotranspiration. The skewness of evapotranspiration was observed to be 74% higher than runoff while 72% greater than precipitation. The skewness of precipitation was found as 10% higher than runoff (Table 4).

**Table 4.** Statistical indices for annual variations in precipitation, evapotranspiration, and runoff of the KRB.


The results of the MK test on the precipitation, evapotranspiration, and runoff time series are presented in Table 4, and are also illustrated in Supplementary Figure S1 (for the entire data record), and Figure S2 (for the pre- and post-change time-periods separately). Precipitation and evapotranspiration exhibited statically significantly increasing trends, whereas no statistically significant trends were observed in the runoff time series.

The results of the ITA supports the trend analyses conducted using the MK test, particularly for precipitation, where all data points are above the 1:1 line (Table 5, Figure 7).

**Table 5.** Results of the MK trend test and ITA on the annual for the hydro-climatic time series of the KRB.


Note: 'Sig.' stands for statistically significant and'+.' Refers to the 0.1 level of significance.

**Figure 7.** ITA of precipitation, evapotranspiration, and runoff time series of the KRB. (**a**) Precipitation; (**b**) Evapotranspiration; (**c**) Runoff.

The results of the recent trend analyses by Latif et al. [73] which was conducted over the Jhelum River Basin (JRB) also located in Northern Pakistan and extending into India, were found to be consistent with current study as they detected increasing trends in annual precipitation and evapotranspiration. Similarly, [74] examined the flow regime of the JRB and found trends similar to those in our study in the runoff of a few stations, whereas trends contrasting those of our study were identified at a few weather stations. Furthermore, [75] found an increase in the frequency of wet precipitation days from the northeast to the middle parts of the northern highlands of Pakistan. The trends in temperature and precipitation of our previous study i.e., [45] were found consistent with our current study but runoff in the upper part of the KRB was found to increase in comparison to the lower part (downstream) of the basin. Moreover, [3] also identified the lack of a trend in the runoff of the upper and middle stream of the Syr Darya River basin, and their results were found consistent with our current study.

#### *3.2. Abrupt Changes in the Hydrological Time Series*

The abrupt change or change point in the annual precipitation, evapotranspiration, and runoff of the KRB is displayed in Figure 8. The accumulative difference curve and double mass curve methods were used to identify the abrupt change in the runoff. An abrupt change in all three-time series was identified in 1996 (Figure 8). A study conducted by [34] used the accumulative anomaly method to determine the abrupt change in the runoff of the Huangfuchuan River, China, and their results were found to be consistent. Moreover, a study performed by [76] also observed an abrupt change in the runoff of the Wei He River basin in the year (i.e., 1993) as similar to our study. A study conducted by [74] also determined the change point in streamflow which was consistent with our results.

Table 6 shows the smaller Cv of during the post-change period in comparison to the pre-change period. The Cv of evapotranspiration, for its part, is higher during the post change than in the pre-change period, while for runoff, it is also higher during the later period.

**Figure 8.** Change point analysis on the annual precipitation (ADP), evapotranspiration (ADET), and runoff (ADR) time series.



Regarding the pre-change period, the precipitation (−1%) and runoff (−2.59%) decreased whereas evapotranspiration (+3.25%) increased during the post-change period (Table 6, Figure 9). The ratio of Cv for precipitation to runoff was observed to be 20% higher in the post-change period whereas the ratio of Cv for evapotranspiration to runoff was 18% less during the post-change period in relation to the pre-change period.

The abrupt change was also determined by plotting the data of cumulative precipitation and cumulative runoff in double mass curve analysis (Figure 10). Figure 10a shows the deviation of the double mass curve from the linear relation.

The MK test for the pre (1971–1996) and post-change (1997–2010) periods for precipitation, evapotranspiration, and runoff are displayed in Figure S3 (Supplementary Material) and Table 7 From the results of the MK trend test and ITA method, one can see increasing trends in precipitation and runoff but a decreasing trend, and lacking statistical significance, in evapotranspiration during the pre-change period. After the change point, both the MK test and ITA show statistically significant increasing trend for runoff and an insignificant increasing trend for precipitation and no trend for evapotranspiration (Table 7). The status of evapotranspiration changed from decreasing trend (pre-change) to a trendless (post-change) period which showed that evapotranspiration has increased during the post-change period in comparison with the pre-change period.

Liang [77] investigated the trends in precipitation, temperature, and runoff time series upstream of the Minjiang River Basin, China and observed an increasing trend in precipitation for a similar period to that of our study. They also observed an increase in runoff during the post-change period. Another study conducted by Latif [74] also identified increasing trends in several hydrological stations of the Indus, Jhelum, and Kabul River basins. They observed decreasing trends in precipitation in recent decades as compared to the reference period, which is consistent with the trends of precipitation for the post-change period of our study. Moreover, Khattak [78] observed increasing trends in temperature and runoff but noted an inconsistent pattern for precipitation trend in UIB, Pakistan. Their precipitation and temperature innovative trends were found consistent, whereas runoff innovative trends were found conflicting with the findings of [45].

**Figure 9.** Temporal variation in precipitation, evapotranspiration, and runoff over the KRB.

**Figure 10.** The determination of abrupt change in a time series using (**a**) a double mass curve and; (**b**) a flow duration curve analysis.


**Table 7.** Results of the MK trend test and ITA on the hydro-climatic variables of the KRB during the pre and post-change periods.

Note: \*\* refers to the 0.01 significance level and \* to the 0.05 significance level.

#### *3.3. Runoff Response to Climate Change and Anthropogenic Activities using Different Methods*

In the eco-hydrological framework, the variations in runoff are assessed by considering similar precipitation and evapotranspiration of the basin in different pairs of years. This framework was used to quantify the contribution of climate change and anthropogenic activities in response to variations in the runoff. In this analysis, the abrupt change was also observed in the year 1996 (Figure 11). Moreover, in this framework, Pex was found to decrease from 0.88 to 0.84 while Eex was found to decrease from 0.45 to 0.37 for pre and post change period as shown in Figure 11. The decrease in both excess energy and excess water implies that anthropogenic activities are pronounced in this region.

**Figure 11.** Eco-hydrological approach applied to the KRB.

In the modeling approach, the SWAT model was successfully calibrated (1972–1981) and validated (1983–1996) on the KRB on a monthly time scale as shown in Figures 12 and 13. Overall, the SWAT model simulated both low and high flows very well during both calibration and validation. Table 8 shows the results of the calibration and validation of the SWAT model, which confirm that the SWAT model can be used with confidence in simulating the streamflow of the KRB.

**Figure 12.** Model performance during the calibration stage.

**Figure 13.** Model performance during the validation stage.

**Table 8.** Statistical evaluation of the calibration and validation of the SWAT model on the KRB at a monthly time scale.


Moreover, the calibrated SWAT model was used to reconstruct natural flows for the post-change period without changing the calibrated parameters, however, only temperature and precipitation data were changed according to the period (Figure 14). Figure 14 depicts that the SWAT model simulated the streamflow of the post-change period as satisfactorily. The method of reconstructing natural runoff was applied to the SWAT model to find out the impacts of climate change and anthropogenic activities on the streamflow (Figure 14). The results obtained from the SWAT model suggest that the impacts of anthropogenic activities on streamflow variations of the KRB are evident as compared to climate change (Table 9).

**Figure 14.** The graph between observed and reconstructed flow using SWAT model for the postchange period.

**Table 9.** Differentiating the impacts of climate change and anthropogenic activities on the runoff of the KRB using different methods.


The statistical model was developed for the baseline period (1971–1996) following Equation (20).

$$R\_{Ref} = 0.69P + 0.81ET - 60.89 \tag{20}$$

The statistical model also suggests that the anthropogenic activities are mainly responsible for streamflow variations in the KRB (Table 9).

Similarly, the climate elasticity model was calibrated for the reference period (1971–1996) as shown in Equation (21). For the development of the climate elasticity model, the model parameters were determined from the natural period. The climate elasticity model for the KRB is given by:

$$\frac{\delta \mathcal{R}\_i}{\overline{\mathcal{R}}} = 0.53 \, \frac{\delta P\_i}{\overline{\mathcal{P}}} + 0.51 \, \frac{\delta ET}{\overline{ET}} \tag{21}$$

From the climate elasticity model one can see that anthropogenic activities played a major role in streamflow variations of the KRB during the post-change period as compared to climate change. A study conducted by [79] also investigated that the streamflow variations in the Soan River basin, Pakistan is mainly attributed to anthropogenic activities. The results of their study were found similar to our findings. Another study performed by [80]

determined that the relative contribution of climate change and land-use change in streamflow variations in the Tarbela catchment was 39.3% and 60.7%, respectively. Similarly, [81] determined that the impacts of human activities on runoff variability of the Dongjiang River basin were obvious and their findings were found consistent with our findings of the current study. The study of [76] also identified human activities as a major factor for runoff variations of the upper Wei He River basin, China by using similar methods as we did. The results of several other studies such as [32,33,39,77] are also in agreement with the results of our study.

#### *3.4. Discussion*

The correlation between runoff coefficient and precipitation, evapotranspiration, and drought index in the KRB is shown in Figure 15. There could be a significant or nonsignificant correlation between the annual runoff coefficient and precipitation. The significant correlation depends on soil type and vegetation characteristics of the watershed area. Runoff also depends on infiltration which is indirectly influenced by the weather condition; geography, soil texture, and land cover [11]. In the current study, the mean annual data were used for correlation analysis. Figure 15a displayed that more than 50% points of the runoff coefficient vs. precipitation occurred in the second quadrant of the coordinate system, whereas few points were found in the first, third, and fourth quadrant, and this condition suggested the weak relationship between runoff coefficients versus precipitation. However, the scatter plot between runoff coefficient and evapotranspiration displayed a strong relationship as compared to runoff coefficient versus precipitation (Figure 15b). Moreover, the correlation between precipitation and runoff was found to be 0.90 whereas it was observed 0.1 for evapotranspiration and runoff for the study period. However, the precipitation was found to be negatively correlated with runoff coefficient and drought index whereas the evapotranspiration was found to be positively correlated with drought index (0.21) and negatively correlated with runoff coefficient (0.37). Furthermore, strong correlation (0.37) was observed between evapotranspiration and drought index during the post-change period, whereas a weak correlation was found during the pre-change period. A weak correlation was found between precipitation and drought index during the post-change period. Moreover, a negative correlation was observed between evapotranspiration and runoff coefficients during both pre-change (−0.46) and post-change (−0.26) periods (Figure 15) whereas precipitation was found to have a positive impact on runoff coefficient during the post-change period. A study performed by [11] determined the correlation of annual time series of streamflow, precipitation, and temperature. The coherency analysis of KRB was also found to be consistent during the pre-change period. However, the coherency analysis displayed a different behavior for the post-change period as compared to the pre-change period. The coherency analysis did not consider the other factors except precipitation and evapotranspiration. Moreover, the results of coherency of streamflow of Colombian Pacific Basins [10] were found consistent with the results of our study. It was also unveiled that the influence of climate was dominant on streamflow in this region.

Past studies [15,82,83] conducted in different parts of the world also explained that runoff is more sensitive to precipitation as compared to evapotranspiration. However, [84] identified that the coefficient for soil water content and maximum soil water holding capacity were key factors influencing the long-term hydrological response. Moreover, it was observed that under different climatic conditions, the impact of precipitation and evapotranspiration on hydrological response was varied [85], which further indicates the different mechanisms for runoff generation. However, the correlation analysis of runoff does not consider the geography, vegetation characteristics, soil type, and other characteristics of KRB. The correlation analysis only considers the relation between hydrometrological variables [10,11].

**Figure 15.** Scatter plots of (**a**) runoff coefficient vs. precipitation; (**b**) runoff coefficient vs. evapotranspiration; and (**c**) runoff coefficient vs. drought index.

Large variations have been observed in the land-cover types of the KRB since 1992 as shown in Figure 16. A large part of the basin had been transformed from natural forest/grassland into agricultural lands and residential areas(Figure 10). A study conducted by Ramírez [86] observed that forests usually evaporate more water than other types of plants (such as agricultural and annual crops), and a decrease in forest area is subjected to increased runoff of the basin. Naran and Kaghan are tourists' spots in the KRB and these areas produce environmental impacts associated with travel, accommodation, and recreational activities. To minimize the environmental degradation associated with tourism and recreation it may require appropriate land-use zoning, regulation and surveillance of access and activities, direct physical protection of particular areas and education both on-site and elsewhere. In addition, it is important to provide incentives to encourage low-impact types of recreation, such as contemplative, naturalist, and wilderness travel activities, and discourage high-impact type recreation such as sports, social activities, motorcycles heavy vehicles, and accommodation involving building and engineering construction. Unmanaged recreational activities and increased encroachment in forest land were found to have an impact on runoff characteristics of the KRB. The population has also increased from 0.77 million (1981) to 1.15 million (1998) in the KRB. The current population of the basin has reached 1.6 million.

**Figure 16.** Land cover variations in KRB during the period 1992–2014.

Moreover, several storage dams and small ponds have been built in the KRB in recent years. The different hydropower projects have been initiated under the United Nations (UN) goal of sustainable clean energy and the Paris agreement. Several hydropower projects have been started by the government of Khyber Pakhtunkhawa such as the Balakot project (190 MW), Kari-Muskhui (446 MW), Naran Dam (210 MW), Batakundi (65 MW), Laspur-Muri Gram (130 MW), Shushghai-Zhendoli (144 MW), Shogo Sin (132 MW), Torkum-Gudubar (409 MW) and Samshel Toren (260 MW) [87].

#### *3.5. Comparison with Other River Basins*

In the last fifty years, a number of studies have related abrupt changes in river runoff to climate change in different parts of the world [75,76] and anthropogenic activities [13]. Several researchers found that the watershed characteristics were the main driving factors that influenced the hydrological response [82]. Another study conducted by [88] determined the runoff characteristics in different environmental conditions of the United States. Moreover, [89] observed the impacts of climate change on the hydro-climatology of Lake Tana Basin, Ethiopia. Another study performed by [90] observed that climate change and anthropogenic activities are mainly responsible for influencing the runoff characteristics of the Tualatin River basin, Oregon. Similarly, [91] observed the remarkable changes in the hydrological response of the Swedish Rivers due to the impacts of anthropogenic activities and climate change. The grassland plains in Russia showed a greater delay in runoff due to watershed characteristics observed by [92].

Moreover, few Asian rivers are largely attributed to watershed characteristics, such as the Yellow River [39,93], Ganga River [94], Huifa River [32], Kofarnihon River [95], Bagmati River [2], and Pearl River [96], and the results of these studies were found consistent with our study. Another study conducted by [79] observed that just after the change point, reduction in Soan River runoff was due to climate change and land cover changes. Similarly, a study conducted by [80] on the Tarbela catchment found that watershed characteristics play a major role in changes in runoff and climate change relatively contributes less as compared to watershed characteristics (anthropogenic activities). Moreover, [97] observed that watershed characteristics were evident for the increased runoff in the Simly watershed. Under the different climatic zone, there is a need to plan watershed-scale water resources management to face global climate change and frequent extreme precipitation as well as land use degradation. The national self-regulation ability is not strong enough to obtain the key object of hydrologic protection and management for climate-sensitive areas.

This study examined the impact of climate change and anthropogenic activities on the hydrology of a river basin. It is necessary to determine the impact of individual anthropogenic activities on the hydrological response. This research also assumed that climate variability and anthropogenic activities are independent, which is not the case.

#### **4. Conclusions**

This study investigated the impacts of climate change and anthropogenic activities on the runoff characteristics of the KRB during the period 1971–2010 using empirical, statistical, and modeling techniques. Potential causes for the observed changes in streamflow were identified and, the following conclusions were drawn:


post-change period, only an increasing trend in runoff was found to be statistically significant; no trend was seen in the evapotranspiration time series, while the increasing trend in precipitation was not found to be statistically significant.


**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/w13223163/s1, Figure S1: MK trend analysis of Runoff (R), Precipitation (Pr), and Evapotranspiration (Et) for whole period (1971–2010); Figure S2: MK trend analysis of Precipitation (Pr), Evapotranspiration (Et), and Runoff (R) during pre-change (1971–1996) and post-change (1997–2010) periods; Figure S3: ITA analysis of precipitation (a,b), evapotranspiration (c,d), and runoff (e,f) during pre-change (1971–1996), and post-change (1997–2010) periods.

**Author Contributions:** Conceptualization and data analysis, M.S.; Hydrological Modeling and abstract, M.S. & M.A.; Introduction, M.S., A.S.G., M.Z. & M.I.K.; S.M., Methodology, M.S. & A.W.; Idea, supervision and funding, S.L. & M.S; Results, M.S. & M.A.; Discussion, M.A., M.S. &; S.L. Paper formatting and conclusion M.S., M.A. & A.S.G. All the authors contributed to the finalization of this manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is supported by the NSFC-ICIMOD joint project (Grant no. 41761144075), Yunnan University grant (2018M643541 and 209071) and Higher Education Commission (HEC) of Pakistan (HEC/STR/279).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study is available on request from the corresponding author.

**Acknowledgments:** We are thankful to the Surface Water Hydrology Project of the Water and Power Development Authority (SWHP-WAPDA) and the Pakistan Meteorological Department for providing the required hydro-meteorological data to conduct this study. We are thankful to editor and associate editor, as well as reviewers for their valuable comments and suggestions, which have improved the final version of this manuscript.

**Conflicts of Interest:** The authors declare that there is no conflict of interest.

#### **Abbreviations**


#### **References**

