Precipitation projections over the Indus River Basin of Pakistan for the 21st century using a statistical downscaling framework

We estimate future changes in precipitation over the entire Indus basin of Pakistan with a particular focus on the high‐elevation upper indus basin (UIB). A statistical downscaling approach is used. We consider the spatial variability of observed precipitation on seasonal scales. Large‐scale atmospheric patterns are employed for general circulation model (GCM) selection and subsequent projections. Firstly, we identify the precipitation governing predictors from ERA‐Interim reanalysis. We further quantify the robustness of governing predictors against other reanalysis datasets (ERA5 and NCEP‐NCAR‐II) to assess future projections' fidelity. We perform S‐mode Principal Component Analysis on predictor fields and compare loading patterns using Taylor diagrams to assess predictor correspondence between different reanalysis. Similarly, we compare ERA‐Interim variables with model‐simulated fields during the historical period to select better performing GCMs and quantify model uncertainty. The regional suitability of available GCMs in our study is also demonstrated. Ensemble (median) changes in regional precipitation derived through atmospheric fields show an elevation‐dependent response of the UIB at representative concentration pathway (RCP) scenarios RCP4.5 and RCP8.5, where increased precipitation will mostly fall at high elevations. However, the positive signals are more distinct during the winter and monsoon seasons, particularly over the central Karakoram. Meanwhile, a decrease in precipitation is robust during the pre‐monsoon period, particularly over the northwestern regions. These signals intensify and become more robust during 2071–2100 under RCP8.5, and the better‐performing models and signal‐to‐noise ratios further support this finding. The spatial patterns of projected changes suggest stronger (weaker) and further northward penetrating westerly systems during the winter (pre‐monsoon) season. Increased warming will also strengthen monsoon circulations, and these will penetrate further into the northwestern and trans‐Himalayan regions. The Lower Indus shows a mixed seasonal response that is more uncertain. The present analysis provides an alternative perspective to the ongoing research of assessing climate responses in complex regions.


| INTRODUCTION
The Indus River system originates within the massifs of the Hindukush, Karakoram, and Himalayans (HKH), which contains the largest non-polar cryosphere (e.g., Bocchiola et al., 2011;Soncini et al., 2015). The River system provides water, renewable energy, food security, and other ecological services to sustain millions of downstream population (Archer and Fowler, 2004). Complex processes involving an interplay of synopticscale circulations (i.e., the western disturbances and Indian summer monsoon) with HKH topography largely govern precipitation within the upper indus basin (UIB) that varies in space and time, and with elevation (e.g., Hewitt, 2005;Bolch et al., 2012). In contrast, the Lower Indus (LI) has an arid to semi-arid climate and depends heavily upon the meltwater from the UIB (Immerzeel et al., 2015).
Projected global warming and changes in large-scale circulations will impact the Indus hydrology by altering input precipitation and glacier contributions (e.g., De Souza et al., 2015). An increase in water demand due to rapidly growing population (UN, 2019) and future environmental conditions will also disrupt the hydrological balance at the basin-level (e.g., Lutz et al., 2016b). Therefore, a pragmatic assessment of the basin's climate response towards projected global warming is essential to support integrated water management.
Precipitation projections are fundamental to assess future water availability, cryosphere stability, and irrigation demand across the basin. Presently, such projections can only be derived by downscaling, either statistically or dynamically, the output of general circulation models (GCMs, e.g., Wilby et al., 2000). Both downscaling methods have been adopted within the UIB to estimate future precipitation (e.g., Akhtar et al., 2008;Lutz et al., 2016b;Khan and Koch, 2018). Nonetheless, uncertainty about projected precipitation (signal strength and direction, seasonality, spatial patterns, temporal evolution, and cryosphere stability) remains high (e.g., Gebre and Ludwig, 2014;Ali et al., 2015;Khan et al., 2015;Hasson, 2016;Su et al., 2016). Such uncertain scientific feedback can undermine the policy response to minimize the projected vulnerability of a large population.
Most uncertainty in downscaling studies stems from the choice of models (e.g., Heo et al., 2014), how precipitation is derived from these models (Pomee et al., 2020), and adopted observations (e.g., Palazzi et al., 2013). For example, the latest GCMs in the Coupled Model Intercomparison Project Phase 5-CMIP5 (Taylor et al., 2012) still show major limitations in representing critical dynamic and thermodynamic processes over this complex region (e.g., Sperber et al., 2013;Ashfaq et al., 2017). Regional Climate Models (RCMs) often provide improved simulations, but the evaluation of the CORDEX-SA experiments (e.g., Kulkarni et al., 2013;Mishra, 2015;Hasson et al., 2019) has shown only limited success over this region. Therefore, using an arbitrarily selected single model to simulate precipitation over topographically heterogeneous UIB (e.g., Akhtar et al., 2008;Mahmood and Babel, 2012;Khan et al., 2015) may lack fidelity for regional adaptations.
Model ensembles encompassing a broader uncertainty (e.g., Gleckler et al., 2008;Sperber et al., 2013) are more useful for such complex regions. Nevertheless, the selection of representative models is a challenge due to rapidly growing climate models. For example, the CMIP3 (Meehl et al., 2007) contains 25 GCMs, whereas more than 60 GCMs are available in the CMIP5 archive. Different model ranking metrics are available for ensemble selections. For instance, the past performance criterion identifies models demonstrating better skills in reproducing past climate (e.g., Christensen et al., 2010;Biemans et al., 2013). The so-called envelope approach only considers the range of projected changes in the variable(s) of interest during ensemble selections (e.g., Sorg et al., 2014;Warszawski et al., 2014). Sometimes both approaches are combined to select models that better simulate the observations and encompass broader future evolutions (e.g., Giorgi and Mearns, 2002;McSweeney et al., 2015). Issues like subjective decisions during the model rankings, the choice and effectiveness of performance metrics, and inter-model similarities can still induce uncertainty in these ensembles (e.g., Knutti et al., 2017). Therefore, a thoughtful model selection to better serve the intended objectives is still an ongoing research challenge.
Within the UIB, some studies used ensemble approaches by either using the past performance (e.g., Hasson, 2016) or a combination of future spread and past simulations (e.g., Immerzeel et al., 2013;Lutz et al., 2016a;Khan and Koch, 2018). However, using precipitation for model ranking and subsequent projections remains common, which may induce significant uncertainties since even the latest GCMs do not reliably simulate precipitation due to complex generating mechanisms and high spatiotemporal variability (Trigo and Palutikof, 2001;Mueller and Seneviratne, 2014). Such model deficiencies further manifest over high-mountain regions like the UIB, where orography causes additional precipitation variability. Without long-term, reliable, and consistent observations to account for regional orography of the UIB (e.g., Immerzeel et al., 2015;Pomee et al., 2020), even the bias-corrected model precipitation may lack reliability.
While additional data availability is still an ongoing issue, the methodological considerations can reduce some uncertainties. For instance, using large-scale atmospheric fields instead of precipitation to develop downscaling models offers one such promising alternative. The latest GCMs show better skills in simulating these patterns than raw precipitation output (e.g., Kaspar-Ott et al., 2019). We used such atmospheric variables for precipitation analysis over the entire Indus basin by focusing on the UIB. Firstly, we identify precipitation-governing variables from ERA-Interim reanalysis within a robust statistical downscaling framework by accounting for seasonality and spatial variability of the observed precipitation. For GCM ranking, we quantify the correspondence between ERA-Interim predictors and corresponding CMIP5-simulations during the historical period using Taylor diagrams (Taylor, 2001). During GCM selection, we also consider predictor influence on precipitation by weighing the GCM predictors with downscaling models' regression coefficients. The sub-regional approach further helps to evaluate model performance over multiple regions. Given the regional complexity, we also evaluate ERA-Interim predictors' robustness against two additional reanalysis datasets (NCEP-NCAR-II and ERA5) to quantify the reference uncertainty. Finally, we use governing predictors to derive ensemble precipitation changes over the study basin for selected radiative-forcing scenarios. We also evaluate the robustness of the change signals by using signal-to-noise-ratio (SNR). The extension of analysis further to the LI helps to analyse the supply-demand perspective for water management at the basin-scale.
In summary, we consider predictor (multiple reanalyses), predictand (homogeneous and long timeseries stations), and GCM (ensemble) levels to provide a more realistic precipitation assessment for the Indus basin during the 21st century. To our knowledge, such fine-scale analysis has never been performed in this region, and hence a different simulation perspective is provided. Our methodology can also facilitate climate assessment studies elsewhere.

| STUDY AREA
The present study focuses on the Indus basin of Pakistan, covering an area of 1.31 million km 2 (FAO, 2011), which ultimately drains into the Arabian Sea (Figure 1a). Although the Indus River system originating within highmountains of the HKH region also shares its UIB (4.03 × 10 5 km 2 , Dahri et al., 2018) with China, India, and Afghanistan (Figure 1b), we could not include these regions in our analysis due to constraints on updated data availability. The Indian summer monsoon, western disturbances, and the Tibetan anticyclone (e.g., Wake, 1989) largely govern the basin-wide precipitation that varies significantly over the UIB. Besides, moisture recycling within the Indian subcontinent and HKH region also plays a significant role in annual precipitation (e.g., Curio and Scherer, 2016). The LI (0.72 million km 2 ) depends heavily on seasonal water supplies from the UIB.

| DATA AND METHODOLOGY
This study builds on previous work (Pomee et al., 2020), where predictor-predictand relationships are used to model observed precipitation at sub-regional scales over the study basin. The present study concentrates on quantifying reference uncertainties, GCM ranking, and precipitation downscaling results under future climate change constraints using previously identified predictors and seasonal precipitation characterization of the basin. In the following sections, we provide a brief description of the data, the adapted regionalization scheme, and the model development process to give the necessary background. More detailed information can be found in Pomee et al. (2020).

| Precipitation: data and regionalization
We use monthly precipitation time series of 58 meteorological stations located across the study basin ( Figure 1b). These stations provide historical (35 stations, 1979-2015) and relatively short-term data (23 stations, 1994-2015) over the lower-elevation and high-elevation (HA) regions of the UIB, respectively. The HA stations (average data length of 17 years) have considerably expanded spatial-altitudinal coverage within the UIB to provide valuable insights about regional orography that governs the basin hydrology. Most study stations are located in the UIB to account for topographic complexity and its significance for Indus flows. More information about the study stations is available in Table S1. F I G U R E 1 (a) The transboundary Indus River basin and demarcation of the study area, where the shading differentiates between the Upper Indus Basin and Lower Indus regions. (b) The locations of the meteorological stations (numbered) used in our study. The circles represent the long time series , and triangles show the recent highaltitude stations with shorter series . Note that, the shading scheme in (b) represents the altitudinal variations (elevations above mean sea level) within the study basin. Table S1 of the Supporting Information provides more information about these stations. Refer to the online version for information on the shading schemes adopted in the figure [Colour figure can be viewed at wileyonlinelibrary.com] Based on observations, we identify three major precipitation seasons for further analysis. These include the winter season-WS (DJFM), pre-monsoon-PMS (AMJ), and monsoon season-MS (JAS), respectively. We group the time series of study stations into these seasons and check data for completeness (Moberg et al., 2006). We use four different statistical procedures to test homogeneity of the time series after Wijngaard et al. (2003).
We employ K-means cluster analysis on all 58 stations to identify precipitation regions with similar covariance using Spearman correlation as a distance measure (Wilks, 2006). We set the objective function during clustering to maximize (minimize) correlation within (across) the regions to define sharp regional boundaries. Considering more uncertainty over the HA, we perform another regionalization experiment covering only elevated parts of the UIB. Our regionalization schemes identify, on average, four precipitation clusters within the UIB, accounting for the high spatial variability. These regions represent precipitation dynamics in the southern Himalayans, trans-Himalayans (including the Karakoram), and northwestern parts of the UIB during each season. Besides, on average, two regions describe LI precipitation variability. Regional representative (RR) stations for each precipitation cluster are selected through multiple considerations (including homogeneity) and serve as predictands. The seasonal precipitation regionalizations for both experiments are shown in Figures 2 and 3.

| Predictors: data and principal component analysis (PCA)
We use variables from ERA-Interim reanalysis (Dee et al., 2011) at different pressure levels to identify major dynamic and thermodynamic drivers of regional precipitation. We consider a larger domain for circulationdynamic variables (10-100 E, 10-60 N, spatial resolution 2 × 2 ) compared to the domain for thermodynamic fields (64-80 E, 22-40 N) to account for large-scale dynamical and more localized thermodynamic influences on precipitation. We did not use ERA5, as its data was not publically available at the beginning of our research during early 2016.
We perform S-mode Varimax-rotated Principal Component Analysis (S-mode PCA) on each predictor field for dimension reduction (e.g., Preisendorfer, 1988). Following the modified dominance criterion (Philipp, 2003) with some additional constraints (Pomee et al., 2020), we retain a maximum of 20 PCs to explain the predictor variance adequately. The resulting PC scores serve as predictor time series, and corresponding PC loadings define the centres of predictor variations.

| Generalized linear models and statistical downscaling framework
We adopt a generalized linear model (GLM) framework (Mc Cullagh and Nelder, 1989) to model predictor-predictand relationships within a robust crossvalidation setting. For those predictand cases containing exact zeros in their time series, Tweedie exponential dispersion models (e.g., Dunn, 2004) and otherwise GLM gamma models are used. We use mean squared error skills scores (MSESS) as a performance criterion (Wilks, 2006) and consider multicollinearity among predictors to identify predictor combinations that best resolve the observed precipitation. The information about identified governing predictors and the regression models' statistical performance is provided in Tables 1, S2, and S3, respectively. Appendix S1 provides details of the downscaling framework after Pomee et al. (2020).

| GCMs: data and precipitation downscaling
We first consider all CMIP5 GCMs (Taylor et al., 2012). However, the availability of governing predictors (Table 1) in the historical period restricts this number to 29. Due to high mountains in our study domain, only eight of those 29 GCMs could provide complete spatial coverage of the required predictors. Many GCMs do not provide spatially complete lower-tropospheric predictors over the high mountains due to the intersection of pressure coordinates with mountain elevations. These spatial inconsistencies may emerge as many modelling centres avoid employing any vertical interpolation (or extrapolation) algorithms over the mountain regions when transforming data from the model to pressure levels. These data gaps are too large to be effectively filled with an interpolation scheme at the end-user level and restrict the computation of spatially consistent historical and future predictor PCAs. Therefore, we could only analyse the output of eight models during the historical  and two future time-periods, covering the mid (2041-2070) and end-of-the 21st century (2071-2100) for precipitation downscaling. We could not consider model independence (e.g.,  in detail due to the smaller ensemble but investigated its influence over a sample region in Section 4.3. We consider RCP4.5 and RCP8.5 scenarios for precipitation projections. RCP4.5 portrays a future where technological advancements will help stabilizing greenhousegas emissions after 2,100. Conversely, RCP8.5 depicts a populous world without abatement efforts where F I G U R E 2 Precipitation regionalization of the entire Indus Basin of Pakistan using K-means cluster analysis on seasonal scales. Different colours represent the identified regions. The circles with the same colour show similar precipitation variability and thus belong to one precipitation region. The numbered circles indicate the location of regional representative (RR) stations of the respective precipitation regions. Triangles represent those stations that could not be clustered to any of the identified precipitation regions. The  radiative forcing reaches the maximum of 8.5 watts/m 2 in 2100 (Van Vuuren et al., 2011). These RCPs cover a plausible range of radiative forcing (Sanford et al., 2014) to support adaptations and comparison with previous studies. We only consider the first realization ('r1i1p1') of these GCMs during the historical and future periods. Table 2 provides information on the CMIP5 models used in our study.
Before downscaling, the GCM data is conservatively re-gridded (2 × 2 ) to match the adopted ERA-Interim resolution. The modelled predictors are standardized over the historical and future time-slices (separately for each scenario) and then projected onto the corresponding loading patterns of the ERA-Interim variables to generate new predictor time series (more details on this method are available in Von Storch and Zwiers, 1999). These new predictor time series are used in the regression models (-Tables S2 and S3) to derive downscaled historical and projected precipitation totals. The difference between downscaled precipitation during the historical and two future time-slices (separately for each RCP and time slice) is used to compute median precipitation changes across the basin.
We further evaluate the robustness of projected change signals in light of observational uncertainty by computing a signal-to-noise ratio (SNR). The ratio uses median precipitation changes (signal) simulated by the individual models and their ensembles and corresponding standard deviations of the historical period (noise). SNR >1 indicates that the change signal exceeds the internal climate variability.

| GCM ranking process
A stepwise procedure is used to rank the GCMs according to their ability to simulate precipitation-governing variables of the ERA-Interim reanalysis during the historical period. ERA-Interim reanalysis provides simulations from 1979 onwards; hence, we could only compare the reference-model predictors during the overlapping historical period .
(I) Initially, S-mode PCA is performed (Section 3.2) for every governing predictor (Table 1) of each individual GCMs to extract the same number of PCs as from ERA-Interim.
(II) Subsequently, the model PC loadings are compared with corresponding ERA-Interim loadings (separately for each GCM) using Taylor diagrams. A simple performance score (PS) derived using two of the three summary statistics of Taylor diagrams is computed to quantify the predictor correspondence. Mathematically, the PS is: where PS = performance score. For a perfect predictor agreement, PS = 1. CR = pattern correlation between the reference (ERA-Interim) and model (GCM) loadings. For a perfect phase match, CR = 1. NSD = normalized ratio of variance (standard deviation of the reference and model loadings). Ideally, the NSD should also take the value 1. Note: Column 1 provides the list of large-scale predictors used in this study where zg, va, ua, hur, hus, and psl denote geopotential heights, meridional wind, zonal wind, relative humidity, specific humidity and, mean sea level pressure fields, respectively. The number after each predictor reflects the atmospheric level (pressure level in hPa). Columns 2-4 represent the seasonal frequencies of different predictors chosen to resolve observed precipitation dynamics over the study basin. The last column shows the average predictor frequencies over different seasons. Under ideal conditions, the PS will attain its highest unit value due to maximization of the phase correspondence (i.e., CR = 1) and the same magnitude of predictor spread (i.e., the term NSD − 1 becomes zero) between the reference and model simulations. Similarly, a smaller PS value will show a weaker predictor correspondence. The PS magnitude will also intuitively influence the third summary statistics (i.e., standardized RMSE), where its maximum value (PS = 1) will ensure zero error. Conversely, the smaller values (PS <1) will reflect higher errors, though not following a clear linear trend due to the typical relationship among these three summary statistics (see Taylor, 2001). Thus, the PS contains useful information about the strength of correspondence between the reference and model-simulated fields and can be used to identify the best-matching pairs for every governing predictor.
(III) We draw two separate sets of Taylor diagrams for each precipitation region and season. In the first set of diagrams, a given reference PC is compared with all modelled PCs (separately for each GCM) to compute the corresponding PS. The reference-model pair demonstrating the maximum PS among all PCs of a GCM is selected as the best GCM-PC for that particular reference. Thus, each GCM has one best corresponding PC for a given reference.
This process is repeated for all other PCs and predictors of the regression models. Subsequently, all best-matching (individual) PCs of different predictors are grouped into the second set of Taylor diagrams (separately for each GCM) to assess the ability of individual GCMs in representing ERA-Interim precipitation-governing predictors over a region. The summary statistics of the second Taylor diagrams are used to compute the average unweighted PS for each GCM. The PS is termed unweighted due to equal PC weighting in the computation.
(IV) Given that each PC has a different influence in a regression model, we further adopted (absolute) regression coefficients of the PCs as weights and computed the weighted PS. Thus, a model with the highest (lowest) weighted PS score can be identified as the best (worst) GCM due to its improved (poor) simulations for more important predictors.
This process (Step I-IV) is repeated for all sub-regions to identify the best regional GCM in different seasons.
(V) Finally, we consider GCM performance over multiple regions to identify models that show superior simulations over the whole spatial scales of the UIB and LI, respectively. We prefer a GCM demonstrating better simulations over multiple sub-regions. Spatial consideration is important since an outlier may strongly influence a model's overall PS (e.g., very high PS just over one subregion).
3.5.1 | GCM ranking: an example We demonstrate the above procedure over a sample WS region (R6). Six PCs from two different predictors (hur1000 and zg700) were selected to skillfully resolve observed precipitation over this region (Table S2). We represent these as reanalysis predictors 1 (R1) and 2 (R2), respectively, and construct three different sets of Taylor diagrams to illustrate each selection step. Firstly, the Taylor diagrams (Figure 4) visually outline the best-matching pairs of the reanalysis PC1 of the first predictor, that is, hur1000 ('R1.1'), and modelled PCs of this predictor for every GCM. We use these Taylor diagrams' statistics to identify the best PC match (for every GCM) using the PS and plot them through the second set of Taylor diagrams ( Figure 5). We repeat this process for the remaining five PCs to identify bestmatching pairs for each GCM (not shown). Subsequently, we group all six best-matching pairs of each GCM into a third set of Taylor diagrams ( Figure 6). Finally, we compute the weighted PS (from the statistics of Figure 6) using regression coefficients (absolute) as weights to identify the best-performing model. In this example, CMCC-CM appears as the best model due to its highest average PS of 0.78 (see Table 3).

| Quantification of reference uncertainty
We further consider governing predictors (Table 1) from the latest ERA5 (Hersbach et al., 2020) and from NCEP-NCAR-II (Kalnay et al., 1996) reanalysis datasets to evaluate the usefulness of ERA-Interim for precipitation modelling. After re-griding to a common spatial resolution (2 × 2 ), the predictors from these additional reanalysis datasets were subjected to a similar S-mode PCA. The resulting PC loadings are compared (separately for each reanalysis) with ERA-Interim using the procedure described in Section 3.5 (without Step V). Thus, the weighted PS of these other two reanalyses can define the range and average magnitude of reference uncertainty. Table 3 quantifies the strength of ERA-Interim predictors' correspondence with the other two reanalysis datasets (ERA5 and NCEP-NCAR-II) and each available GCM in different precipitation seasons. The regions are grouped into the UIB and LI scales to estimate ERA-Interim predictors' spatial effectiveness (uncertainty) against each dataset by computing average PS over these two regions.

| Reference uncertainty
First, we use the output of Table 3 to quantify reference uncertainty. Generally, the three reanalysis datasets indicate more robust simulations of the various dynamic drivers compared to the thermodynamic drivers of regional precipitation (Table 1). The reliability of ERA-Interim predictors varies with season, reanalysis, and over the regions. For example, ERA-Interim predictors show more robust correspondence with ERA5 (NCEP-F I G U R E 4 Identification of the best-matching pairs of the reference (ERA-Interim) and model (GCMs) simulated PCs of a given predictor using Taylor diagrams. Each Taylor diagram elaborates the PC matching process separately for each GCM by comparing all modelled PCs of a predictor field with a given reference PC. In this example, we use reanalysis PC1 of hur1000 to identify the best matching PC of hur1000 as simulated by different GCMs. The letter R (M) and associated numbers represent the reference (model) related PC information for a given predictor field in the legend key. For instance, in the text 'R1.1', the letter 'R1' stands for the first reference predictor (i.e., hur1000), and the following number ('.1') reveals the identification of its PC (i.e., PC1 of hur1000). Similarly, in the text 'M1.1', the letter 'M1' and the following number ('.1') show the first GCM (CMCC_CM) and the identification of its simulated PC (i.e., PC1) of the predictor field (i.e., hur1000), respectively. Likewise, The letter 'M2.6' shows the second GCM (CMCC_CMS) and its simulated PC (i.e., PC6). Thus, identifying the modelled PC that best corresponds to PC1 of the reanalysis predictor (hur1000) can be seen in each of the above Taylor diagrams. In the case of GCM4 (CanESM2), among all six modelled PCs, its first PC (M4.1) better corresponded to the given reanalysis PC (R1.1) [Colour figure can be viewed at wileyonlinelibrary.com] NCAR-II) over five different UIB regions during the MS, as indicated by a higher average PS of 0.84 (0.72). Considering the mountainous terrain and MS complexity, represented by a set of diverse predictors located across the troposphere (see Table S1), the strong predictor correspondence justifies ERA-Interim use for precipitation modelling over the UIB.
Similar predictor robustness is noticed over two LI regions (ERA5 = 0.81 and NCEP-NCAR-II = 0.70) and demonstrates the ERA-Interim predictors' reliability for F I G U R E 5 Effectiveness of the PS in identifying the best reference-model PC combinations of Figure 4. The other details are the same as in Figure 4. The statistics of Figure 4 are used to compute the PS for all modelled PCs (separately for each GCM). The model simulated PC that shows the maximum PS is identified and plotted through separate Taylor diagrams of Figure 5. As the best PCs in Figures 4 and 5 are precisely the same; therefore, the PS can be used to identify such pairs numerically [Colour figure can be viewed at wileyonlinelibrary.com] the basin-wide MS precipitation analysis. ERA-Interim also shows strong predictor similarities with ERA5 and NCEP-NCAR-II during the WS (0.83 and 0.74) and PMS (0.76 and 0.71) over multiple UIB regions. A similar predictor correspondence over different LI regions during these two seasons further highlights the effectiveness of ERA-Interim for westerly-dominated regimes.
Although governing predictors are robust among different reanalyses datasets, still some interesting patterns emerge. For example, free atmospheric predictors that largely dominate the seasonal precipitation models (e.g., MS-R3, PMS-R1, and WS-R1 in Table S2) show more similarities among reanalysis datasets, partly because the atmospheric boundary layer exerts lesser F I G U R E 6 Final Taylor diagrams representing all six modelled PCs (separately for each GCM) of the two governing predictors (i.e., R1 = hur1000 and R2 = zg700) that show maximum correspondence with ERA-interim simulated PCs. The statistics of these Taylor diagrams are used to compute the PS, and the model that shows the maximum weighted PS (across all six PCs) can be identified as the best model. In this example, GCM2 (CMCC-CM) shows the highest PS (see Table 3) and hence selected as the best GCM for this particular WS region. Visual inspection of these Taylor diagrams also confirms this selection, where its modelled PCs show higher correspondence with reanalysis data than other GCMs [Colour figure can be viewed at wileyonlinelibrary.com] T A B L E 3 Quantification of the reference and model predictor correspondence on seasonal scales across the Indus Basin of Pakistan using weighted PS [Colour influence over these predictors, and most of the predictor PCs are located outside the mountain region (not shown). Therefore, their simulations are more robust across different numerical models. Conversely, the simulations of near-surface predictors (notably hur1000) show larger differences (e.g., PMS-R6 and WS-R6), which may arise due to the stronger influence of the atmospheric boundary layer. Note that, all thermodynamic PCs are located over and in the surrounding of high mountains due to a reduced predictor domain for these variables. Differences in the numerical models and interpolation issues for the near-surface variables over the mountains (e.g., Palazzi et al., 2013) may further add to such simulation disparities. The contribution of hur1000 in precipitation models is low, and the weaker correspondence is mainly limited to some LI regions, particularly during the PMS and WS. Such region-specific discrepancies may not substantially influence the overall regional suitability of ERA-Interim, particularly over the UIB that controls river flows during all seasons. Another pattern relates to the seasonal outperformance (higher PS) of ERA5 over NCEP-NCAR-II (except the WS-LI regions due to more differences in simulations of hur1000). That seems logical as both ERA-Interim and ERA5 are ECMWF reanalysis and share more similarities in their simulation schemes. However, we purposefully select NCEP-NCAR-II that shares the same temporal resolution with ERA-Interim (from 1979 onward), maximizes the use of post-70s satellite measurements, and provides the opportunity to assess ERA-Interim variables more independently.
The difference between a perfect predictor match (PS = 1) and actual correspondence among reanalysis data can define the reference uncertainty range. ERA5 mostly defines the lower, and NCEP-NCAR-II outlines the upper bound of such uncertainty. Table 3 provides reference uncertainty during the MS (16-28%), PMS (24-29%), and WS (17-26%) over the UIB. The corresponding LI uncertainty ranges from 19% to 30%, 24% to 25%, and 18% to 27%, respectively. Considering regional heterogeneity, such strong predictor correspondence (low uncertainty) of ERA-Interim with different reanalysis datasets fully justifies its use for constructing downscaling models.

| GCM ranking
Similarly, we identify better-performing GCMs using weighted PS. During the MS, all GCMs show relatively low correspondence with ERA-Interim predictors over different UIB regions (Table 3). Such poor predictor similarities (except for R3) highlight the MS complexity and its highly uncertain representation in most GCMs (e.g., Ashfaq et al., 2017). Among available models, CMCC-CMS turns up as the best single model (PS = 0.50 over five UIB regions). This model outperforms over two larger sub-regions (R1 and R3; see Figure 3c for regional identification) and appears the second-best model for the northwestern region (R4). Moreover, its performance over the remaining two UIB regions (R5 and R7) is comparable with other GCMs. Therefore, its use for MS projections over the entire UIB seems more justified among available models. Although this model only shows 50% predictor correspondence, a high reference uncertainty (up to 28%) also needs to be considered for judging its true fidelity. In contrast, Can-ESM2 and MPI-ESM-LR show the least predictor correspondence over the UIB. The MS predictors for two LI regions are better simulated by MPI-ESM-LR (PS = 0.52).
All GCMs showed relatively high predictor correspondence with ERA-Interim during the PMS. However, CNRM-CM5 offers more consistent performance over four UIB regions (PS = 0.61). In contrast, MPI-ESM-LR better represents the LI predictors. Both models provide improved simulations of hur1000, a primary predictor to simulate precipitation over R5 and R6 regions (Tables S2  and S3). Consequently, the basin-wide compatibility of these models with ERA-Interim predictors improves significantly.
During the WS, most ensemble members demonstrate the highest PS. Among these, MPI-ESM-LR shows maximum predictor correspondence over the UIB (PS = 0.67) by outperforming predictor simulations over the two large and cryosphere-dominated regions in trans-Himalayans (R1 and R5). It also shows comparable performance with other GCMs over the third UIB region (R3). As before, a different model (CMCC-CM) provides better simulations of the governing predictors over two LI regions.
In summary, no single model can effectively simulate precipitation dynamics at the basin-level without making significant compromises. However, our model ranking process can identify models demonstrating spatially consistent performance over the UIB (LI) with tolerable uncertainties. Significantly improved representation of the westerly circulations in our ensemble may increase our understanding and confidence about projected cryosphere dynamics that largely influence basin sustainability.
It is important to note that our model rankings only rely on predictor correspondence during the historical period and require a stationarity assumption for their future validity (e.g., Lanzante et al., 2018). The stationarity considerations may induce some additional uncertainties, as the predictor-predictand relationships and the models' ability for their simulations may alter in the future (e.g., Hertig et al., 2017). Stationarity violation may particularly happen in highmountain regions, where numerous feedback mechanisms exert influence on regional climate, and their future evolution may substantially differ from the observations.

| Diversity of the GCM ensemble
Only eight GCMs in our study may not adequately represent the entire CMIP5 simulation diversity. Using literature review and exploratory quantitative analysis, we evaluate our ensemble's regional efficacy. For instance, Khan and Koch (2018) assess the entire CMIP5 dataset under RCP4.5 and RCP8.5 scenarios to shortlist the GCMs defining a so-called full spectrum of future climate over the UIB. Their final models explaining Warm-Wet (Can-ESM2), Dry-Cold (MPI-ESM-LR), and mean climate (Nor-ESM1-ME) are available in our ensemble. Besides, our ensemble also includes their secondary considerations for Dry-Warm (CMCC-CMS) and Wet-Cold (CNRM-CM5) regional future (see Table 2). Similarly, the 14-model CMIP5 ensemble of Ali et al. (2020) also contains almost all our models. Such ensemble similarities justify the regional relevance of our models, despite its smaller sample size.
We further performed an additional analysis over a sample MS-UIB region (R3) whose precipitation predictor (va200) was available in most CMIP5 GCMs. We selected 15 GCMs belonging to 13 different modelling centres (Table 4) to cover most of the remaining institutional spread in the CMIP5. The historical uncertainties of these additional GCMs were similarly computed (see Section 3.5) and compared with the 8-model ensemble using the box and whisker plots (Figure 7). The ensemble comparison shows that nearly all 15 GCMs lie within the uncertainty spectrum defined by the 8-models; the medians are similar, and the best GCM (lower outlier of the 8-model ensemble) remains the same. The reduced uncertainty among the new GCM ensemble suggests close similarities between these models for MS simulations instead of providing a significantly different perspective. Therefore, both literature review and sample analysis demonstrate the regional suitability of our 8-model ensemble.
We also investigate the influence of model uniqueness  over the sample MS region. We computed historical uncertainties using 23 GCMs (model democracy) and compare them with 18 GCMs belonging to different modelling centres (model independence; see Tables 1 and 4). Using the box and whisker plots (not shown), we could not find significant differences in the two distributions. Besides, similar uncertainty estimates of the T A B L E 4 Details of the additional GCMs used to evaluate the CMIP5 representativeness of the 8-model ensemble ( independent models and the 8-model ensemble over the sample region also induce confidence in our ensemble simulations.

| Downscaled future precipitation changes
We use available GCMs to compute the median precipitation changes for two future periods (2041-2071 and 2071-2100) relative to the historical period  under RCP4.5 and RCP8.5. Figure 8 shows sub-regional multi-model ensemble (MME) and individual GCM simulated precipitation changes during the end of the 21st century (2071-2100) under both RCPs. The inter-model spread defines the range of uncertainty around MME signals. The corresponding changes during 2041-2071 are quantitatively different but show similar spatial patterns (not shown).

| WS projections
The WS projections under both RCPs (Figure 8a,b) show considerable positive changes (MME) over large parts of the UIB, covering mainly the HA regions of the trans-Himalayan and northwestern Hindukush (R5, R1; refer to Figure 2a for regional location). The positive signals are robust (all models project positive changes) and much stronger over the central Karakoram (R1) than the northwestern and eastern regions (R5) under both scenarios.
Some precipitation decrease appears along the lower elevations of the northwestern and southern Himalayans (R3) under RCP4.5 (Figure 8a), but it will eventually stabilize under RCP8.5 forcing. The precipitation changes suggest that increased warming (RCP8.5) will promote an elevation-dependent response in the UIB, where HA regions will receive more precipitation. The best seasonal model (MPI-ESM-LR) demonstrating high predictor correspondence during the historical period further supports these ensemble signals, particularly for RCP8.5. The likelihood of decreasing precipitation over the UIB remains low, as mainly negative signals are projected by two Norwegian models that demonstrate the lowest historical performance (Table 3). The large inter-model spread reflects a highly uncertain future over the UIB. The westerlies approaching the UIB bifurcate into northern and southern branches along the Himalayans (e.g., Pang et al., 2014). Increased strength and more northward penetrating westerlies at the end of the 21st century may increase dynamic forcing over the UIB to explain the typical spatial changes under warming scenarios. Note that circulation-dynamic predictors (Table S1) mainly govern the regional precipitation. The strong and northward located westerlies may continue to support the regional cryosphere and anomalous behaviour of the Karakoram.
Conversely, the precipitation changes are subtler over the two LI regions (R6 and R4) under both RCPs. For instance, the MME change over the upper irrigated plains (R6) shows some robust decrease (many models show negative signals), which confirms the weakening of the southern westerly limb. Meanwhile, the lower irrigated plains (R4) indicate strongly positive and robust changes under both RCPs. Local feedback mechanisms may strongly influence the future precipitation changes over this region. Overall, some decrease in LI precipitation seems plausible (though with more uncertainty) due to quantitatively more negative changes over a large region (R6) than the positive signals over R4.
F I G U R E 7 Comparison of the two GCM ensembles for simulating precipitation predictors (seven PCs of va200) over a sample MS region (R3) during the overlapping historical period . The uncertainties along the y-axis reflect the degree of mismatch between the ERA-Interim reanalysis's predictor and corresponding simulations of every available GCM (black dots) in the two ensembles. The GCM uncertainty is computed using the mathematical relationship (1 − PS)*100 [Colour figure can be viewed at wileyonlinelibrary.com] 4.4.2 | PMS projections PMS changes (MME) show a decrease in precipitation over all four UIB regions (R1, R5, R3, R7; see also Figure 3b) that further intensify with increased warming (Figure 8c,d). These changes exhibit spatial variability, and strongly negative signals (up to 25%) are noticeable in the lower northwestern regions (R5 and R7). Most individual models also project negative changes for these regions except Can-ESM2, but it shows the lowest predictor correspondence in the historical period. However, a large HA part of the trans-Himalayans region (R1) indicates a lesser decrease and more uncertainty (some models even project similar positive signals) under F I G U R E 8 Downscaled unweighted seasonal precipitation changes in % over the Indus Basin of Pakistan during the future period (2071-2100) relative to the historical period  under RCP4.5 and RCP8.5. The y-axis shows the identified precipitation regions and the corresponding range of regional altitudes that decreases from top to bottom. The regional altitudes are expressed as the elevations above mean sea level (m-amsl). The WS, PMS, and MS precipitation changes under RCP4.5 and RCP8.5scenarioes are shown by the subplots (a,b), (c,d), and (e,f), respectively. The coloured circles (triangles) show the individual GCM (MME) simulated precipitation changes. The solid blue line indicates no change compared to the historical period. Note that, the range of the x-axis (i.e., precipitation changes) is different in these panels [Colour figure can be viewed at wileyonlinelibrary.com] both scenarios. Moreover, a wetter northeastern region (R3) shows nearly neutral changes, and the best seasonal model (CNRM-CM5) instead projects a slight increase in regional precipitation. Generally, a precipitation decrease is more robust over the UIB, though its intensity and reliability reduce with elevation and over eastern parts. The decreasing trend is also visible over large parts of the southwestern highlands and the LI plains (R6). In contrast, the upper irrigated plains and northeastern rainfed regions (R4) show some precipitation increase that becomes stronger under RCP8.5. Generally, a smaller inter-model spread suggests more robustness of the sub-regional signals during the PMS.
Weaker and more northward oriented westerlies may reduce the moisture advection to reduce future precipitation, particularly over the lower northwestern regions, where relative humidity has a strong impact (Table S2). Alternatively, the projected seasonal warming (e.g., Ashfaq et al., 2020) may reduce the saturated atmospheric moisture content (i.e., relative humidity) to delay precipitation onset.

| MS projections
We assess an overall precipitation increase over the entire Indus basin under both RCPs (Figure 8e,f). Most individual models, including the best seasonal models, further support these sub-regional trends that intensify further under RCP8.5. However, the change signals significantly vary over different UIB regions (R7, R6, R1, R3, R5; see Figure 3c). For example, the highest elevations towards the Karakoram-West (R7) show some consistent decrease (MME = approx. −8%). In contrast, a northwestern region (R4) and a larger area around the central Karakoram (R1) depict an up to 12% increase in future precipitation. Similarly, a large and the wettest part along the southern Himalayans (R3) also shows increased precipitation (MME = 4%) in both scenarios. Another northwestern region representing the lower elevation of the Hindukush (R5) also indicates some positive signals (up to 2%). These distinct precipitation changes (predominantly positive) support the monsoon system's strengthening and further penetration into the Northwestern and trans-Himalayans regions under increased warming scenarios. Such intense MS circulations may continue to support the regional cryosphere and downstream water needs in the future. More uncertainty over the magnitude and the direction of precipitation signals highlights the complex interplay between the MS currents and UIB topography.
Similarly, two LI regions (R6 and R2), covering spate irrigation in the southwestern mountains, irrigated plains, and coastal areas, mainly project positive changes under both RCPs. Therefore, increased water development potential in spate-regions and a slight decrease in net irrigation over the plains are likely under future warming. Compared to the UIB, the inter-model agreement over the LI precipitation changes is high.
4.5 | Impact of model weighting on ensemble changes Figure 9 shows the impact of model weighting (Table 3) on sub-regional ensemble changes during 2071-2100 under both RCPs. Generally, better-performing models (models with higher weights) have a relatively small impact on change signals due to the adopted weighting scheme (maximum value of 1) and intermodel similarities in our ensemble. Still, the weighting has demonstrated refinement of the change signals by altering magnitude in many cases. For instance, the model weighting shows maximum influence during the MS season ( Figure 9c) by mainly strengthening positive changes over the UIB under RCP8.5. A considerable increase in precipitation over the central Karakoram and northwestern HA, including the wettest part (R3) in the UIB, is apparent. This pattern also continues into the southern highlands (R6) and irrigated plains (R2) in the LI. The strengthening of MS circulations appears more in betterperforming models under the intense warming (RCP8.5) than RCP4.5 forcing, where a mixed response mostly prevails.
In contrast, the model weighting increases PMS dryness over most of the UIB for both scenarios (Figure 9b). Thus, the weakening of the westerly system during PMS is more consistent among better-performing GCMs, particularly under RCP8.5.
During the WS, the model weighting has mainly enhanced the positive signals over the UIB (Figure 9a) to indicate a slight strengthening of the westerly circulations in better performing models. Since HA regions receive much higher precipitation (e.g., Immerzeel et al., 2015;Pomee et al., 2020), even a slight increase in precipitation signals during the main seasons (MS and WS) can substantially increase the actual precipitation over the UIB under RCP8.5.
Using model weights in an ensemble setting can explain additional insights about projected changes in complex regions to assess confidence in projected signals despite the lesser impact.

| Robustness of the change signals
We use SNR to evaluate the strength of projected signals over the observational uncertainty, which remains very high over this complex region. Figure 10 shows the WS distribution of SNR as a function of regional elevations during 2071-2100 under RCP scenarios. Generally, the sub-regional ratios are higher under RCP8.5 forcing than under the RCP4.5 scenario, confirming the positive role of increased warming on regional precipitation changes. These ratios also demonstrate elevation dependency of future changes by depicting positive (strong) signals over HA regions compared to negative (weak) values over low elevations. The spatial patterns of SNR further support future strengthening and the northward oriented westerly regimes that will induce distinct positive changes over the HA-UIB.
The other seasons also show similar patterns for both RCPs (not shown).

| Downscaled precipitation: HA-UIB regionalization scenario
We also separately compute precipitation changes, model weighting, and SNR over HA regions of the UIB. Figure 11 only shows the precipitation changes during 2071-2100 under RCP8.5. These seasonal changes further confirm our previous findings (Section 4.7) about the strengthing and northward penetrating westerlies during F I G U R E 9 Impact of the model weighting on ensemble precipitation changes (i.e., MME) during 2071-2100 under RCP4.5 and RCP8.5 scenarios. (a) The WS, (b) the PMS, and (c) the MS sub-regional changes under both RCPs. In the legend, the 'unw' stands for unweighted, and 'w' represents the weighted ensemble changes for each RCP simulation [Colour figure can be viewed at wileyonlinelibrary.com] the WS and easterlies during the MS to promote positive precipitation changes over the high-elevation UIB. The PMS dryness over large parts of the UIB is also verified, indicating a future weakening of the westerly circulations during this season. The magnitude of changes differs, but the spatial patterns are similar for RCP4.5 (not shown).

| Further discussions
Our main findings support the results of some earlier studies. For example, Ali et al. (2015) and Khan and Koch (2018) also reported elevation-dependent precipitation changes over the UIB. Similarly, Lutz et al. (2016b) concluded a drying (increased precipitation) during the PMS (MS and WS) months. Although their changes were different in details, the spatial patterns of increasing precipitation from west-to-east confirm our findings. Bokhari et al. (2018) also estimated a reduction in future precipitation over the Kabul river basin, located in the northwestern UIB. Studies of Archer and Fowler (2004), Khattak et al. (2011), andForsythe et al. (2014) also projected increased precipitation during the main seasons.
However, our findings contrast to Palazzi et al. (2013Palazzi et al. ( , 2014, who found insignificant WS changes over the UIB. Using direct GCM precipitation, which contains significant wet biases over the region, in these studies may reduce signal strength. A bias-corrected analysis, despite some criticism (e.g., Ehret et al., 2012;Mezghani et al., 2017), may sharpen or even reverse the signal direction (e.g., Hagemann et al., 2011;Navarro-Racines et al., 2020).
Similarly, the intensity of drying (Khan and Koch, 2018) and a remarkable increase of projected precipitation over the UIB, particularly during PMS months, as shown by Hasson (2016) under RCP8.5, are not in agreement with our analysis. Surprisingly, both studies used the output of CORDEX-SA experiments but concluded a contrasting regional outlook, mainly due to differences in the reference datasets. For example, Hasson (2016) only used three short-term HA stations located along the Karakoram for scaling, which might induce a wet bias. In contrast, the adopted climatology of Khan and Koch (2018) could produce a dry tendency, as some recent studies (e.g., Immerzeel et al., 2015;Dahri et al., 2018) showed much higher precipitation over the UIB compared to the precipitation amounts used in their study.
Using the so-called model democracy (e.g., Knutti, 2010) and appreciating the complexity and future uncertainty of F I G U R E 1 0 WS distribution of the SNR to evaluate the robustness of projected precipitation changes under RCP4.5 (a) and RCP8.5 (b) scenarios over Pakistan's Indus Basin. The ratios are computed by dividing median changes (signal) during 2071-2100 with standard deviations (noise) of the historical period . The absolute value of the SNR along the Y-axis indicates the magnitude of signal strength. The blue horizontal line serves as a reference and indicates no strength of the projected signal. The precipitation regions are shown along the x-axis and arranged in decreasing altitudes from left to the right [Colour figure can be viewed at wileyonlinelibrary.com] the climate system, some studies criticize model rankings (e.g., Chiew et al., 2009;Hasson et al., 2019). We argue that some known model performance may reduce future climate uncertainties.
We also analysed future precipitation changes over the sample MS region (R3) during 2071-2100 under selected RCPs using the 8-and 14 model ensembles ( Table 4). The two ensembles showed strong similarities (not shown). Outliers in the 14-model ensemble were mostly those GCMs that demonstrated the lowest historical performance (e.g., IPSL-CM5A-MR, GFDL-ESM2, FGOALS-g2, and INMCM4). These models are known for their poor MS simulations over the study region (e.g., McSweeney et al., 2015) and could have been eliminated during the model ranking process. In summary, the 8-model ensemble adequately captures the future simulated by a larger CMIP ensemble and stands useful for precipitation analysis over the region for practical purposes.
Model-uniqueness influence on future precipitation signals over the sample region also remains negligible (not shown).

| SUMMARY AND CONCLUSIONS
We assessed future precipitation changes over Pakistan's Indus basin that primarily derives runoff from cryospheredominated watersheds in the UIB. Large-scale atmospheric patterns were used as predictors for constructing downscaling models, uncertainty quantification, GCM selection, and F I G U R E 1 1 Downscaled unweighted precipitation changes from 1976-2005 to 2071-2100, using the HA-UIB regionalization experiment's output under RCP8.5. Subplots (a), (b), and (c) represent sub-regional precipitation changes over the UIB during the WS, PMS, and MS, respectively [Colour figure can be viewed at wileyonlinelibrary.com] subsequent projections, as raw precipitation of the GCMs still lacks reliability. Such GCM limitations manifest over high-mountain regions like the UIB, where complex processes govern precipitation variation at sub-grid levels.
Firstly, we identify precipitation-governing predictors from ERA-Interim reanalysis within a robust statistical downscaling by accounting for spatial variability of the observed precipitation on seasonal scales. K-means cluster analysis was used for basin characterization. We further evaluated ERA-Interim predictors' robustness against two other reanalysis datasets (ERA5 and NCEP-NCAR-II) to demonstrate their usefulness for precipitation modelling.
We also identified better-performing GCMs by comparing model-simulated predictors with ERA-Interim variables during the overlapping historical period. On a seasonal scale, the MS governing predictors were more diverse, indicating their poor representation by most GCMs. Some earlier studies (e.g., Ashfaq et al., 2017) also reported major shortcomings of the GCMs in simulating MS dynamics over this region. We argue that a high reference uncertainty may restrict an accurate model assessment during the MS. In contrast, the available GCMs better simulated westerly-dominated seasons that account for more than two-thirds of the HA precipitation (Hewitt et al., 1989). Our analysis further showed that no single model could effectively simulate the basin-wide precipitation during different seasons. However, within a season, our model ranking process could identify models that provide improved simulations for influential predictors over multiple sub-regions. Such better-performing models can also guide the selection of driving GCMs in the scope of dynamic downscaling in this region.
Concerning future changes, the ensemble medians showed an elevation-dependent response of the UIB towards projected warming, where HA regions will mostly receive more precipitation. These positive signals were distinct during the main precipitation seasons (i.e., WS and MS), particularly over the central Karakoram. A decreasing precipitation pattern also emerged during the PMS, particularly towards the northwestern regions. These positive and negative signals further intensify under RCP8.5 at the end of the 21st century.
The projected spatial changes suggest a strengthening (weakening) and further northward penetration of the westerly system during the WS (PMS). Similarly, the future MS circulations will also intensify and penetrate further into the northwestern and trans-Himalayan regions. However, the LI regions showed a mixed seasonal response, where MS precipitation will primarily increase.
We also evaluated the impact of model weighting on ensemble signals, which generally was small and more prominent during the MS. In many cases, the model weighting refined the change signals by assigning more importance to better-performing models. For example, the weighting indicated a further strengthening of the WS westerlies and MS circulations and increased dryness during the PMS under RCP8.5. As the GCMs that perform better during the historical period may also provide reliable projections (e.g., Shukla et al., 2006), using historical weights in an ensemble setting can be advantageous (e.g., Kaspar-Ott et al., 2019).
We used SNR to assess the robustness of change signals and demonstrated the significance of positive changes over observational uncertainty, particularly at HA regions. The HA-UIB analysis further supported elevation-dependent precipitation changes. Our precipitation projections support augmentation of the existing cryosphere and continuing anomalous behaviour of the Karakoram (e.g., Bashir et al., 2017) even under the extreme warming scenario during the end of the 21st century.
Although we used predictor, predictand, and model level considerations to improve the quality of projected precipitation, some issues may still affect our analysis. For example, the absence of a large number of CMIP5 models and the existence of inter-model similarities in our ensemble may induce additional uncertainties. Although we demonstrated the usefulness of our 8-model ensemble through literature review and sample analysis, additional GCMs may differently model other precipitation predictors. Stationarity assumption may induce uncertainties not estimated in our analysis. Quantifying such uncertainties is difficult, but some studies (e.g., Merkenschlager et al., 2017) suggested mechanisms of non-stationarity considerations in statistical downscaling. Such analysis is, however, beyond the scope of the present work.
We also assume our regionalization scheme's effectiveness to extend precipitation inferences beyond the observations and draw conclusions over the transboundary basin regions not covered in our analysis (Pomee et al., 2020). Since our projection patterns are primarily in line with those studies that employ glacial and transboundary information (e.g., Lutz et al., 2016b), our assumptions regarding inferences beyond the observations seem justified.
Despite some limitations, our study presented an alternative and realistic perspective to assess precipitation changes over a complex, highly uncertain, and yet enormously important river basin. Our approach will open new avenues of regional research and has the potential for replication in other regions. However, for the correct assessment of future water availability and cryosphere stability, temperature modelling is similarly required.
ACKNOWLEDGEMENT PARC mainly supported our research work through the Himalayan Adaptation, Water, and Resilience (HI-AWARE) project, funded by the CARIAA consortium. The German Research Foundation supported EH under project number 408057478, and DAAD through Augsburg University provided valuable financial support for Muhammad Saleem Pomee. The authors also acknowledge the World Climate Research Programmer's Working group on the Coupled Modelling, ECMWF, and NOAA-ESRL Physical Sciences Division for publically sharing CMIP5 and reanalysis datasets.
Open access funding enabled and organized by Projekt DEAL.

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.