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

Abstract

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.

1. 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:

• 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;
• To devise procedures/methodologies for evaluating the skills of climate models in simulating the annual climatic cycle of the recent past;
• 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
• To downscale and/or bias correct the selected GCMs (or GCM-RCM chains, using the selected GCMs as boundary conditions) through appropriate methods.

2. Study Area and Data Used

2.1. 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 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].

2.2. Data Used

2.2.1. 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. Representative concentration pathways (RCPs), their radiative forcing, emissions (CO₂ equivalent and growth rate %) and temperature increase.

RCP	Radiative Forcing	CO2 Equiv. (ppm)	Temperature Increase (°C)	Pathway	CO2 Growth Rate (%)
RCP8.5	8.5 Wm−2 in 2100	1370	4.9	Rising	≈2.5
RCP6.0	6 Wm−2 post 2100	850	3.0	Stabilization without overshoot	≈1
RCP4.5	4.5 Wm−2 post 2100	650	2.4	Stabilization without overshoot	≈1.5
RCP2.6 (RCP3PD)	3 Wm−2 before 2100, declining to 2.6 Wm−2 by 2100	490	1.5	Peak and decline	≈1.6
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₂ 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₂ 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₂ 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₂ 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.

2.2.2. 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.

2.2.3. 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].

2.2.4. RCM Outputs

Five CORDEX-SA experiments (Table 6), including IPSL-CM5A-MR_RCA4, MPI-ESM-LR_RCA4, NorESM1-M_RCA4, Can ESM2_RegCM4-4, and GFDL-ESM2M_RCA4, were downscaled and bias corrected. These five GCMs have been dynamically downscaled by CORDEX, using two different RCMs (RCA4 and RegCM4). Their RCM outputs are at considerably finer scale (0.44°) then the source GCMs.

3. Methods

3.1. 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.

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 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.

3.1.1. 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 (1976–2005) 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 Dry-Warm corner, represented by the 10th percentile ΔP but the 90th percentile value of ΔT;
• the Wet-Cold corner, represented by the 90th percentile ΔP and the 10th percentile value of ΔT;
• the Wet-Warm corner, represented by the 90th percentile values for both ΔP as well as ΔT; and finally
• the median projected future climate, represented by the 50th percentile values of both ΔP and Δ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.

3.1.2. 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. List of ETCCDI extreme indices used during the GCM selection procedure.

Climate Variable	ETCCDI Index	Description of the ETCCDI Index
Precipitation	R99pTOT	Precipitation due to extremely wet days (>99th percentile)
	CDD	Consecutive dry days: maximum length of dry spell (P < 1 mm)
Air Temperature	WSDI	Warm spell duration index: count of days in a span of at least 6 days where TX > 90th percentile
	CSDI	Cold spell duration index: count of days in a span of at least 6 days where TN < 10th percentile

) 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 (1976–2005) 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.,

$$S{k}_{EI}=1-\frac{E{I}_h-E{I}_t}{E{I}_h}$$

(1)

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,

$$S{k}_m=1-\frac{{(\Delta Tor\Delta P)}_{{10,50\mathrm{or}90}^{\mathrm{th}}\mathrm{percentile}}-{(\Delta Tor\Delta P)}_{t\arg et}}{{(\Delta Tor\Delta P)}_{{10,50\mathrm{or}90}^{\mathrm{th}}\mathrm{percentile}}}$$

(2)

3.1.3. Ranking Based on Skill in Reproducing the Reference Climate

The models were also evaluated with respect to their skill at simulating the past climate during the reference period (1976–2005). 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 ($S{k}_{Tmp}$) can be expressed as follows:

$$S{k}_{Tmp}={\displaystyle \sum _1^n\mathrm{minimum}(}{Z}_{CM},Z_{Obs})$$

(3)

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:

$$f_{1j}=1-{\left(\frac{|A_{CMj}-A_{Obs}|}{2\cdot A_{Obs}}\right)}^{0.5}$$

(4)

$$f_{2j}=1-{\left(\frac{|A_{CMj}^{+}-A_{Obs}^{+}|}{2\cdot A_{Obs}^{+}}\right)}^{0.5}$$

(5)

$$f_{3j}=1-{\left(\frac{|A_{CMj}^{-}-A_{Obs}^{-}|}{2\cdot A_{Obs}^{-}}\right)}^{0.5}$$

(6)

$$f_{4j}=1-{\left(\frac{|\overline{P_{CMj}}-\overline{P_{Obs}}|}{2\cdot \overline{P_{Obs}}}\right)}^{0.5}$$

(7)

$$f_{5j}=1-{\left(\frac{|{\sigma }_{CMj}-{\sigma }_{Obs}|}{2\cdot {\sigma }_{Obs}}\right)}^{0.5}$$

(8)

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:

$$S{k}_{Prec}=f_{1j}\cdot f_{2j}\cdot f_{3j}\cdot f_{4j}\cdot f_{5j}$$

(9)

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:

$$FinalSkillScore=S{k}_{E{I}_1}\cdot S{k}_{E{I}_2}\cdot S{k}_{\Delta T}\cdot S{k}_{\Delta P}\cdot S{k}_{Temp}\cdot S{k}_{\mathrm{Pr}ec}$$

(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.

3.2. 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); (2) Local Intensity Scaling (LIS); (3) Power Transformation (PT); (4) Distribution Mapping (DM); and (5) Distribution Mapping followed by Intensity and frequency Scaling (DM-IS).

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).

4. Results

4.1. Selection of Climate Models

4.1.1. Shortlisting of Models: Changes in Climatic Means

The results of the initial shortlisting of the GCM model runs are given in Figure 3 and Figure 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.

4.1.2. 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. Percentage change in ETCCDI indices (R99pTOT, CDD, WSDI, and CSDI) along with changes in mean precipitation (ΔP) and temperature (ΔT), for all corners/scenarios (Warm-Wet, Warm-Dry, Cold-Wet and Cold-Dry) and both RCPs (RCP4.5 and RCP8.5).

Projection	Model	Δ R99pTOT (%)	Δ CDD (%)	Δ WSDI (%)	Δ CSDI (%)	ΔP (%) (%)	ΔT (°C)
RCP 4.5
Wet-Warm	CanESM2	29.0	−7.2	814	−96.2	13.0	3.6
	HadGEM2-ES	28.6	12.5	1002	−98.7	3.7	3.6
	MIROC5	76.4	−8.8	938	−96.3	12.1	4.0
	MIROC-ESM-CHEM	19.8	2.2	611	−89.9	6.7	3.8
Wet-Cold	bcc-csm1-1-m	45.3	−1.0	298	−87.6	5.0	2.2
	GFDL-ESM2M	42.4	−4.9	202	−61.6	4.9	1.8
	IPSL-CM5B-LR	32.2	−11.7	293	−81.6	5.2	2.0
	MRI-CGCM3	59.6	−7.5	471	−89.8	9.0	2.2
Dry-Warm	ACCESS1-0	46.4	0.9	656	−92.1	3.47	3.5
	CMCC-CMS	61.9	7.1	454	−89.8	−3.35	3.6
	IPSL-CM5A-MR	54.5	12.0	604	−90.2	1.28	3.9
	MIROC-ESM	26.8	1.8	718	−97.0	2.41	4.2
Dry-Cold	CCSM4	4.8	−0.8	323	−92.0	4.54	2.4
	GFDL-ESM2G	16.2	−0.1	373	−70.9	2.14	2.2
	inmcm4	2.0	4.3	216	−48.9	−5.29	4.2
	MPI-ESM-LR	42.3	17.7	406	−89.1	−5.76	2.8
RCP 8.5
Wet-Warm	CanESM2	101.7	−12.3	1181	−97.3	18.5	6.7
	GFDL-CM3	9.7	−5.0	1426	−100.0	9.0	8.6
	MIROC5	257.2	−13.4	1640	−98.5	35.1	6.2
	MIROC-ESM-CHEM	28.5	14.5	1314	−100.0	6.6	7.6
Wet-Cold	GFDL-ESM2G	95.9	−1.0	668	−99.0	12.6	4.7
	GFDL-ESM2M	72.9	−3.1	1696	−95.5	13.5	4.2
	CNRM-CM5	68.6	−3.5	638	−96.1	14.1	4.1
	MRI-CGCM3	195.5	−12.6	1309	−98.4	24.1	4.6
Dry-Warm	IPSL-CM5A-LR	94.5	23.3	1022	−97.6	−12.01	7.0
	IPSL-CM5A-MR	194.6	9.3	1358	−99.1	−2.95	7.1
	MIROC-ESM	9.5	4.4	1521	−100.0	0.06	7.3
	CMCC-CMS	143.9	18.8	985	−99.9	−2.26	6.0
Dry-Cold	MPI-ESM-LR	136.0	29.1	1067	−98.2	−4.49	5.2
	CCSM4	48.3	7.0	871	−99.5	0.86	4.6
	inmcm4	61.3	4.7	849	−85.9	1.48	4.1
	NorESM1-M	107.1	3.5	1010	−98.6	6.01	4.6

. 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

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).

Projection	Climate Model	a	b	c	d	e	f	g	h	Final Skill Score (a*b*c*d*e*f*g*h*10)	Final Rank
		Weighted Rank Δ R99pTOT (%)	Weighted Rank Δ CDD (%)	Weighted Rank Δ WSDI (%)	Weighted Rank Δ CSDI (%)	Weighted Rank Δ T (°C)	Weighted Rank Δ P (%)	Skill Score for Temperature (SkTmp)	Skill Score for Precipitation (SkPerc)
Wet-Warm	CanESM2	0.38		0.85		0.97	0.63	0.79	0.36	0.57	2
	HadGEM2-ES	0.37		1.00		0.94	0.39	0.73	0.29	0.29	3
	MIROC5	1.00		0.94		0.95	0.72	0.81	0.38	1.93	1
	MIROC-ESM-CHEM	0.26		0.61		1.00	0.71	0.71	0.25	0.20	4
Wet-Cold	bcc-csm1-1-m	0.74			0.93	0.77	0.53	0.71	0.40	0.80	3
	GFDL-ESM2M	0.71			0.75	0.97	0.52	0.79	0.41	0.88	2
	IPSL-CM5B-LR	0.54			1.00	0.94	0.55	0.71	0.14	0.28	4
	MRI-CGCM3	1.00			0.90	0.82	0.96	0.78	0.35	1.91	1
Dry-Warm	ACCESS1-0		0.07	0.91		0.93	0.78	0.77	0.35	0.14	2
	CMCC-CMS		0.59	0.75		0.97	0.81	0.70	0.31	0.75	4
	IPSL-CM5A-MR		1.00	1.00		0.97	0.59	0.79	0.33	1.52	1
	MIROC-ESM		0.15	0.81		0.87	0.92	0.74	0.22	0.16	3
Dry-Cold	CCSM4		0.05		1.00	0.64	0.34	0.79	0.26	0.02	3
	GFDL-ESM2G		0.00		0.77	0.58	0.56	0.75	0.32	0.00	4
	inmcm4		0.24		0.53	0.89	0.76	0.66	0.35	0.20	2
	MPI-ESM-LR		0.98		0.97	0.75	0.72	0.75	0.32	1.23	1
Mean	NorESM1-M					0.94	0.58	0.79	0.43	1.83	1
	bcc-csm1-1-m					0.76	0.70	0.76	0.44	1.76	2
	GFDL-ESM2G					0.87	0.85	0.75	0.32	1.77	2
	CMCC-CMS					0.56	0.20	0.70	0.31	0.24	4

and

Table 5. Similar to Table 4, but for RCP8.5.

Projection	Climate Model	a	b	c	d	e	f	g	h	Final Skill Score (a*b*c*d*e*f*g*h*10)	Final Rank
		Weighted Rank Δ R99pTOT (%)	Weighted Rank Δ CDD (%)	Weighted Rank Δ WSDI (%)	Weighted Rank Δ CSDI (%)	Weighted Rank Δ T (°C)	Weighted Rank Δ P (%)	Skill Score for Temperature (SkTmp)	Skill Score for Precipitation (SkPerc)
Wet-Warm	CanESM2	0.40		0.72		0.93	1.00	0.79	0.36	0.76	1
	GFDL-CM3	0.04		0.87		0.80	0.48	0.71	0.39	0.04	3
	MIROC5	1.00		1.00		0.86	0.10	0.81	0.38	0.27	2
	MIROC-ESM-CHEM	0.11		0.80		0.94	0.36	0.71	0.25	0.05	3
Wet-Cold	GFDL-ESM2G	0.49			0.97	0.89	0.68	0.75	0.32	0.69	4
	GFDL-ESM2M	0.50			1.00	1.00	0.73	0.79	0.34	0.97	2
	CNRM-CM5	0.35			1.00	0.98	0.76	0.81	0.41	0.87	3
	MRI-CGCM3	1.00			0.98	0.90	0.70	0.78	0.35	1.65	1
Dry-Warm	IPSL-CM5A-LR		1.00	0.67		0.98	0.41	0.79	0.33	0.71	2
	IPSL-CM5A-MR		0.40	0.89		1.00	0.84	0.78	0.34	0.79	1
	MIROC-ESM		0.19	1.00		0.98	0.60	0.74	0.22	0.18	4
	CMCC-CMS		0.81	0.65		0.84	0.78	0.70	0.31	0.74	2
Dry-Cold	MPI-ESM-LR		1.00		0.86	0.73	0.96	0.75	0.31	1.42	1
	NorESM1-M		0.12		0.85	0.64	0.01	0.79	0.50	0.00	2
	CCSM4		0.24		0.84	0.64	0.48	0.79	0.26	0.13	4
	inmcm4		0.16		1.00	0.57	0.42	0.66	0.34	0.09	2
Mean	NorESM1-ME					0.93	0.91	0.79	0.50	3.34	1
	GFDL-ESM2G					0.95	0.09	0.75	0.32	0.21	2
	CCSM4_r1i1p1					0.93	0.13	0.79	0.26	0.25	2
	bcc-csm1-1					0.92	0.86	0.77	0.28	1.70	4

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.

4.1.3. Ranking Based on Skill in Reproducing the Reference Climate

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 Table 4 and Table 5. For most scenarios, the same models performed better than the others for both RCPs in simulating past climate.

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.

4.1.4. Limitations of the Model Selection Procedure

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.

Additionally, the weighting of different skill scores in this study also differed for similar work, such as [21]. In our study, the final skill score was a combination of scores allocated for change in mean, change in extremes, and performance in reproducing past climate. This may have reduced the chances of selecting the climate model with the best past performance, but increased the chance of a better spread of scenarios over the entire range, while still taking past performance as a key factor in selections. Our model selection approach also assumes that all the evaluated model runs are independent of each other, which may not be the case, as some models use the same forcing and validation data or may share similar model codes [44,45].

Despite these limitations, the adopted approach made it possible for us to identify a limited number of model runs, representative of the full range of future projected means and extremes while giving due preference to models which perform better at simulating the reference climate.

4.2. Bias Correction and Downscaling of Future Climate Scenarios

The future climate projections of the selected climate models needed to be downscaled and corrected for biases, before further use in hydrological model simulation. Therefore, as a first option, all the dynamically downscaled climate projections available for UIB were checked for whether any Regional Climate Model (RCM) projections were available that have dynamically downscaled projections for the already shortlisted GCM at first or second position in the ranking. We found that, for both RCPs, the outputs of at least three (3) CORDEX-SA experiments were based on the GCMs ranked 1st in our study (IPSL-CM5A-MR_RCA4, MPI-ESM-LR_RCA4 and NorESM1-M_RCA4). The GSM CanESM2, which is at 2nd rank for RCP4.5 and at 1st for RCP8.5, has dynamically downscaled projections under CORDEX-SA experiments, CanESM2_RegCM4-4. The output of one CORDEX-SA experiment (GFDL-ESM2M_RCA4) was based on GFDL-ESM2M, which is ranked at 2nd position for both the RCPs in our study (Table 4 and Table 5).

It was decided to utilize the available dynamically downscaled data for our selected GCMs, ranked at 1st or 2nd positions. Therefore, five CORDEX-SA experiments (

Table 6. List of CORDEX South Asia experiments for RCP8.5 and their CMIP5-GCM forcing.

Nr	Scenario	Experiment		Driving AOGCM	RCM	RCM Description
		Name	Short Form
1	Wet-Warm	CanESM2_ RegCM4-4	CAN	CCCma-CanESM2 (1st)	RegCM4	Abdus Salam International Centre for Theoretical Physics (ICTP) Regional Climatic Model version 4 (RegCM4; [46])
2	Wet-Cold	GFDL-ESM2M_ RCA4	GFDL	NOAA-GFDL-GFDL-ESM2M (2nd)	RCA4	Rossby Centre regional atmospheric model version 4 (RCA4; [47])
3	Mean	NorESM1-M_ RCA4	NOR	Nor-ESM1-M (1st)
4	Dry-Cold	MPI-ESM-LR_ RCA4	MPI	MPI-ESM-LR (1st)
5	Dry-Warm	IPSL-CM5A-MR_ RCA4	IPSL	IPSL-CM5A-MR (1st)

), including IPSL-CM5A-MR_RCA4, MPI-ESM-LR_RCA4, NorESM1-M_RCA4, Can ESM2_RegCM4-4, and GFDL-ESM2M_RCA4 were selected for further processing and bias correction. These five GCMs were dynamically downscaled by CORDEX using two different RCMs (RCA4 and RegCM4). Their RCM outputs are at a considerably finer scale (0.44°) than the source GCMs.

4.3. 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. Future precipitation and temperature projections from 5 GCM models, 2 RCPs and 2 periods.

Model	Duration	Precipitation (mm)								Temperature (°C)
		RCP 4.5 Values and (change %)				RCP 8.5 Values and (change %)				RCP 4.5 Values and change			RCP 8.5 Values and change
		PCP (mm) (Av-An)	90th percentile (mm)	Probability-Wet Days (days)	Intensity-Wet Days (mm)	PCP (mm) (Av-An)	90th percentile (mm)	Probability-Wet Days (days)	Intensity-Wet Days (mm)	TMP (C°) (mean)	90th Percentile (C°)	10th Percentile (C°)	TMP (C°) (mean)	90th Percentile (C°)	10th Percentile (C°)
IPSL-CM5A-MR_ RCA4	41–70	539 (2.9%)	19.7 (36.8%)	107.4 (−7.8%)	5.0 (13.2%)	532 (1.7%)	19.5 (35.4%)	106.9 (−8.3%)	4.9 (11.8%)	5.52 (4.12)	17.6 (3.85)	−6.1 (4.93)	6.3 (4.91)	18.4 (4.66)	−5.3 (5.70)
	71–00	557 (6.2%)	20.7 (42.1%)	107.4 (−7.8%)	5.2 (17%)	502 (−4.2%)	20.5 (18%)	94.4 (−1.9%)	5.4 (22%)	7.7 (6.33)	19.1 (5.34)	−3.6 (7.44)	10.4 (9.0)	21.9 (8.17)	−1.2 (9.82)
MPI-ESM-LR_ RCA4	41–70	536 (2.3%)	18.6 (29.1%)	110.9 (−4.8%)	4.6 (5.2%)	535 (2.0%)	18.6 (29.1%)	108.5 (−6.9%)	4.7 (7.3%)	4.0 (2.64)	16.2 (2.44)	−8.1 (2.85)	4.5 (3.08)	16.7 (2.99)	−7.7 (3.27)
	71–00	537 (2.4%)	18.7 (41.1%)	109.1 (−6.4%)	4.7 (7.7%)	559 (6.7%)	20.3 (41.1%)	106.1 (−8.9%)	5.1 (15.5%)	5.5 (4.11)	17.4 (3.67)	−6.5 (4.48)	7.3 (5.86)	19.2 (5.51)	−4.9 (6.09)
NorESM1-M_ RCA4	41–70	536 (2.4%)	20.3 (46.1%)	109.0 (−6.4%)	4.9 (10.7%)	555 (6.0%)	21.1 (46.1%)	111.3 (−4.5%)	4.9 (12.3%)	3.8 (2.36)	16.8 (3.03)	−8.0 (3.02)	4.3 (2.92)	17.9 (4.2)	−7.5 (3.53)
	71–00	537 (2.5%)	20.3 (54.4%)	109.0 (−6.4%)	4.9 (12%)	548 (4.6%)	22.2 (54.4%)	107.0 (−8.2%)	5.2 (17.5%)	4.9 (3.50)	17.7 (3.99)	−6.76 (4.23)	6.6 (5.23)	19.9 (6.16)	−5.1 (5.93)
GFDL-ESM2M_ RCA4	41–70	540 (3.1%)	17.9 (42.6%)	111.9 (−4%)	4.7 (7%)	578 (10.4%)	20.5 (42.6%)	114.9 (−1.4%)	4.9 (11.1%)	3.8 (2.41)	16.0 (2.31)	−7.8 (3.22)	4.1 (2.73)	16.2 (2.43)	−7.0 (4.02)
	71–00	536 (2.2%)	19.4 (52.8%)	112.8 (−3.2%)	4.7 (5.7%)	612 (16.8%)	22.0 (52.8%)	114.7 (−1.5%)	5.2 (18.6%)	5.1 (3.70)	17.14 (3.42)	−6.4 (4.57)	6.6 (5.22)	18.8 (5.03)	−4.8 (6.17)
CanESM2_ RegCM4-4	41–70	560 (6.9%)	21.1 (43.8%)	119.6 (2.7%)	4.7 (6.4%)	557 (6.3%)	20.7 (43.8%)	115.6 (−0.8%)	4.6 (5.5%)	4.5 (3.14)	16.8 (3.08)	−7.8 (3.20)	4.9 (3.51)	16.9 (3.16)	−7.2 (3.75)
	71–00	607 (15.9%)	23.2 (51.6%)	117.2 (0.6%)	5.05 (14.8%)	590 (12.5%)	21.8 (51.6%)	114.9 (−1.3%)	5.0 (13.2%)	5.6 (4.24)	17.5 (3.8)	−6.5 (4.47)	7.4 (6.03)	20.0 (6.24)	−5.1 (5.89)
Observed	1976–2005	524	14.4	116.5	4.4	524.1	14.4	116.5	4.4	1.4	13.7	−11.0	1.4	13.7	−11.0

.

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).

4.4. 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.

5. 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.