This study investigates changes in the hydrologic regime of the Chitral River, Hindukush–Karakoram–Himalayan (HKH) region, Pakistan. Different statistically based methods were used to assess climate change-induced hydrologic alterations that can possibly impact aquatic habitat in the study region. The hydrological model Hydrologiska Byråns Vattenbalansavdelning (HBV) was calibrated, validated, and applied to predict streamflow in the Chitral River basin. The HBV model was forced with the ensemble of four general circulation models under different representative concentration pathway emission scenarios to generate future streamflow under climate change conditions in the basin for the mid-twenty-first century. The results of this study show that hydrologic regimes in the study area, expressed by the magnitude, duration, frequency, timing, and rate of streamflow, are likely to alter in the future. Positive (i.e., with increased frequency) hydrologic alteration is projected for most flow parameters under all scenarios for the 2021–2050 period compared with values observed during the historical period (1976–2005). These hydrologic alterations might have impacts on fish and migratory bird species in the study area. This research can be helpful in providing practical information for more effective water resources and aquatic ecosystem management in the HKH region.
Investigation of changes in the hydrologic regime of the Chitral River.
Employment of IHA and RVA methods to evaluate riverine ecosystem health.
Hydrological modeling using HBV-light to generate future streamflow.
Assessment of climate change impacts on different streamflow characteristics.
Analysis of streamflow changes and their relationship with aquatic habitat.
Climate has substantially changed compared to pre-industrial times and is projected to further alter considerably worldwide in the twenty-first century at the global scale (Hoegh-Guldberg et al., 2018). Changes have also been observed in Pakistan, and extreme climatic events are predicted to increase especially in the Hindukush–Karakoram–Himalayan (HKH) region of Pakistan. The HKH region has the most overstressed aquifers in the world and is a climate change hotspot; therefore, the predicted extreme events will pose additional threats to the region, already exposed to numerous and severe flood events (Lutz et al. 2016). Increase in precipitation and average summer temperature is projected for catchments of the Upper Indus Basin (Archer & Fowler 2004; Fowler & Archer 2006; Bocchiola & Diolaiuti 2013; Ahmad & Hussain 2017), often associated with alterations of the hydrologic cycle (Bell et al. 2016). Climate change impacts are stronger in magnitude and outweigh all other direct impacts (such as dams and water withdrawals) on the water cycle (Finke et al. 1998; Döll & Zhang 2010). Several studies have been conducted at various spatio-temporal scales on some major rivers in the HKH region revealing significant impacts of climate change on catchment streamflow (Wijngaard et al. 2017; Hashmi et al. 2019; Li et al. 2019; Shrestha & Nepal 2019).
The natural ﬂow regime of rivers is closely associated with the biodiversity and integrity of the riverine ecosystems (Poff et al. 1997; Petts et al. 2006; Roy et al. 2006) and it controls many physio-ecological processes (Arthington et al. 2006; Naiman et al. 2008; Guastini et al. 2019). These processes include the exchange of nutrients and the transport of sediments, which in turn inﬂuence factors such as velocity and depth of ﬂow, which are critical elements of river habitats (Poff et al. 1997; Benda et al. 2004; Shiau & Wu 2007). Therefore, the flow regime of a river is important for sustaining the surrounding environment and for the functionality of aquatic, riparian, and wetland ecosystems (Poff et al. 1997; Petts et al. 2006; Roy et al. 2006). Alterations in ﬂow regime can have negative impacts on stream ecosystems, including failed fish and micro-invertebrate recruitment, extinction of local species, and the successful invasion of exotic species (Bunn & Arthington 2002; Poff et al. 2007, 2010; Poff & Zimmerman 2010).
Fish are one of the most important constituents of biological diversity in the HKH region (Ahmed & Joyia 2003). Approximately 193 fish species have been identified in the Pakistani rivers, of which 31 species have economic importance (Rafique & Khan 2012). Eight species are widely distributed in the HKH region, including Oncorhynchusmykiss, Salmo trutta fario, Schizothorax plagiostomus, Diptychus maculates, Ptychobarbus conirostris, Racoma labiata, and Schizopyge esocinus (Rafique & Khan 2012). Flora represents another important source of biodiversity in the HKH region and plays a vital role in the catchment water cycle in terms of soil conservation and evapotranspiration fluxes. Both fish and plant biodiversity in the Chitral River are seriously exposed to ecological threats due to hydrological alteration. For instance, the Chitral River shows poor fish diversity associated with changes in water temperature and high turbidity (Rafique 2001; Bari et al. 2014). Moreover, deterioration of water quality in the Chitral River has been observed and is associated with physical factors depending on climate change (Nafees et al. 2012).
Despite the important ecological and societal implications of the climate change impacts in the Chitral River basin, previous studies so far have focused, as far as we know, only on hydrologic alterations related to changes in mean and extreme streamflow values but more detailed analysis of hydrological modifications induced by climate change and their effects on the riverine ecosystem in this region is missing and is highly needed for operational purposes. Therefore, overall, this study aims to identify and quantify the main sources of hydrologic alteration that might cause pressing effects on aquatic organisms in the Chitral River basin.
MATERIALS AND METHODS
Study region and data availability
The Chitral River has a drainage area of 11,396 km2 and is a transboundary river basin between Pakistan and Afghanistan (Figure 1). Large parts of the basin are at high elevations and host wide glaciers. More than 50% of the catchment area is between 4,000 and 5,000 m a.s.l. (Hayat et al. 2019). Glacier and snowmelt waters (Burhan et al. 2020a) make significant flow contributions to the Chitral River and serve as a critical resource for irrigation and hydropower generation. In addition, the river water is crucial for supporting domestic activities, playing therefore a vital role in the socio-economic development of the region.
Hydrometeorological data in the Chitral River basin were available for the period 1994–2012 (Table 1). A 90-m resolution digital elevation model was also available for the study area. Potential evapotranspiration was computed using the approach by Irmak et al. (2003).
The results from different General Circulation Models (GCMs, Table 2) have been averaged in order to have more reliable precipitation and temperature data inputs in the study region and to carry out more solid interpretation of climate evolution and land surface-atmosphere interactions in this large-scale basin (He et al. 2016). As the outputs of GCMs are typically too coarse to properly represent local sub grid-scale features and dynamics (Coulibaly et al. 2005), they were downscaled using the statistical downscaling model suggested by Su et al. (2016). As GCMs do not always accurately simulate the climate variables, especially in large regions with complex topography as in our study area, biases might be present and should be corrected (Su et al. 2016). The bias of the GCM outputs was corrected by employing the cumulative distribution functions matching method introduced by Li et al. (2010) that was successfully applied in other mountain regions (Su et al. 2016). Data of the GCMs were forced with four respective representative concentration pathways (RCPs), a set of scenarios that include time series of emissions and concentrations of greenhouse gases and aerosols and chemically active gases, as well as land use/land cover (Van Vuuren et al. 2011). The RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 emission scenarios were adopted: the RCP 2.6 represents a low emission scenario, with a radiative forcing of 2.6 w/m2 towards the end of the twenty-first century, the RCP 4.5 is an intermediate (average) emission scenario with a radiative forcing of 4.5 w/m2 towards the end of the century, the RCP 6.0 represents a high emission scenario with radiative forcing of 6.0 w/m2 towards the end of the century, and the RCP 8.5 corresponds to the highest (extreme) emission scenario, with a radiative forcing of 8.5 w/m2 towards the end of the century (Van Vuuren et al. 2011).
The HBV-light model was used to simulate daily streamflow in the Chitral River basin. HBV-light is a modified version of HBV (Seibert & Vis 2012), a lumped, conceptual, and semi-distributed hydrological model developed by the Swedish Meteorological and Hydrological Institute. The modified version is based on the degree-day method to simulate snowmelt based on a linear function of the temperature difference between average air temperature and the threshold temperature for snowmelt and is well suited to describe the runoff variability (Seibert 2001; Radchenko et al. 2013; Bhattarai et al. 2018; Usman et al. 2019; Burhan et al. 2020b). The model has different routines: snow routine, soil routine, response routine, and routing routine. The time series of temperature, precipitation, and potential evaporation were used as inputs to the model to simulate streamflow, as follows: in the snow routine, precipitation was first utilized as input and then simulated by the model either as snow or rain depending on a threshold temperature. In the soil routine, the simulation of groundwater recharge and actual evaporation depended on the actual water storage. Actual soil evaporation matched potential evaporation when the availability of water was not limiting evaporation. When water availability limited evaporation, a linear reduction was utilized. In the response, routine streamflow was calculated as a function of water storage. Finally, in the routing routine, a triangular weighting function was employed for the purpose of simulating the route of the streamflow at the basin outlet. For a comprehensive description of the model including all the associated characteristics and some applications, the reader is referred to the following studies (Bergström 1976; Seibert 2005; Seibert & Vis 2012; Bhattarai et al. 2018; Usman et al. 2019).
Model development and performance assessment
The HBV model was calibrated by employing an automatic calibration method, i.e., the genetic algorithm and Powell optimization method (Seibert 2000). This method is based on a recombination of parameter sets and the best set of parameters resulting in the highest objective function (e.g. Nash–Sutcliffe Efficiency index (NSE)) was selected. This process was repeated several times to obtain the best-optimized set of parameters and improve the model efficiency. Calibration was performed by dividing the available daily time series (1994–2012) into two periods. The first period (1994–2007) was used for calibration of the model, using the first years as the warm-up period. The calibrated model was then validated against the different observed data from the second period (2008–2012). The model was applied for the generation of future streamflow using the average weather data obtained by the four GCMs forced with the different RCPs.
Analysis of streamflow alteration
Potential streamflow alterations caused by changing climate conditions in the study basin were first assessed using the indicators of hydrologic alteration (IHA) method (Richter et al. 1996). The IHA method includes 33 parameters that determine how flow regime is perturbed by different disturbances (Table S1 in Supplementary Materials). Out of 33 IHA indicators, one indicator (number of zero-flow days) was not considered in this study, as no zero-flow days were found in the monitoring period. The remaining 32 parameters were distributed in ﬁve different groups based on magnitude, timing, frequency, duration, and rate of change of streamflow. The ﬁve groups were deﬁned as follows.
Group 1 includes parameters related to the magnitude of water condition (monthly mean ﬂows); group 2 includes parameters representing the magnitude and duration of annual extreme ﬂows; group 3 includes parameters describing the frequency of occurrence of drought or flood conditions; group 4 includes parameters representing the frequency and duration of high and low pulses (deﬁned as the annual periods when the daily ﬂows are above the 75th percentile and below the 25th percentile of daily streamﬂow); and group 5 includes parameters related to the rate of change in water flow.
The computation of the different indices involves four steps: (1) definition of the streamflow data series for pre-impact (1976–2005) and post-impact (2021–2050) periods; (2) calculation of hydrologic attributes for the entire monitoring period (1976–2050); (3) computation of inter-annual statistics (central tendency and dispersion measure); (4) comparison of pre-impact and post-impact period streamflow data and calculation of the percent change for each of the IHA parameters (3rd column in Table S1). Thirty years were selected as the pre-impact period and an additional 30 years were used as the post-impact period; these values were chosen based on the recommendation that at least 20 years should be used to account for natural climatic variability (Richter et al. 1997; Pfeiffer & Ionita 2017).
The degree to which the RVA target ranges are not attained is accepted as a measure of hydrologic alteration. A positive value indicates that the respective parameter values fell within the target range more frequently in the post-impact period than expected, and a negative value indicates that the respective parameter values fell within the target range less frequently in the post-impact period than expected. A hydrologic alteration is zero when the observed frequency of post-impact annual values that fall within the RVA target range equals the expected frequency.
The HBV-light, IHA software version 7.1 (The Nature Conservancy 2009), Microsoft Excel, and OriginPro have been used for generating the models and the statistical analysis.
RESULTS AND DISCUSSION
Calibration and validation
The time series during the model calibration and validation periods and the metrics of model of performance (Figures 2 and 3 and Table 3) indicate a good agreement between the observed and the simulated daily streamflow values, which means the HBV-light performed well both during calibration and validation. The NSE values greater than 0.9 during the calibration and greater than 0.8 during the validation were achieved. The R2 value was above 0.9 for both calibration and validation, PBIAS was 3.7 and −2.0% for the calibration and validation period, IoAd was more than 0.95 for both calibration and validation, and VE was 0.81 for calibration and 0.77 for validation, respectively. All these observations reveal that the HBV model was efficient in replicating catchment processes in the Chitral River basin with very small prediction uncertainties. However, the PBIAS values indicate that the model slightly underestimated streamflow during the calibration and overestimated streamflow during the validation period (Figure 3).
Measure of hydrological alteration and ecological implications
All 32 ecologically relevant (IHA) parameters adequately represented the streamflow. Changes in streamflow parameters projected in future periods by the model indicate climate change impacts on streamflow, with different extents and implications according to the different IHA groups.
Alteration of monthly streamflow
The values of the 12 parameters of the IHA group 1 (monthly streamflow) for the Chitral River under the RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 emission scenarios are listed in Table S2. Ten out of 12 parameters of group 1 returned high values ( > 67%), with especially high values for June, January, and November under the RCP 2.6 emission scenario. Under the RCP 4.5 and 6.0 emission scenarios, nine and six parameters, respectively, of group 1 returned high values ( > 67%). Under the extreme emission scenario, eight parameters of group 1 returned high values. The maximum values of group 1 under each emission scenario are given in Figure 4.
The IHA group 1 is composed of parameters related to the magnitude of monthly streamflow and representing seasonal variations. Most of the parameters of group 1 projected a significant high value of , revealing an increased frequency of high flows in the post-impact period compared to the pre-impact period under all emissions scenarios. From mid-spring to early autumn, snow and glacier melt are major contributors to streamflow in the Chitral River basin (Hasson 2016). However, in the dry season (October to February), most of the streamflow is composed of groundwater, and liquid precipitation does not contribute significantly to streamflow. This increase in values implies that the frequency of these parameters has increased in the post-impact period compared to that in the pre-impact period. Generally, most of the parameters of group 1 indicate an increased frequency of streamflow alteration in the post-impact period compared to the pre-impact period under all emissions scenarios. However, this increase in frequency is more significant in winter, early spring, and early summer in the post-impact period compared to the pre-impact period under all emissions scenarios, suggesting that more values are likely to fall within the target range in the future during these seasons (Figure 4). Furthermore, the magnitude of this increase in frequency is much stronger under the RCP 8.5 emissions scenario, indicating that extreme emissions might be suitable for attaining the targeted range. According to the improved overall degree of hydrologic alteration, the parameters of group 1 (monthly median flows) are likely to be altered most under the extreme emission scenario and least under the high emission scenario.
The most abundant species of fish in the Chitral River (Schizothorax plagiostomus) spawns twice a year from March to April and from September to October (Qadri et al. 1982; Jan et al. 2014). Therefore, significant hydrologic alteration in monthly streamflow (increase frequency of high streamflow values especially in winter) projected for these months potentially affects spawning of these fish species.
Alteration of annual extreme flows
The values of the 12 parameters of the IHA group 2 (annual extreme flow conditions) and the two parameters of the IHA group 3 (annual extreme flow conditions) for the Chitral River under the RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 emissions scenarios are listed in Table S3. Ten out of 15 parameters of the IHA groups 2 and 3 returned high values (>67%), with the highest values under the RCP 2.6 emission scenario. The number of parameters that returned high values (>67%), under RCP 4.5 emissions scenario was 40% less than that of the low emissions scenario. Under RCP 6.0 and RCP 8.5 emissions scenarios, seven parameters returned high values. The maximum values for parameters of groups 2 and 3 under each emissions scenario are given in Figures 5 and 6.
According to the improved overall degree of hydrologic alteration, parameters of group 2 (magnitude and duration of annual extreme flows) are likely to be altered most under extreme emissions scenario and least under average emissions scenario.
There are two parameters in the IHA group 3 (timing of annual extreme flows), namely, the date of minimum flow and the date of maximum flow. The date of minimum flow is the only parameter that projected on average negative low values (<-27%), but this change is not significant compared to other changes (average positive values >194%). On the contrary, the date of maximum flow projected an enormous increase in its frequencies in the post-impact period compared to the pre-impact period, indicating that it is more likely to lie within the targeted range in the future. The change in flow timings can impair ecosystem functioning (Baron et al. 2002) by effecting fish migration and spawning, as well as egg hatching, rearing, movement onto the floodplains for feeding and reproduction, and migration upstream and downstream (Poff et al. 1997; Gao et al. 2018). If the water level reaches too low or too high values during spawning or reproduction months, it can affect the propagation of fish species and likely of other aquatic organisms. Moreover, most of the birds enter and exit Pakistan through the Chitral River, which makes its ecosystem more vulnerable to changes in flow timing, and the projected significant increase in the date of maximum flow might have an impact on the natural routes of different migratory species. These results depict an increase in the frequency of parameters that are relevant to low flows in the post-impact period compared to the pre-impact period under all emissions scenarios; additionally, the results indicate that the quantity of parameters undergoing of high (positive) category is reduced as radiative forcing in climate increases. However, the magnitude of change in the low flow indicators strengthens under the extreme emissions scenarios.
Alteration of magnitude, frequency, and duration of annual extreme flows and flow changes
The values of the four parameters of group 4 (frequency and duration) and the three parameters of group 5 (frequency and duration) for the Chitral River, under all RCPs, are listed in Table S4. Four of the seven parameters in groups 4 and 5 under the RCP 2.6 emission scenario have very high values (>67%). Four parameters returned high values also for the RCP 4.5 emissions scenario, although the magnitude of change is much stronger here than for RCP 2.6. For RCPs 6.0 and 8.5, five and three parameters returned high values, respectively. The maximum values of groups 4 and 5 under each emissions scenario are given in Figures 7 and 8.
Most of the parameters in the IHA group 4 (frequency and duration) have positive projected values. An of up to 200% in low-pulse counts under the RCP 8.5 emission scenario is likely to considerably alter the ecological environment of the Chitral River floodplain by foisting basic restrictions on the aquatic communities, and it may significantly affect the population and diversity of fish species and other organisms.
Hydrologic alteration of all parameters for the four RCPs
The ranked values of all 32 parameters of the five groups for the RCP 2.6, 4.5, 6.0, and 8.5 emissions scenarios are illustrated in Table S5. The maximum was projected for the monthly streamflow indicator in June of group 1, and the minimum was projected for the date of the minimum streamflow parameter of the third IHA group for the RCP 2.6 emission scenario. Under the RCP 4.5 emission scenario, the maximum was projected for high-pulse durations and minimum was projected for the date of the minimum streamflow parameter of group 3. High-pulse duration experienced the maximum and parameter date of the minimum streamflow was the least hydrologically altered parameter for the RCP 6.0 emissions scenario. These results indicate that the days in the year when the minimum streamflow was observed will likely be similar in the future to those in the past, under different emissions scenarios.
The fall rate of streamflow (the rate of change of flow, group 5) remained unchanged in the post-impact period under all respective emissions scenarios. The number of reversals flow parameters (transitions between positive and negative streamﬂow, i.e., the instants when streamflow values start to vary from high to low values (fall rate) and from low to high values (rise rate)) projected a high increase in the post-impact period compared to the pre-impact period under the RCP 6.0 emissions scenario, which is an indicator of enhanced frequency in alternating periods (Pfeiffer & Ionita 2017). An increase in the number of reversals under the high emissions scenario suggests that climate change is likely to result in the modification of streamﬂow characteristics, e.g., high flow variability, of the Chitral River. It is interesting to note that nearly 60% of parameters in the IHA groups 4 and 5 projected a significant increase in their frequency in the post-impact period compared to the pre-impact period, and the average to high emissions scenarios projected stronger magnitudes of change.
Improved overall degree of hydrologic alteration
The improved overall degree of hydrologic alteration () for the five IHA groups and the four RCPs is reported in Figure 9, calculated using Equation (5). The highest overall X values were projected for group 4 (magnitude, frequency, and duration of annual extreme flows) for RCPs 6.0 and 4.5, followed by the first group for RCP 8.5. Consistently lower values were observed for group 5 for all RCPs. Groups 2 and 3 showed relatively similar X values.
Fourteen, 18, 9, and 16 parameters out of 32 indicators returned values of more than 100% in the post-impact period for the low (RCP 2.6), average (RCP 4.5), high (RCP 6.0), and extreme (RCP 8.5) emission scenarios, respectively. The maximum average alteration (323%) was returned for the RCP 6.0 emission scenario and the minimum average alteration (127%) was projected for the RCP 2.6 emission scenario in all IHA groups (Table S6).
In this study, an assessment of hydrological alteration was performed in the flow regime of the Chitral River in the HKH region, Pakistan. The widely used methods of the IHA and the non-parametric RVA were adopted to characterize the degree of alteration and the possible future impacts on the riverine (ecohydrological) system under different future climatic conditions. The HBV-light hydrological model employed in this study showed a good streamflow prediction efficiency. The average results from climatic multi-model ensemble projections in this large basin are used to feed the HBV-light and predict future streamflow regimes more reliably.
The IHA/RVA analysis reflected a range of streamflow characteristics. Natural streamflow conditions were projected to change in the mid-twenty-first century. Our results suggest an increase in the alterations of all five IHA groups, representing the magnitude, duration, frequency, timing of daily streamflow, and the rate of change of streamflow. These alterations might have impacts on the most abundant fish species in the Chitral River, especially during their spawning and development stages, as well as on the migration processes of migratory bird species that are common in the study area.
Different previous studies (e.g., Finke et al. 1998; Döll & Zhang 2010; Chen et al. 2019) revealed that the impacts of climate change are stronger in magnitude and outweigh all other direct impacts on the streamflow (e.g., those related to unsustainable use of surface water) and need to be adequately considered in planning and management of water resources (Serrat-Capdevila et al. 2007). However, it might be difficult to explicitly identify the individual impacts of climate change on hydrologic alteration (Malmqvist & Rundle 2002; Pfeiffer & Ionita 2017) and this limitation applies also to this study. Despite this limitation, the work adequately describes the overall streamflow regime alteration due to climate change in the Chitral River basin, offering potentially useful information to different stakeholders responsible for policy and decision making for defining the vulnerability of aquatic ecosystems to hydrologic alteration and for more effective and sustainable water resources and aquatic ecosystem management. This study improves our knowledge and understanding of the relationship among aquatic habitat, streamflow alteration, and climate change, and can provide references for further studies in this region. More research on the quantification of climate change impacts is necessary, possibly adopting more integrated approaches that consider the combined effect of multiple sources of alteration on river streamflow.
Funding for this work was provided by Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA20060603) and the National Natural Science Foundation of China (Grant No. 41471292). Authors are grateful to Strategic Priority Research Program of the Chinese Academy of Sciences and the National Natural Science Foundation of China.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.