Selecting and Downscaling a Set of Climate Models for Projecting Climatic Change for Impact Assessment in the Upper Indus Basin ( UIB )

This study focusses on identifying a set of representative climate model projections for the Upper Indus Basin (UIB). Although a large number of General Circulation Models (GCM) predictor sets are available nowadays in the CMIP5 archive, the issue of their reliability for specific regions must still be confronted. This situation makes it imperative to sort out the most appropriate single or small-ensemble set of GCMs for the assessment of climate change impacts in a region. Here a set of different approaches is adopted and applied for the step-wise shortlisting and selection of appropriate climate models for the UIB under two RCPs: RCP 4.5 and RCP 8.5, based on: (a) range of projected mean changes, (b) range of projected extreme changes, and (c) skill in reproducing the past climate. Furthermore, because of higher uncertainties in climate projection for high mountainous regions like the UIB, a wider range of future GCM climate projections is considered by using all possible extreme future scenarios (wet-warm, wet-cold, dry-warm, dry-cold). Based on this two-fold procedure, a limited number of climate models is pre-selected, from of which the final selection is done by assigning ranks to the weighted score for each of the mentioned selection criteria. The dynamically downscaled climate projections from the Coordinated Regional Downscaling Experiment (CORDEX) available for the top-ranked GCMs are further statistically downscaled (bias-corrected) over the UIB. The downscaled projections up to the year 2100 indicate temperature increases ranging between 2.3 ◦C and 9.0 ◦C and precipitation changes that range from a slight annual increase of 2.2% under the drier scenarios to as high as 15.9% in the wet scenarios. Moreover, for all scenarios, future precipitation will be more extreme, as the probability of wet days will decrease, while, at the same time, precipitation intensities will increase. The spatial distribution of the downscaled predictors across the UIB also shows similar patterns for all scenarios, with a distinct precipitation decrease over the south-eastern parts of the basin, but an increase in the northeastern parts. These two features are particularly intense for the “Dry-Warm” and the “Median” scenarios over the late 21st century.


Introduction
Future climate projections provided by general circulation models (GCMs) can serve as the basic input for climate change impact studies on water resources.As the outputs from these general circulation models (GCMs) have only coarse spatial resolution, and so are often not suitable as direct input to distributed or semi-distributed hydrologic models, they have to be downscaled in most cases to appropriate (higher) resolutions.Such a downscaling can be done either through applying statistical downscaling or through dynamical downscaling via use of a regional climate model (RCM) embedded in a larger GCM.
Despite the availability of a large number of GCM outputs in the CMIP5 archive, and the on-going improvements in their process representations, issues of large uncertainties with regard to the future climate are not yet avoidable.The inherent uncertainties, along with other factors such as time limitations, human resource availability, or computational constraints, make it imperative to sort out the most appropriate individual GCM or small ensemble of GCMs suitable for downscaling and subsequent use in the assessment of climate change impacts.
This aforementioned selection of GCMs is not simple or straightforward, as there can be nearly an unlimited number of criteria and approaches through which climate models can be evaluated for their skill and suitability for specific purposes and regions.In most cases, though, the selection can be based either on a single criterion or a whole set of criteria.One approach may be to consider the total change projected by the GCMs, in the means and/or extremes of a climate variable and its location on the overall spectrum of the future projected by all GCMs.Another approach may place more emphasis on the success of GCMs in simulating past climate for either the means, extremes, or seasonality [1,2] of the study region.Additionally, there may be approaches based on some combination of the aforementioned approaches.The first approach, which considers all the possible projected futures (stretching from warm and wet to cold and dry, or opting for the middle path of all possible futures) is becoming more relevant, especially in regions such as the Hindu Kush Himalayas (HKH) and UIB, where GCMs/RCMs have been reported to struggle in simulating the past climate [3][4][5][6].As no individual model can be separated out as superior in simulating the past climate in the HKH region, it is therefore important to consider the full range of possible projected futures when focusing on assessments of climate change impacts.
The criteria to be used for selecting the most appropriate model runs are also defined based on their intended purpose or the region.Both of these factors are important, as a different intended uses my require consideration of assessment based on totally different skills or variables, while the importance of a specific selection criteria may differ for different locations and topographically contrasting areas.Additionally, as not all the available models may be equally good for specific locations, regions or topographies, the need for the assessment of the ability of climate models to reproduce important processes in the study region is vital and essential.
In the current study, we consider a combination of these approaches to shortlist climate model runs, along with utilizing new and improved data for the past climate in the UIB [7] for assessment of model skill in simulating the seasonal cycles in the region.The main aim of the study was to select a set of GCM simulations that can represent the full spectrum of the future climate, as projected by the entire pool of climate models, in term of both means and extremes, and which can be subsequently used as climate forcing for hydrological modelling to assess a wider range of possible climate change hydrological impacts, especially for the expected changes in water yield, annual cycle, high and low flows, and floods.
The specific objectives of the study included: 1.
To devise a procedure for the identification/filtering of a limited number of climate model runs that can represent the full spectrum of future climate as projected by the entire pool of climate models, in term of both means and extremes; 2.
To devise procedures/methodologies for evaluating the skills of climate models in simulating the annual climatic cycle of the recent past; 3.
To select suitable climate model runs using the devised methodologies, based on their skills in simulating past climate, as well as on their ability to represent specific parts of the full spectrum of climate model projections; and 4.
To downscale and/or bias correct the selected GCMs (or GCM-RCM chains, using the selected GCMs as boundary conditions) through appropriate methods.

Study Area
In the current study, the climate change model selection procedure was carried out for the UIB, which is spread over the Hindu-Kush, Karakorum and Himalayan ranges, and feeds the largest canal system in the world (Figure 1).This river basin is very important due to two main reasons: first, the irrigated agriculture of Pakistan overwhelmingly depends on the inputs from this river basin; and second, the region is probably a climate change hot-spot [8,9], with an extremely uncertain future hydro-climatology.The future scenario data from the selected models are intended to be used, after downscaling and bias correction, as input to the SWAT hydrological model [10] for quantifying possible climate change impacts on the hydrological dynamics of the basin.
Climatic variables are usually strongly influenced by topographic altitude.Thus, the northern valley floors of the UIB are arid and warm, with an annual precipitation of only 100-200 mm.These totals increase to 600 mm at 4400 m altitude, and glaciological studies suggest annual accumulation rates of 1500-2000 mm at height of 5500 m [11].The UIB draws more than 50% of its water from melting of seasonal and permanent snow cover in the Himalaya, Karakoram and the Hindu Kush (HKH) mountains [5,[12][13][14][15].A rise in temperature in the UIB will, therefore, result in elevated melt rates with huge impacts on the timing and magnitude of the generated flows.This will not only lead to a higher average stream flow, but also to an increase in the occurrence and magnitude of extremes, especially during high-precipitation events [16].There is also the possibility that the peak flows may shift to earlier months or other seasons, with a rise in temperature [5] in the UIB.

Study Area
In the current study, the climate change model selection procedure was carried out for the UIB, which is spread over the Hindu-Kush, Karakorum and Himalayan ranges, and feeds the largest canal system in the world (Figure 1).This river basin is very important due to two main reasons: first, the irrigated agriculture of Pakistan overwhelmingly depends on the inputs from this river basin; and second, the region is probably a climate change hot-spot [8,9], with an extremely uncertain future hydro-climatology.The future scenario data from the selected models are intended to be used, after downscaling and bias correction, as input to the SWAT hydrological model [10] for quantifying possible climate change impacts on the hydrological dynamics of the basin.
Climatic variables are usually strongly influenced by topographic altitude.Thus, the northern valley floors of the UIB are arid and warm, with an annual precipitation of only 100-200 mm.These totals increase to 600 mm at 4400 m altitude, and glaciological studies suggest annual accumulation rates of 1500-2000 mm at height of 5500 m [11].The UIB draws more than 50% of its water from melting of seasonal and permanent snow cover in the Himalaya, Karakoram and the Hindu Kush (HKH) mountains [5,[12][13][14][15].A rise in temperature in the UIB will, therefore, result in elevated melt rates with huge impacts on the timing and magnitude of the generated flows.This will not only lead to a higher average stream flow, but also to an increase in the occurrence and magnitude of extremes, especially during high-precipitation events [16].There is also the possibility that the peak flows may shift to earlier months or other seasons, with a rise in temperature [5] in the UIB.All these facts make UIB a very sensitive region to possible climate change, and even, according to some [17], a climate-change "hotspot".However, despite the necessity of intensified investigations on different aspects of climate change and its possible implications, the task is hindered by the harshness of the environment and the unavailability of representative data.The climatic data available in the UIB lacks suitable coverage, since the in situ meteorological observations in the UIB are sparse and mostly taken at valley stations.Furthermore, the complex orography of the UIB region also affects the amounts, spatial patterns and seasonality of the precipitation.Therefore, neither the All these facts make UIB a very sensitive region to possible climate change, and even, according to some [17], a climate-change "hotspot".However, despite the necessity of intensified investigations on different aspects of climate change and its possible implications, the task is hindered by the harshness of the environment and the unavailability of representative data.The climatic data available in the UIB lacks suitable coverage, since the in situ meteorological observations in the UIB are sparse and mostly taken at valley stations.Furthermore, the complex orography of the UIB region also affects the amounts, spatial patterns and seasonality of the precipitation.Therefore, neither the sparsely observed station data and gridded data products based on them, nor the sensors-based data, fully represent the precipitation regime of the region [6].

GCM Outputs
In the IPCC 5th assessment report, four representative concentration pathways (RCPs) are normally used as a basis for future climate modelling: one very high baseline emission scenario (RCP8.5),two medium stabilization scenarios (RCP4.5 and RCP6) and one mitigated scenario (RCP2.6)(Table 1).Source: [18][19][20] The current study intended to include emission scenarios, covering a wider range of Radiative forcing and future temperature anomaly, while remaining close to the reality and considering RCP's showing minimum differences with the 2005 onwards actual observed CO 2 emission trend and growth rates.Keeping these prerequisites in mind, out of the four options: RCP2.6 was not considered in the current selection as it seemed to be the least likely [18,21] and the mitigation effort implied by this RCP, is unfeasible in the current circumstances [22,23], because it needs a sustained global CO 2 mitigation rate of around 3% per year, not a likely prospect, at least in the near future [20].
Out of the remaining three RCPs, the high baseline emission scenario (RCP8.5)and one medium-stabilization scenario (RCP4.5)were selected for the current study.The RCP8.5 was included because it covers the higher end of radiative forcing, as well as the temperature change, and it is also in line with the observed trend of around 3% in the average annual CO 2 emission growth rates for 2005-2012 [20,23].
For the medium-stabilization scenarios, both RCP4.5 and RCP6 are equally acceptable, but due to time constraints, and because RCP4.5 shows a better match (≈1.5%) of the trends of the average annual CO 2 emission growth rates for the period of 2005-2012 than RCP6 (≈1.0%) [20,23], RCP4.5 was picked along RCP8.5 for the GCM selection procedure.
Additionally, in the current study, only the available GCM runs for the ensemble member r1p1i1 in the CMIP5 repository [24] are included in the initial list.This is done so as to keep open the possibility of using dynamically downscaled projections (driven by the selected GCMs as boundary conditions) by Regional Climate Models (RCMs), which in most cases have utilized boundary conditions from the ensemble member r1p1i1 of the GCMs.
In the current study, a total number of 42 available model runs (ensemble member r1p1i1) are evaluated for RCP4.5 and of 39 for RCP8.5.

Extremes Indices
For the assessment of model runs for extremes, the ETCCDI extremes indices are utilized.The annual extremes of the daily CMIP5 data were acquired from the ETCCDI extremes indices archive [25,26], provided at the "Canadian Centre for Climate Modelling and Analysis".This data was indirectly obtained and downloaded through the "KNMI Climate Explorer", which is a web-based research tool to investigate climate and climate change.

Observed Data
The climate station network in the UIB has historically been comprised of only a few low-altitude, valley-based stations.Although the number of in situ observational points has increased since the mid-nineties, with the installations of a few higher altitude automatic weather stations, the coverage is still very thin, and the data is often not very representative, especially for different elevation zones.Similarly, while most of the weather stations have become operational after the mid-nineties, long-term data is a rare commodity and is only available at limited locations.
Similarly, owing to the complex orography of the UIB region and to the co-action of different hydro-climatic regimes, neither the sparse observed station data or the gridded data products based on them, nor the sensor-based climatic datasets fully represent the precipitation regime of the region [6,16,27,28].Several studies have pointed out that precipitation and other climatic variables in the HKH region exhibit large changes over short distances and considerable vertical gradients [11,[29][30][31][32][33][34].
In the absence of long-term climate data with acceptable representation of the UIB climate, most climate-change studies have relied on either the very thin climatic observation network records or the gridded datasets based on them.In all these cases, either the data have acceptable quality, but shorter duration, or they have huge biases, especially, in the case of precipitation in regions with higher altitudes.These biases are further amplified when this data is used as a reference for bias correction or downscaling of climate projections, making the results questionable.
In the current study, therefore, a new long-term climate dataset was prepared (Figure 2).The work related to this new long-term gridded data product [7] is not included in this paper, but we utilized this new dataset instead of the readily available global or regional gridded historical climate datasets, for bias correction, downscaling and assessment of the reliability of climate models for the simulation of the past climate in the region.
These gridded precipitation and temperature data are derived, based on all the available in situ observations available in the UIB, through reconstruction for the periods before the mid-nineties, interpolation and correction for the orography and elevation-induced effects guided by available data for runoff, actual evapotranspiration and glacier mass-balance [7].

Selection and Shortlisting of GCMs/RCMs
The full spectrum of GCM projections is wide, with large uncertainties attached [35][36][37], and it cascades to even a larger spectrum when downscaled or translated into possible impacts.Furthermore, the available future projections differ vastly from each other and may range from very wet to drier or very warm to colder future climates, so that the models can be categorized as representing either Warm-Wet, Warm-Dry, Cold-Wet and Cold-Dry corners of the full spectrum, in addition to the projections which are around the median tendency of future model projections.

Selection and Shortlisting of GCMs/RCMs
The full spectrum of GCM projections is wide, with large uncertainties attached [35][36][37], and it cascades to even a larger spectrum when downscaled or translated into possible impacts.Furthermore, the available future projections differ vastly from each other and may range from very wet to drier or very warm to colder future climates, so that the models can be categorized as representing either Warm-Wet, Warm-Dry, Cold-Wet and Cold-Dry corners of the full spectrum, in addition to the projections which are around the median tendency of future model projections.
These issues have led to diverse views on how to select or use these climate model projections, or even whether these climate models or their downscaled outputs should explicitly be used at all, or should only be indirectly used, instead, as guides to generate a range of plausible scenarios more suited for targeted impact studies and practical adaptation planning [38].
In mountainous regions, such as the Upper Indus Basin (UIB), the issue of how to proceed with the climate change impact studies becomes more complicated, because not only may the uncertainties shown by the climate models for these regions be even greater [4,39], but also because of the lower margin for error, as the lives and livelihood of millions of people depend purely on the water resources generated in these basins.These issues have led to diverse views on how to select or use these climate model projections, or even whether these climate models or their downscaled outputs should explicitly be used at all, or should only be indirectly used, instead, as guides to generate a range of plausible scenarios more suited for targeted impact studies and practical adaptation planning [38].
In mountainous regions, such as the Upper Indus Basin (UIB), the issue of how to proceed with the climate change impact studies becomes more complicated, because not only may the uncertainties shown by the climate models for these regions be even greater [4,39], but also because of the lower margin for error, as the lives and livelihood of millions of people depend purely on the water resources generated in these basins.
The usual approach of selecting results of a certain model or group of models or opting for a scenario with the mean trend of future projections may not be practical, as the full range of possible future climatic conditions needs to be covered in order to assess the full range of expected impacts required for climate adaptation needs.
As mentioned earlier, the selection of GCMs can be done following different approaches and may be based on a single criterion or a set of criteria.These approaches may include criteria such as: the total amount of change in the mean and/or an extreme of a projected climate variable; the success of a GCM in simulating the past climate for means or extremes; or maybe the skill in presenting the same pattern of tele-connections that drive the climate of the study region, and so on.
The current study adopted a combination of some of these approaches and applied a step-wise shortlisting of climate models based on a range of projected change in the (a) mean, (b) extremes, and (c) skill in reproducing the past climate.As the aim was to arrive at a limited number of models that Climate 2018, 6, 89 7 of 24 can represent not only all the possible futures as projected by the entire pool of climate models, but also changes in climatic extremes, so that the selected model runs can provide representation of the full spectrum of future climate projections by GCMs in terms of change in mean, as well as extremes.In other words, for each selected RCP, we intended to filter and select five climate model runs, each representing one the four corners of the spectrum or the median tendencies.

Shortlisting Based on Changes in the Means
As a first step, the total number of available model runs (ensemble member r1p1i1) for RCP4.5 (42) and RCP8.5, (39) were evaluated and shortlisted based on the change presented by them, in terms of the mean annual precipitation sum (∆P) and the mean air temperature (∆T), averaged across the UIB, between the simulated reference period historical data  and the late 21st-century projected data (2071-2100).The calculations were done using the web-based application "Climate Explorer" managed by the Royal Netherlands Meteorological Institute (KNMI) (http://climexp.knmi.nl).
As our intention was to identify fewer model runs that best represent the four corners of the full spectrum, as well as the central and middle tendencies, we first determined the 10th, 50th and 90th percentile values of ∆P and ∆T for the entire ensemble considered for each RCP, to explore the extent of the full spectrum of the projected changes in temperature and precipitation under that RCP.This was followed by determining the four (4) closest projections to each of the corners, as well as the center of the spectrum.The total number of shortlisted model runs for each of the two RCPs then amounted to 20.
Details of the different parts of the full spectrum considered during this study are as follows: • the Dry-Cold corner, represented by the 10th percentile ∆P as well as 10th percentile value of ∆T; The identification of the closest model runs to any corner point was done according to the procedure suggested by [19].It should be noted that 10th and 90th percentiles were selected as the central points of the corners, rather than the maximum or minimum values, in order to avoid selection of any outlier projections.

Ranking Based on Changes in Climate Extremes
To ascertain that preference will be given to those climate model runs that represent the full range of projected change in extremes, all 20 shortlisted model runs for each of RCP4.5 and RCP8.5, were further scrutinized and ranked based on their projected changes in climatic extremes.To that end, the ETCCDI indices [25] (Table 2) were used to evaluate changes in climatic extremes for air temperature, as well as precipitation.For the former, changes in the extremes were ranked and evaluated based on two indices-the warm spell duration index (WSDI), and the cold spell duration index (CSDI)-while for the latter, consecutive dry days (CDD) and the precipitation due to extremely wet days (R99pTOT) were considered.
To keep the work manageable, we only analyzed four indices in total, two indices to represent changes in precipitation extremes and two for changes in temperature extremes.Furthermore, as the intended use of the selected climate model ensemble was to force the hydrological model for assessing climate change impacts on both flows and extremes, we chose the four most obvious indicators of precipitation and temperature extremes.
The R99pTOT (precipitation due to extremely wet days (>99th percentile)) and CDD (consecutive dry days: maximum length of dry spell (P < 1 mm)) are appropriate indicators for precipitation extremes and suitable for assessment of associated hydrological extremes.R99pTOT is an important indicator for wet spells in terms of their length and magnitude, which are both key influencing factors in shaping extreme hydrological events (floods and high flows).Similarly, CDD is an important indicator for dry spells that can provide a good opportunity for assessment of the associated low flow episodes.The two temperature-related extreme indices used, i.e., WSDI (count of days in a span of at least 6 days where TX > 90th percentile) and CSDI (count of days in a span of at least 6 days where TN < 10th percentile), seemed best suited for their effects on evapotranspiration and cryospheric processes.The snow and glacier melt/accumulation, as well as evapotranspiration dynamics, are very important in the highly glacierized study area.
The changes in these indices, averaged over the UIB and over 30 years, between the reference period  and the late 21st century projections (2071-2100), were calculated using the database available at the ETCCDI extremes indices archive (http://climexp.knmi.nl),constructed by [26,40].
Only the relevant index for the air temperature or for the precipitation was considered for each of the previously selected group of models (a set of four) initially shortlisted models for each corner or the center), so that for the models in the Wet-Warm corner, only the R99pTOT index for precipitation and the WSDI index for temperature were considered, because they were the only relevant indices, as R99pTOT indicates extreme precipitation events, while WSDI indicates warm spells (Table 2).The other two indices, i.e., CDD and CSDI, were not considered in this case; however, they were the only indices considered for models in the Dry-Cold corner.For each corner, the relevant indices were given scores based on the ratio of the extreme index to the mean of that index, for all four models in a corner.For example, in the Wet-Warm corner, the % change in R99pTOT for a single model is divided by the mean of the % change in R99pTOT for all four models in that corner.The same procedure was applied for WSDI and, finally, both scores were averaged to obtain a final score.For each of the extreme indices, a weighted rank/skill score (Sk EI ) was calculated, with the highest value among the group getting the highest weighted rank/skill score of 1, and the others getting a rank according to their difference from this highest value, i.e., where Sk is the weighted rank for the specific extreme index EI, h denotes the highest index value in a group, and t denotes the target index to be ranked.Similarly, in the case of the change in means, i.e., ∆T ( • C) and ∆P (%), the ranking (Sk m ) was done based on the difference ∆T ( • C) or ∆P (%) shown by each member with the percentile value relevant to that group, The models were also evaluated with respect to their skill at simulating the past climate during the reference period .The selected climate model simulations were compared to the reference temperature and the precipitation gridded dataset [7] and were assigned skill scores.We did not use the same method for assigning skill score to temperature and precipitation.For assessing the performance of models in simulating past temperature, the method we applied was adopted from Perkins et al. [41].In this method, the skill score for temperature is calculated based on the identification of similarities between PDFs of modelled data and the observed reference data.A metric is generated to calculate the cumulative minimum value of each binned value for the two distributions, which represent the common area between two PDFs.This skill score (Sk Tmp ) can be expressed as follows: where n is the number of bins used to calculate the PDF, Z CM is the frequency of values in a given bin from the model while Z CM is the frequency of values in a given bin from the observed data.This skill score is 1, when there is a perfect match between simulated and the observed data, while a score of 0 means no similarities at all.The number of bins used in this study to generate the PDFs was 50.
In the case of precipitation, the skill score is calculated by a method proposed by [42] as the product of five skill functions, each assessing similarities between modelled and observed data, while covering different aspects of precipitation behavior.These five skill score functions for a particular model j are listed below: where A CMj and A Obs are the areas below the simulated (climate model j) and the observed precipitation cumulative density function (PDF) curves, respectively, and A+ and A− are the fractional areas over (+) and under (−) the 50th percentile.P denotes the average annual precipitation over UIB and σ is the standard deviation of the probability distribution function.Each of the above factors is intended to cover different aspects of probability distribution characteristics of the climate models, so that the distribution as a whole is taken into account through the mean and the total area (Equations ( 4) and ( 7)), the smaller and higher precipitation amounts are accounted for, through the 50th-percentile limit (Equations ( 5) and ( 6)), while the shape of the distribution is defined through the variance (Equation ( 8)).
These five factors are multiplied together to yield a single final skill score (Sk Prec ) for precipitation estimated by each model j: As a final step, all the rankings/scores, based on the changes in the means and in the extremes, as well as the skill scores for reproducing reference temperature and precipitation, are multiplied together to get the final overall skill or rank as follows: (10) Under this skill score, a higher value indicates better performance, while a lower value indicates otherwise.These skill scores can be further translated to a simple ranking of 1 to 4 for each group of climate models.
The climate model selection procedure adopted in this study is in line with the approach and methods suggested by [21,41,42], although with certain modifications in the evaluation criteria.For assessing the performance of models in simulating past temperature, the method applied is adopted from [41], while in the case of precipitation, guidance is taken from [42].A major difference from [21] in assessing model performances in simulating past climate is the use of a new long-term climate data set and an additional evaluation step for assessing model runs for their skill in reproducing the annual cycle of precipitation and temperature as well.

Downscaling and Bias Correction
After the shortlisting and ranking of the GCMs, the next step was to address the two primary issues inhibiting impact studies: firstly, the coarse spatial scales represented by the GCM may not be as fine as required by regional-and local-scale environmental modelling or impact studies; and secondly, the GCM raw outputs, or their downscaled versions, are deemed to contain systematic errors (bias) of certain magnitude, relative to the observational data, and therefore need post-processing by correcting it with and towards observations prior their its use in environmental modelling or impact studies.
The downscaling can be done either by applying statistical downscaling methods or through dynamical downscaling via application of a regional climate model (RCM).As dynamical downscaling was too demanding in terms of time and computational resource requirements, we decided only to explore whether any dynamically downscaled RCM projections were available for the already shortlisted and ranked GCMs.The Coordinated Regional Downscaling Experiment (CORDEX) has generated fine-scale climate projections for different regions of the world, of which the CORDEX-South Asia experiments cover the UIB region.We found that CORDEX-RCM model projections were available for four of the selected GCMs at the 1st rank and one GCM at the 2nd rank.These RCM projections provide dynamically downscaled data at a resolution of ~50 km for all of our selected GCMs.The data for the relevant GCM-RCM combinations were downloaded, but needed further downscaling, as the scale was still not fine enough, and also needed to undergo bias correction before further use in hydrological modelling.This downscaling and bias correction was achieved by the Distribution Mapping method (DM) [43], which was selected out of five different bias correction methods for the precipitation climate variable.These methods included: (1) Linear Scaling (LS); For the temperature, the selection was made after evaluating the performance of the following three bias correction methods: (1) Linear Scaling (LS); (2) Variance Scaling (VS); and (3) Distribution Mapping (DM).Further details of these methods can be found in [43].
The calibration and validation statistics, along with brief explanations, are provided as appendices (Supplementary Materials: Appendix A, Tables A1 and A2).

Shortlisting of Models: Changes in Climatic Means
The results of the initial shortlisting of the GCM model runs are given in Figures 3 and 4. In this step, only those GCM runs were retained which showed minimal difference with the 10th, 50th and 90th percentile values of ∆T ( • C) and ∆P (%), so that, for each RCP, we were left with sets of 4 GCM runs at each corner and 4 in the middle, while the remaining model runs were not processed any further.In this way, a total of 20 model runs were selected for each RCP.
It is worth mentioning that the range of projections for ∆T and ∆P for the RCP8.5 model pool was much larger than for the RCP4.5 model pool.For the latter, more extreme RCP, ∆P ranges from −5.42% to 19.56%, and ∆T ranges from 1.26 • C to 5.41 • C; while for the former (RCP4.5),these ranges are much higher, with ∆P ranging between −12.01%and 35.12% and ∆T between 1.48 • C and 8.57 • C.
The shortlisted GCM runs were also ranked according to their differences with the 10th, 50th or 90th percentile values in the respective corner or center.This ranking was intended for use in the final selection step, so that those model runs which show closest representation of the group of models or type of scenarios (Warm-Wet, Warm-Dry, Cold-Wet, Cold-Dry or the Median) get preference during the final selection.
It should be noted that the term Cold used in the "Wet-Cold" and "Dry-Cold" scenarios does not mean that the future temperatures will be colder than those of the reference period, but rather indicates that the warming will be less than that of the Warm scenarios.Similarly, the term Dry in the scenarios "Dry-Cold"and "Dry-Warm" is also only indicative of its comparative position relative to other climate models.

Ranking Based on Changes in Climatic Extremes
The 20 shortlisted model runs for each RCP were further scrutinized based on their projected changes in climatic extremes.The details of the projected changes in selected extreme indices are

Ranking Based on Changes in Climatic Extremes
The 20 shortlisted model runs for each RCP were further scrutinized based on their projected changes in climatic extremes.The details of the projected changes in selected extreme indices are given in Table 3.The darker colors indicate the higher values, while the lighter indicates lower values.These indices were given a weighted rank/score based on their difference from the highest value in the group of four model runs in a corner.
Similar to the rank assigned based on changes in the means, this ranking was also intended for use in the final selection step, so that the model runs, which show the largest changes in the extreme indices for each of the corner: Warm-Wet, Warm-Dry, Cold-Wet or Cold-Dry, get preference during the final selection.Unlike the four corners, evaluation based on the extreme indices was not carried out for the central or the mean scenario.
The ranking and scores for means and extreme indices, as well as the skill scores for simulating reference climate, are presented in Tables 4 and 5 for RCP4.5 and RCP8.5, respectively.
In most cases, the model run with the highest or the lowest changes in mean precipitation or temperature coincide with the highest change in relevant extreme index as well.
The index ∆ R99pTOT (%) was evaluated to represent the "Wet" scenarios, while the ∆ CDD (%) represented the "Dry" scenarios.Similarly, ∆ WSDI (%) was considered for the "Warm" scenarios, while ∆ CSDI (%) was considered for the "Cold" scenarios.In this way, a set of two (2) indices out of the four (4) were evaluated for each of the scenarios: Warm-Wet, Warm-Dry, Cold-Wet, and Cold-Dry.After checking the model runs for their projected changes in means and extreme indices, they were finally evaluated for their skill at reproducing the reference precipitation and temperature data.
The ranking for past performance utilized a new set of reference precipitation and temperature data [7], averaged over the UIB.The skill scores were calculated following the procedure of Section 4.1.3,and are presented in columns g and h in Tables 4 and 5.For most scenarios, the same models performed better than the others for both RCPs in simulating past climate.
Climate 2018, 6, 89 14 of 24 After allocating the skill score based on the past performance, the final skill scores and ranks were calculated by multiplying all the relevant skill scores allocated to each model run.The final ranks were allocated to each scenario, with the highest rank allotted to the model run with highest final skill score, and so on.
It is interesting to note that for the 4 scenarios, Warm-Dry, Cold-Wet, Cold-Dry and Median, for both RCPs, the same GCMs get the highest skill scores and ranks.The only exception is the Warm-Wet scenario, where different models top the ranking.In this scenario, for RCP4.5, the GCM "MIROC5" is in the top rank, followed by "CanESM2", while for RCP8.5, the ranking of these two GCMs is reversed.In the previous section, the step-wise shortlisting of the various climate models was based on the range of projected change in (a) mean, (b) extremes, and (c) skill in reproducing the past climate.Although the main aim of this approach was to combine the strengths of two different methodologies, i.e., the selection of the GCMs based on the properties of the full range of projections and the selection procedures based on past performance, certain limitations are unavoidable and need to be discussed.
First of all, the analysis considered only models selected based on the changes in the means and only the ensemble member r1p1i1, resulting in a reduced number of GCM runs for evaluation and possibly a smaller range of climatic extremes.This may also have led to possibly screening out models which may have had better past performance.
Similarly, another issue is concerned with the scale at which the method was applied.During the shortlisting step, and also during the evaluation of extreme indices, the projected changes or ETCCDI extremes indices were averaged over the entire UIB, which has the possibility of decreasing the spatial variation in projected changes.

Projected Changes in Temperature and Precipitation
The five ( 5) selected (CORDEX-SA) RCM outputs were further bias-corrected using the "distribution mapping technique" [43] for RCP4.5 and RCP8.5 for two sets of durations, i.e., mid-century (2041-2070) and end-century (2071-2100).Major properties of the downscaled projections are given in Table 7.The downscaled projections show changes in temperature ranging from 2.3 • C to 6.33 • C for RCP4.5 and of 2.92 • C to 9.0 • C for RCP8.5.The downscaled and bias-corrected precipitation ranges from a minor increase of 2.2% for the drier scenarios to as high as 15.9% for the wet scenarios.Thus, both temperature and precipitation show increases, as do the extremes, since the probabilities of the wet days are projected to decrease, while the precipitation intensities are projected to increase unanimously by both RCPs.
The spatial distribution of the projected future changes for precipitation and temperature across the UIB also show certain distinct trends.Thus, the precipitation (Figure 5) over the mid-century (2041-2070), as well as the late century (2071-2100), reveals for all scenarios a remarkable decrease in the southeastern parts of the basin, but an increase in the northeastern parts.This decrease/increase is particularly intense for the "Dry-Warm" and the "Median" scenarios over the late 21st century.The spatial distribution of the projected changes for in temperature (Figure 6) also shows similarities across all scenarios, with the northern and northwestern parts of the basin exhibiting higher increases, while the eastern and southern parts experience a comparatively smaller temperature increase.
For RCP8.5, the projected temperature changes appear to be very high over the late 21st century and this occurs under all scenarios, especially for the "Warm" scenarios, with an almost uniform spread across the whole UIB.The projected temperature changes range for all RCPs and the two 20th-century periods from a minimum increase of 3.76 • C (NorESM1-M_RCA4, RCP4.5, Period: 2041-2070) to a maximum increase as high as 10.4 • C (IPSL-CM5A-MR_RCA4, RCP8.5 and period: 2071-2100).

Limitations of the Downscaling and Bias Correction Approach Adopted
The downscaling and bias correction approaches described in the previous sections, although adopted as the best available option considering the availability of time, resources and data, are still not without limitations.
Firstly, we opted for using GCM-RCM chains for the already shortlisted climate model runs where the RCM add further detail to global climate simulations and may provide regional-to local-scale information, but their outputs are still subject to inherited or new systematic errors and may therefore require a bias correction or further downscaling to a higher resolution, along with the fact that they may also produce quite a different response in terms of the temperature and precipitation change than the forcing GCMs.In fact, many authors recognized this affliction of GCM-RCM model chains with systematic errors (biases) [48][49][50][51][52][53][54][55], but still the representation of explicit atmospheric and surface processes and the level of spatial details provided by the RCMs [56] make them a better option for regions with high topographic variability, The selection procedure may in future, though, directly include the RCM runs instead of GCM, when an appropriate amount of RCM runs are available.
The systematic errors (biases) in GCM/RCM outputs make them unsuitable for certain uses, and application of different bias-correction methods has increasingly become a standard procedure, especially in climate change impact studies [57].
Although these bias correction approaches improve the agreement of climate model output with observations, and therefore narrows the uncertainty range of predictions and simulations, they do so without a sound physical basis [57].At the same time, a growing number of authors are showing reservations on the use of bias correction methods for a variety of reasons.Some argue that the main assumption behind bias correction approaches is the stationarity of the correction parameters, which is not realistic and may not be the case, especially under climate change [57,58].While others believe that as bias correction cannot overcome major climate model errors, inexperienced application might result in ill-informed adaptation decisions [59].Despite being critical of bias correction methods, many of these authors e.g., [57,59,60] acknowledge the need for some type of bias correction and recommend different precautions to avoid any pitfalls.

Conclusions
It is essential to have representative future climate projections of appropriate quality for climate change impact studies, especially in the water resource sector.Despite the availability of an increasing number of GCM outputs in the CMIP5 archive and the on-going improvements in their process representations, issues of large uncertainties in their future climate predictions cannot be avoided.This situation, along with other factors, such as time, human resources or computational constraints, make it imperative to sort out the most appropriate individual GCM or small ensemble of GCMs for a more reliable assessment of climate change impacts.
The approach presented in the present study seeks the most suitable set of climate model runs, while considering not only the full ranges of projected changes in terms of means and extremes by different climate models, but also their skills in simulating the past climate in a reference period.
This selection procedure was applied for future climate projections over the Upper Indus Basin for two representative concentration pathways (RCPs), the RCP 4.5 and RCP 8.5.All available model runs for the r1p1i1 ensemble member of each GCM in the CMIP5 repository were included in the initial list.The total number of model runs available for RCP4.5 was 42, and 39 for RCP8.5.
Based on the huge uncertainties reported in the GCM runs for the UIB, all possible extreme future scenarios (Wet-Warm, Wet-Cold, Dry-Warm, Dry-Cold) were considered, in addition to the selection of GCMs representing the mean future climate change, with respect to both changes in the projected means and the extremes.This procedure made it possible to arrive at a limited number of climate models, from which the final selection was performed by assigning ranks based on the weighted score for each of the mentioned selection criteria.
Finally, the precipitation and temperature time series of the selected GCM model runs were bias corrected and further downscaled to the scale of the reference data by means of a distribution mapping technique.The ensembles of the selected GCM runs for RCP4.5 and RCP8.5 scenarios show that the uncertainty of future climate in the study region is very large for the raw data, as well as their downscaled versions.
The downscaled projections indicate increases of temperature ranging between 2.3 • C and 9.0 • C and changes in precipitation that range from a slight annual increase of 2.2% under the drier scenarios, to as high as 15.9% for the wet scenarios.Thus, for both temperature and precipitation, the future projections under all scenarios and both RCP's only show increases in the mean annual values, with no negative trend.Moreover, for all scenarios, the future precipitation is projected to be more extreme, as the probability of wet days decreases, while at the same time, the precipitation intensities will increase.
The spatial distribution of the downscaled predictors, namely, the precipitation, also shows distinct patterns across the UIB, such that this variable shows for all time periods/scenarios considered a distinct decrease in the southeastern parts, but an increase in the northeastern parts of the basin.This decrease/increase is particularly intense for the "Dry-Warm" and the "Median" scenarios over the late 21st century.
Overall, the future climate of the UIB region remains very uncertain, which justifies the selection procedure proposed here to arrive at a wider range of possible climate scenarios that can then be further utilized and translated into a wider spectrum of climate change impact scenarios.

Climate 2018, 6 , 25 Figure 3 .
Figure 3. Projected changes in mean air temperature (ΔT) and annual precipitation sum (ΔP) between 2071 and 2100 and 1971 and 2000 for all included RCP4.5 GCM runs.Blue crosses indicate the 10th, 50th and 90th percentile values for ΔT and ΔP.The model runs shortlisted during this step are indicated in red color.

Figure 3 .
Figure 3. Projected changes in mean air temperature (∆T) and annual precipitation sum (∆P) between 2071 and 2100 and 1971 and 2000 for all included RCP4.5 GCM runs.Blue crosses indicate the 10th, 50th and 90th percentile values for ∆T and ∆P.The model runs shortlisted during this step are indicated in red color.

Figure 3 .
Figure 3. Projected changes in mean air temperature (ΔT) and annual precipitation sum (ΔP) between 2071 and 2100 and 1971 and 2000 for all included RCP4.5 GCM runs.Blue crosses indicate the 10th, 50th and 90th percentile values for ΔT and ΔP.The model runs shortlisted during this step are indicated in red color.

Table 4 .
Weighted ranks for all shortlisted RCP4.5 GCM runs based on change in means (e and f), change in extremes (a, b, c and d) and their skill scores for simulating reference precipitation and air temperature(g and h).

Climate 2018, 6 , 25 Figure 5 .
Figure 5. Spatial distribution of projected precipitation change across the UIB over the mid (2014-2070) and the late-(2071-2100) 21st century for 5 models and 2 RCPs.The figure is arranged in a tabular form where the 1st and 2nd column represent projected change in precipitation for RCP 4.5, for the mid-century (2014-2070) and the late-century (2071-2100), respectively, while the 3rd and 4th columns show the projected change in the mid-century and the late-century precipitation for RCP 8.5, respectively.The rows represent the climate models used.

Figure 5 .
Figure 5. Spatial distribution of projected precipitation change across the UIB over the mid (2014-2070) and the late-(2071-2100) 21st century for 5 models and 2 RCPs.The figure is arranged in a tabular form where the 1st and 2nd column represent projected change in precipitation for RCP 4.5, for the mid-century (2014-2070) and the late-century (2071-2100), respectively, while the 3rd and 4th columns show the projected change in the mid-century and the late-century precipitation for RCP 8.5, respectively.The rows represent the climate models used.

Figure 5 .
Figure 5. Spatial distribution of projected precipitation change across the UIB over the mid (2014-2070) and the late-(2071-2100) 21st century for 5 models and 2 RCPs.The figure is arranged in a tabular form where the 1st and 2nd column represent projected change in precipitation for RCP 4.5, for the mid-century (2014-2070) and the late-century (2071-2100), respectively, while the 3rd and 4th columns show the projected change in the mid-century and the late-century precipitation for RCP 8.5, respectively.The rows represent the climate models used.

Figure 6 .
Figure 6.Similar to Figure 5, but for temperature changes.Figure 6.Similar to Figure 5, but for temperature changes.

Figure 6 .
Figure 6.Similar to Figure 5, but for temperature changes.Figure 6.Similar to Figure 5, but for temperature changes.

Table 2 .
List of ETCCDI extreme indices used during the GCM selection procedure.

Table 7 .
Future precipitation and temperature projections from 5 GCM models, 2 RCPs and 2 periods.