Application of Machine Learning Techniques to Delineate Homogeneous Climate Zones in River Basins of Pakistan for Hydro-Climatic Change Impact Studies

: Climatic data archives, including grid-based remote-sensing and general circulation model (GCM) data, are used to identify future climate change trends. The performances of climate models vary in regions with spatio-temporal climatic heterogeneities because of uncertainties in model equations, anthropogenic forcing or climate variability. Hence, GCMs should be selected from climatically homogeneous zones. This study presents a framework for selecting GCMs and detecting future climate change trends after regionalizing the Indus river sub-basins in three basic steps: (1) regionalization of large river basins, based on spatial climate homogeneities, for four seasons using different machine learning algorithms and daily gridded precipitation data for 1975–2004; (2) selection of GCMs in each homogeneous climate region based on performance to simulate past climate and its temporal distribution pattern; (3) detecting future precipitation change trends using projected data (2006–2099) from the selected model for two future scenarios. The comprehensive framework, subject to some limitations and assumptions, provides divisional boundaries for the climatic zones in the study area, suitable GCMs for climate change impact projections for adaptation studies and spatially mapped precipitation change trend projections for four seasons. Thus, the importance of machine learning techniques for di ﬀ erent types of analyses and managing long-term data is highlighted.


Introduction
Climate change has emerged as a major driving force behind Pakistan's flooding over the past several decades, mainly by affecting snow melting and interfering with summer monsoon patterns [1]. Previous studies have also shown a growing interest in an accurate and calculated evaluation of climate change that can provide a rational system to adapt to a rapidly changing environment. Climate change impact analyses use long records of climate data from numerous gauging sites to estimate the historical and projected climate change trends. There is a growing institutional and social need to predict the potential impacts of climate change and the damage due to floods and droughts owing to spatio-temporal variability in the climate, specifically the rainfall. Previous studies have established that spatio-temporal dynamics in atmospheric covariates are strongly linked to the precipitation characteristics on a regional scale. The atmospheric circulation covariates are important in defining the precipitation patterns across the region [2]. Different atmospheric oscillations are responsible 1. comparing the spatio-temporal average of historical/baseline observational and GCM data of different climate variables [22,25,26], using several performance indicators [27][28][29]; however, changes in the spatio-temporal average variable values may not represent the change in extreme values [30]. 2.
comparing the observational and GCM data at each grid point [21,31,32] and then selecting the GCM in the entire study area using different filtering approaches, such as clustering hierarchy [33,34], Bayesian weighting [35], weighted skill score [15,36] and spectral analysis [37]; however, such selection processes do not cater to the area-specific efficiencies of the climate models because this efficiency varies amongst models and regions [4,27]. Several studies [21,31,32] have shown that the evaluated GCMs only account for some grid points over the entire study region. However, these studies have overlooked the ability of GCM to reproduce spatial climatic patterns. The drivers of this spatial variability are atmospheric circulation patterns that depend on the geographical location of the area. The capacity of the GCMs to reproduce the patterns of climate parameters spatially, in the baseline/historical period is, therefore, imperative [7,21]. The GCM or ensemble of GCMs selected, via the aforementioned comparative analysis of spatio-temporal averages and filtering techniques, may be appropriate for simulating the spatial climatic patterns for some percentage of the entire study area but not all, thus enhancing the uncertainties in projecting the climate. To narrow this uncertainty in projecting the climate variables (precipitation) trends, we propose a different method of the selection of GCMs by first identifying the climatologically homogeneous zones and then selecting the GCM in each zone through past performance assessment. These homogeneous sub-regions are based on similar spatial and temporal climate patterns, as depicted by the highly resolved observation data.
Following the philosophy of McSweeny et al. [4] for the regional suitability of the GCMs, we selected the GCMs in every homogeneous precipitation region and validated our selection through the comparative analysis with Simple Composite Method (SCM). The SCM is broadly used for generating the multi-model ensemble that is, to obtain the equally weighted mean of all the ensemble members data at each grid point [7].
With the support of machine learning techniques, long records of climate data, from numerous gauging sites and web sources, can be easily analyzed and used to determine historical and projected trends of climate change. In the present study, the modules using machine learning techniques were developed throughout the various steps of the framework. The framework comprises five important steps: (1) developing homogeneous precipitation zones using highly resolved observation data; (2) the selection of the best suitable GCM for each zone based on the estimated correlation coefficient using the GCM data and highly resolved observed data in every homogeneous climatic region; (3) comparison and validation of the selected GCM; (4) sampling the daily precipitation data for the projected period ; based on the forcing scenarios and (5) performing future precipitation trend detections. The outcomes of this framework are the estimated divisional boundaries of the climate zones, suitable GCMs for the climate impact studies and seasonal rainfall trend projections. It is necessary to mention that the results obtained using the present framework are subject to specific assumptions: (1) daily and 0.25 • temporal and spatial scales of the data have been used for climate zoning (2) seasonal and 0.25 • , temporal and spatial scales of the data have been used for the climate projection trends (3) uncertainty of the results depends on the uncertainty involved in the climate variable outputs from the GCMs. It is very important to understand the limitations and uncertainties involved in the simulation outputs of different GCMs. The present study does not aim to quantify the uncertainties present in the GCMs' simulations but only uses the GCMs' outputs for their selection in each climate zone and trend projections. Nonetheless, the framework can be replicated, using other versions of ensembles of GMCs' with finer resolution and accuracy. The scientific data associated with the present framework can be shared on request.
This paper is structured as follows-Section 2 describes the study area, stations and data details; Section 3 provides a step-by-step presentation of the methodology; Section 4 presents the results and discussion; and Section 5 presents the conclusions.

Study Area
The basin areas of the Jhelum and Chenab Rivers are 33,330 and 67,515 km 2 , respectively. The elevations in the two basins vary between 146 and 6915 m. Figure 1 shows a topographic map with APHRODITE grid points of the basins of the Jhelum and Chenab river in Pakistan. The Southwest Monsoon triggers approximately 40% of the annual precipitation, whereas the remaining~60% of precipitation is due to Westerly aggravations [24]. Figure 2 presents the average annual seasonal precipitation in various parts of the two basins. The spatially and temporally variable climatic conditions, atmospheric circulation, advected moisture and topography, with high elevations in the north and flat plains in the south, demand a flexible framework for the evaluation of the climate change effects and future predictions. Diminishing water resources due to high hydro-climatic variability and extreme climate events in Pakistan are some of the factors that inspired the research community to investigate sustainable solutions after analyzing future needs and availability.

APHRODITE Data
To delineate the climate zones in the study area, a dense network of data stations was required. There has been an increasing trend of using easily available gridded precipitation datasets for hydrologic and climatic assessments [38][39][40]. Previous studies have identified the performance of the APHRODITE dataset (version V1101) as the best-gridded product over the high mountainous regions of Asia [41]. To justify the use of the gridded dataset of APHRODITE, we performed the comparative analysis of the APHRODITE dataset and European Reanalysis gridded dataset (ERA5)

APHRODITE Data
To delineate the climate zones in the study area, a dense network of data stations was required. There has been an increasing trend of using easily available gridded precipitation datasets for hydrologic and climatic assessments [38][39][40]. Previous studies have identified the performance of the APHRODITE dataset (version V1101) as the best-gridded product over the high mountainous regions of Asia [41]. To justify the use of the gridded dataset of APHRODITE, we performed the comparative analysis of the APHRODITE dataset and European Reanalysis gridded dataset (ERA5)

APHRODITE Data
To delineate the climate zones in the study area, a dense network of data stations was required. There has been an increasing trend of using easily available gridded precipitation datasets for hydrologic and climatic assessments [38][39][40]. Previous studies have identified the performance of the APHRODITE dataset (version V1101) as the best-gridded product over the high mountainous regions of Asia [41].
To justify the use of the gridded dataset of APHRODITE, we performed the comparative analysis of the APHRODITE dataset and European Reanalysis gridded dataset (ERA5) [42] and Global Metrological forcing dataset for land surface Modeling (GMFD) [43] with the observed dataset and monthly temporal scale, at 11 gauging stations at different altitudes, as shown in Figure 1. ERA5 is also popular for its accuracy. The results of the comparative analysis have been shown in Table 1. The Pearson correlation coefficient and Kolmogorov Simirnov (KS) test (the method is discussed in Section 3.3.2) were used to compare the monthly precipitation datasets. The results tend to show agreement with the APHRODITE dataset with higher correlation coefficients at all the stations. For the KS test, p-values greater than 0.05 have been shaded, depicting the null hypothesis of similar probability distribution is not rejected. The results of comparative analysis suggests that the APHRODITE datasets have higher correlation with the observed data that is, more than 0.8 in nearly all the observed gauging station as shown in Table 1, so the APHRODITE dataset [40] of 0.25 • × 0.25 • spatial resolution for the period of 1970-2005 has been used at a daily temporal scale, for the delineation of climate zones. Nevertheless, the regionalization framework is equally applicable to the other datasets having comparatively higher resolution and accuracy.  Figure 2 presents the total seasonal average precipitation sampled from APHRODITE. The daily gridded precipitation long term data from 1951 and onwards, is provided by the APHRODITE product (V1101). It provides the data at the continental scale, including a dense rain gauge data network of Asia, comprising the Middle East, South and Southeast Asia with valid stations of 5000 to 12,000 [40].

NEX-GDDP-GCMs-CMIP5 Data
For future climate change analysis, this study used newly developed NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset. The experiments included in Coupled Model Intercomparison Project 5 (CMIP5) were devised to address the research interrogations in AR4 of the IPCC [44]. The original resolution of most GCMs in the CMIP5 is >100 km, where such coarse resolution cannot provide fine-scale information for impact studies and decision-making at a local scale. The NEX-GDDP dataset is a collection of 21 GCMs. These datasets are bias-corrected and downscaled to a finer resolution (0.25 • × 0.25 • ) using the bias-corrected spatial disaggregation (BCSD) method [45]. The Global Meteorological Forcing Dataset (GMFD), developed by Princeton University, provides the gridded observed climate data and the NEX-GDDP after bias correction and downscaling [46,47]. Table S1 lists the 21 GCMs in the NEX-GDDP, along with information on the research centers that produce the GCMs.
NEX-GDDP data consists of daily precipitation and minimum and maximum temperatures from historical periods  and future projections . The future projections are available for two representative concentration pathways (RCPs), which are global greenhouse gas emissions scenarios. These scenarios were employed by the IPCC fifth assessment report (AR5) [44]. The NEX-GDDP dataset is publicly available. The comprehensive details to use these scenarios dataset for climate change impact evaluations and adaptations have been provided by IPCC [48]. The CMIP5 data provides multi-model datasets of climate variability, which was used to develop AR5 [44].

Methodology
Machine learning algorithms are becoming more popular owing to their ability to identify the patterns and variances in large datasets composed of multivariate atmospheric covariates, location parameters, meteo-climatic information, tele-connection indices and attributes that influence precipitation [2,49,50] or the precipitation and temperature statistics [51,52]. In this study, we used the Python module scikit-learn [53], which integrates a broad range of machine learning algorithms for supervised as well as unsupervised learning. As aforementioned, this framework consists of five steps: delineation of homogeneous climate zones; selection of GCMs; validation of selected GCMs; selection of forcing Scenarios to be used for the predictions of the climate change; and climate change projections and trend detection. Figure 3 presents the methodology flow diagram adopted in this study. Further details of each step are discussed in the subsequent sections.
(BCSD) method [45]. The Global Meteorological Forcing Dataset (GMFD), developed by Princeton University, provides the gridded observed climate data and the NEX-GDDP after bias correction and downscaling [46,47]. Table S1 lists the 21 GCMs in the NEX-GDDP, along with information on the research centers that produce the GCMs.
NEX-GDDP data consists of daily precipitation and minimum and maximum temperatures from historical periods  and future projections . The future projections are available for two representative concentration pathways (RCPs), which are global greenhouse gas emissions scenarios. These scenarios were employed by the IPCC fifth assessment report (AR5) [44]. The NEX-GDDP dataset is publicly available. The comprehensive details to use these scenarios dataset for climate change impact evaluations and adaptations have been provided by IPCC [48]. The CMIP5 data provides multi-model datasets of climate variability, which was used to develop AR5 [44].

Methodology
Machine learning algorithms are becoming more popular owing to their ability to identify the patterns and variances in large datasets composed of multivariate atmospheric covariates, location parameters, meteo-climatic information, tele-connection indices and attributes that influence precipitation [2,49,50] or the precipitation and temperature statistics [51,52]. In this study, we used the Python module scikit-learn [53], which integrates a broad range of machine learning algorithms for supervised as well as unsupervised learning. As aforementioned, this framework consists of five steps: delineation of homogeneous climate zones; selection of GCMs; validation of selected GCMs; selection of forcing Scenarios to be used for the predictions of the climate change; and climate change projections and trend detection. Figure 3 presents the methodology flow diagram adopted in this study. Further details of each step are discussed in the subsequent sections.

Climate Zoning
Climate zoning/regionalization is the most important step in this framework and the basis for the subsequent steps. For regionalization, solid evaluations of different precipitation attributes require accessibility to long records of chronicled estimations at various stations inside an area [2]. 1

Climate Zoning
Climate zoning/regionalization is the most important step in this framework and the basis for the subsequent steps. For regionalization, solid evaluations of different precipitation attributes require accessibility to long records of chronicled estimations at various stations inside an area [2]. For climate zoning, 35 years' (2006-2075) daily precipitation data of APHRODITE were used. These data were resampled seasonally to develop the climate zones for each of the seasons (defined in the subsequent sections). Regionalization was performed by grouping the rain gauges with homogeneous precipitation statistics in an area [54]. This step enables partitioning of the entire area into several climate zones for every season.
The regionalization of the precipitation statistics has numerous applications in various water resource management fields, such as agriculture practices, spatial and temporal rainfall patterns, hydrological analysis, extreme event forecasting [55] and watershed management. There are several methods available to delineate climate regions, for example, subjective and objective partitioning, geographical convenience and multivariate analysis [54,56,57]. The arbitrary and slightly misleading demarcation approach for a region, which is based on administrative boundaries and physical and geographical groupings, is referred to as the geographical convenience method. The subjective partitioning method is based on homogeneous statistical characteristics of the rainfall and previous knowledge of the region. Objective partitioning demarcates a region by grouping the sites of similar climates. Principal component analysis (PCA), clustering techniques and correlation analysis are examples of multivariate analysis techniques, which are widely used for climate zone delineation [51,58,59]. These analysis techniques require the development of a large matrix that defines the characteristics of the climate in the region. We adopted the method of PCA and agglomerative hierarchical clustering (AHC) to group the sites of homogeneous precipitation. These groups of stations were validated using different cluster validity indices.

Seasonal Data Resampling
The daily precipitation data in the APHRODITE time series, at 138 Grid stations in the study area from 1975-2005, were resampled as hydrological seasons: warm-wet ("July, August and September"), cold-dry (October, November and December), cold-wet (January, February and March) and warm-dry (April, May and June). The daily precipitation data of 35 years were resampled in a seasonal cycle using simple Python code.

Principal Component Analysis (PCA)
The benefits associated with the PCA are as follows-(1) identification of spatially coherent variations that can improve the signals, where the leading components represent the maximum variance in the patterns; (2) explanation of the covariance between the variables at various places; and (3) reduction of dimensionality, thereby providing resources to describe the variability in large spatially dimensional datasets with a reduced number of principal components [60]. PCA can present all the data in a smaller matrix and attempts to retain as much information about the data as possible. We used the dimensions from 138 APHRODITE grid point datasets for PCA. The aim was to project the data onto different orthogonal axes known as principal components. A linear transformation is performed using a symmetric covariance matrix, developed from the dataset, into the principal orthogonal components. The direction of principal components is represented by eigenvectors and the eigenvalues represent the magnitude of the stretch of the axis. This gives the direction (eigenvector) and magnitude (eigenvalue) of the spread of the dataset. The leading principal component, which presents the maximum spread of the dataset, is used for further analysis. Several principal components were selected for use in hierarchical clustering using the scree plot, which shows the variance explained in percentage by each of the principal components. Each principal component has component scores for every station, based on the eigenvectors and eigenvalues. These component scores depict the climate change pattern at the respective grid point/station and can be treated as metrological parameters that are stochastically independent of each other [61]. This analysis has been performed using the PCA function of scikit-learn [53].

Agglomerative Hierarchical Clustering
The purpose of this step is to obtain the groups of sites/ stations with similar climate change patterns depicted by the component scores of the leading principal components obtained in the previous step. This is achieved using different clustering algorithms. Different behaviors in the clustering algorithm can be expected based on the data features, dimensions and input variable values. Traditionally, different clustering algorithms are used in previous studies; there is extensive literature regarding these clustering techniques [51,[61][62][63]. Recently, the agglomerative method of hierarchical clustering [61,64] has received increased attention in climate literature. In this method, smaller clusters are developed according to a bottom-up approach, followed by a sequential combination with larger clusters depending on the Euclidian distance [65] between the clusters.
In the algorithm, each evaluation (climate change pattern) was allocated to its self-cluster. Subsequently, the iterations in the algorithm were identified and joined with the closest cluster. This agglomeration continued until the formation of one cluster. A tree-like structure (dendrogram) of similarity was formed because of the hierarchical clustering, which provides meaningful information regarding the correlations/distance among different clusters. The validity of the number of clusters (CN) Appl. Sci. 2020, 10, 6878 8 of 26 was determined prior to the identification of the maximum valid Euclidian distance. The agglomerative clustering algorithm of scikit-learn [53] was utilized for grouping the sites for each season, in this study.

Optimal Clustering and Climate zone Formation
In this step, the optimal CN of stations for each season is identified. Optimal clustering validity methods must be employed for the real grouping of datasets. Comprehensive studies have been performed for the comparative analyses of the cluster validity indices. The silhouette index (S) [66], S Dbw index (S Dbw) [67] and Calinski-Harabasz index (CH) [68] have been used in this study following the recommendations on these comparative evaluations [67,69].
Using the average of the three values identified through the aforementioned validity indices, the optimal CN was identified. The numbering and identification of stations in each cluster were performed via the truncation of the dendrogram. The truncation (cut-off) bar is kept at the valid maximum Euclidian distance corresponding to the validated value of CN [20]. This was performed using a supervised learning algorithm in scikit-learn [53].

Silhouette Score
The silhouette score method [66] is employed in the analysis to obtain the optimum CN. The silhouette score is the mean of the distance between clusters. In this study, the silhouette score yielded a maximum of 14 clusters for climate zoning for all seasons. The CN corresponding to the maximum silhouette score is the basis for the decision. The Silhouette score (S) is obtained as and where the number of clusters is denoted by CN; C i represents the ith cluster; n i represents the number of objects in C i ; c i denotes the center of C i ; and d(r,s) is the distance between r and s [66].

S Dbw Validity Index
S Dbw considers the inter-cluster density to calculate inter-cluster separation [67]. One of the densities of each cluster pair must be greater than the midpoint density for the cluster centers. This index is the summation of the scatter in the clusters and the density between the clusters. Previous studies have investigated the validation properties of different cluster validity indices using synthetic experimental data [70]. They conclude that S Dbw performed best concerning the aspects of noise, monotonicity, density, skewed distribution and sub-clusters. Equations (3)-(5) can be used to calculate the density between clusters, the scatter in the clusters and the S Dbw, respectively: Appl. Sci. 2020, 10, 6878 where D is the dataset; CN is the number of clusters; σ(C i ) is the variance vector of C i , σ(D) is the variance in the dataset, (c i , c j ) is the distance between c i and c j and

Calinski-Harabasz Index
This cluster validation scheme depends on the average of the cluster sum of squares and the average of the sum of squares within clusters and is known as CH [68]. It is obtained by formula as shown in Equation (6) Here, n represents the number of points/objects in D, c denotes the center of D and d (x, c i ) denotes the distance between x and c i .

Climate Zone Polygon Formation
The realistic partition of climatic stations for coherent climate regions is a result of this method. The delineation of these regions depends on the similarity among the statistical climate characteristics of all the stations present in a region.
The clusters of stations for every season were plotted using the ArcGIS (Environmental Systems Research Institute, CA, USA) interface. After the identification of a group of stations in the different clusters, approximate cluster boundaries were drawn on the maps using the Arc-map tool to enable a clearer presentation.

GCM Selection
The criteria used to select the GCMs are the availability of the latest generation of GCMs, good spatial resolution, past performance of GCM to replicate the historical data and the representativeness of the GCM for a wide range of climatic variable (precipitation) projections [7]. The most commonly adopted criteria are the assessment and selection of GCMs based on their capability to simulate the historical and present climate, which have been adopted by numerous studies [8][9][10]. An efficient and sensible selection of GCMs is required to generate reliable and diverse meteorological inputs for the impact models. Most previous studies show that a past performance assessment is among the most effective methods for the selection of GCMs, as the GCMs thus selected can be better predictors of future climatic conditions [31].
Several selection methods have been reported in the literature. Aghakhani et al. [71] extracted GCM data at four points in the vicinity of observation stations and performed a comparative analysis between the GCM and the observed data. They used the averages of the time-series data at each station for their analysis. Xuan et al. [72] selected GCMs using an average of the climate parameters for the entire study area in the Zhejiang Province of Southeast China. Najeebullah et al. [32] ranked the GCMs at every grid point after re-gridding the GCM data at the same spatial resolution, comparing the data with the gridded data product of the Asian Precipitation -Highly-Resolved Observational Data Integration Towards Evaluation (APHRODITE) for GCM performance evaluations in all of Pakistan. When comparing the best GCM with the highly ranked GCM at each grid station, they found a significant difference between the two with respect to precipitation. Maxino et al. [15] analyzed the climate models from the Fourth Assessment Report (AR4) of the United Nations Intergovernmental Panel on Climate Change (IPCC) for two regions in the Murray Darling Basin, divided based on rainfall classification. For different climate variables, they demonstrated that the models that are flawed in one region may be better for another region; model selection should be area specific. They calculated the skill scores associated with the different climate models in the study region. Lutz et al. [22] averaged the climate data over 2.5 • × 2.5 • grid cells in three river basins (i.e., the Indus, Brahmaputra and the Ganges) and presented their model with the selected GCMs. Latif et al. [73] also suggested certain Coupled Modelled Intercomparison Project 5 (CMIP5) GCMs after evaluating their performance using seasonal rainfall spatial correlations in the Indo-Pak region. Ahmed et al. [7] evaluated the spatial accuracy of the CMIP5 GCMs using different spatial metrics for all of Pakistan. They performed the analysis after re-gridding the GCM data into a 2 • × 2 • grid size. Srinivasa et al. [31] ranked the GCMs for the maximum and minimum temperatures across India. This framework comprised the evaluation of GCM suitability at every grid point and then ranking the GCMs using compromised programming.
In the present study, after defining the clusters of stations with homogeneous precipitation, one representative station was selected in every climate zone, using quota sampling [61]. The selection of the representative station was done based on the average climate signal climate signals (Component Scores) in a respective climate zone. The comparative analysis of GCMs data and APHRODITE data at the representative station of each climate zone was accomplished using the coefficient of determination (R 2 ).
A reduced correlation existed within the daily outputs of the GCMs and the in situ data. However, better correlations were obtained when seasonal data in the GCMs were compared with the seasonal observation data. Therefore, using the seasonal data, the correlations between the GCM data and observed gridded data were evaluated at each representative station in the respective climate zone.

Validation of Selected GCMs
This is the important step of comparative analysis of the climatic data generated through the present selection method of the GCMs and the contemporary method of sampling through simple mean based multi-model Ensemble (MME) data. The KS test [74], a nonparametric test, was used to determine the validity of the selected GCM by detecting the changes in the distribution of the two datasets (dataset generated through present selection method and simple mean based MME data) with the APHRODITE. For every grid station, the observational precipitation data (APHRODITE), selected GCM simulated data of precipitation and multi-model ensemble(MME) mean daily precipitation data of all the 21 GCMs of NASA Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) CMIP5 for the same baseline period from 1970-2005 were used in this test.

Seasonal Data Sampling
The daily precipitation data of the selected GCM, corresponding to each homogeneous climate zone, were seasonally resampled for every grid station (i.e., the grid station within each grid cell). The multi-model ensemble mean data were generated by resampling the data as the average of the GCM daily precipitation data at every grid station.

Kolmogorov-Smirnov Test
In this study, the KS test stats were derived as the highest vertical difference between the two cumulative distribution functions (CDFs) of a time series data. For CDF of APHRODITE data ø n (x) and the CDF of the selected GCM or MME mean at a grid station ø (x), the KS test statistics were obtained as given in Equations (7) and (8). where Here, sup x is the upper bound of the set of distances. The indicator function is represented by I x i ≤x (If x i ≤ x, I x i ≤x equal to 1 or if x i > x I x i ≤x equal to 0). The null hypothesis of a similar distribution is rejected if S n ≥ 0.05 significance level. The KS test was applied at every grid station using scikit-learn.

Selection of Forcing Scenarios
Four climate forcing scenarios are commonly used in climate research for climate simulation studies on both long term and short term scales. There are four common Representative Concentration Pathways (RCPs)-RCP 2.6 is a mitigation scenario; RCP 4.5 and RCP 6.0 are scenarios of medium stabilization and RCP 8.5 is a high baseline emissions scenario [75]. We did not include RCP 2.6 for the ensemble of climate models, as it is not considered a realistic or practical scenario necessary to promote adaptation planning. These rest of the three scenarios represent the complete extent of radiative forcing. In this study, we used the RCP 4.5 and RCP 8.5. The same framework applies to other forcing scenarios [22].
Extraction of Selected GCMs Data for RCP 4.5 and 8.5 The daily precipitation data for the selected GCM in the specific climate zone for the projected period from 2006-2099, as well as forcing scenarios RCP 4.5 and 8.5, were sampled for all seasons at every grid station.

Climate Change Projections and Trend Detection
The highly correlated GCMs, which can replicate the gridded dataset of seasonal precipitation for the period from 1970-2005, were selected in each climate zone for use in the future predictions of precipitation change trends in RCP 4.5 and RCP 8.5. The precipitation change trends were detected using the Mann-Kendall (MK) test during the projected period (2006-2099). The trend was quantified by employing Sen's slope estimator.

Mann-Kendall Test (MK Test)
The MK test [76] was used to evaluate the statistically significant annual seasonal trends in the long term precipitation data. The "no trend" in precipitation with time is assumed in the null hypothesis (H o ) of the test and vice versa for the alternative hypothesis (H a ). The Equations (9)- (12) show the test stats T of the MK test.
where D i and D j are the consecutive time-series observations, which are organized chronologically for the length of data (n); t p represents the data points number for the pth value in a tied group; and q denotes the total number of tied groups σ is the variance. An upward time series trend is represented by a positive Z value and vice versa for the negative trend. If |Z| > Z 1−α/2 , the null hypothesis (H o ) is rejected which is indicative of a statistically significant trend. The critical value of Z 1−α/2 is assumed at the p-value of 0.05.

Sen's Slope Evaluation
The Python module was applied to obtain the slope of the trend, if exists, according to Sen's method [77]. According to the method, linear slope sets can be estimated as shown in Equation (13).
where Ti is the slope, D j and D k are the variables at j and k time steps, respectively and n represents the data points number. Sen's slope is estimated as the median of all slopes.

Results and Discussion
Taking the spatiotemporal average [22,25,73] of the climate dataset at all the grid points of the area at different temporal scales (daily, seasonal or annual) or using various spatial metrics on the results of the individual grid points [21,31,32] are the conventional methods used for the selection of GCMs. However, these methods do not address the spatial coherence of the climate variability patterns in the study area. This poses uncertainty in using these conventional methods. Several studies support the association between atmospheric covariates and precipitation. In GCMs, these atmospheric covariates are the input variables. To address the uncertainty arising from the input variables of the GCMs, it is imperative to assess the GCM performance in reproducing the regional climate pattern. These results show that the climate variability pattern of the study area is spatially heterogeneous in all seasons and a single or an ensemble of GCMs may not represent this spatial climatic heterogeneity. Therefore, we classified/ regionalized the large study in several climate zones founded on the climate variability patterns similarity and the selection of a GCM should be done separately in each of the homogeneous climate zones.
This section presents the detailed results along with the discussion. The results have been visualized using ArcGIS. This software is based on the integration of different machine learning techniques. Section 4.1 describes the results of the analyses performed for the development of climate zones. Section 4.2 shows the results of the GCM selection procedure and Section 4.3 is related to the validation of the method used for the GCM selection. Section 4.4 described the results of the MK test of trend detection and Sen's slope estimation. These tests determined the precipitation change trend projections for the period of 2006 to 2099.

Climate Zones
The Principal Component Analysis (PCA) performed as the first step in climate zoning provided the details of the climate change pattern in the entire study area. A statistically significant climate pattern heterogeneity was observed in the region when PCA was performed on a seasonal basis.

Principal Component Analysis
The application of the Principal Component Analysis (PCA) was performed on the daily precipitation series data from APHRODITE on a network of 138 stations distributed across the Jhelum and Chenab river basins. This first led to 20 primary and physically significant principle components (PCs), which collectively explained approximately 95% of the total variability in the precipitation throughout the study area. These were retained for the cluster analysis. Scree Plot is shown in Figure 4. The PCA algorithm produced the scree plots, allowing us to select the number of PCs for the cluster analysis. For the warm-wet season, 20 PCs present 95% of the data variance from the station network (Figure 4a). For the cold-dry season, 20 PCs present 94% of the variance in the data from the station network (Figure 4b). For the cold-wet season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4c). For the warm-dry season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4d). The spatial distribution of the component scores is shown in Figure 5, for the first two PCs, explaining the highest percentage of variance in the data.
Jhelum and Chenab river basins. This first led to 20 primary and physically significant principle components (PCs), which collectively explained approximately 95% of the total variability in the precipitation throughout the study area. These were retained for the cluster analysis. Scree Plot is shown in Figure 4. The PCA algorithm produced the scree plots, allowing us to select the number of PCs for the cluster analysis. For the warm-wet season, 20 PCs present 95% of the data variance from the station network (Figure 4a). For the cold-dry season, 20 PCs present 94% of the variance in the data from the station network (Figure 4b). For the cold-wet season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4c). For the warm-dry season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4d). The spatial distribution of the component scores is shown in Figure 5, for the first two PCs, explaining the highest percentage of variance in the data.   components (PCs), which collectively explained approximately 95% of the total variability in the precipitation throughout the study area. These were retained for the cluster analysis. Scree Plot is shown in Figure 4. The PCA algorithm produced the scree plots, allowing us to select the number of PCs for the cluster analysis. For the warm-wet season, 20 PCs present 95% of the data variance from the station network (Figure 4a). For the cold-dry season, 20 PCs present 94% of the variance in the data from the station network (Figure 4b). For the cold-wet season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4c). For the warm-dry season, 20 PCs explained 95% of the variance in the data from the station network (Figure 4d). The spatial distribution of the component scores is shown in Figure 5, for the first two PCs, explaining the highest percentage of variance in the data.   These two PCs cumulatively produced variance of 44.8% in the warm-wet season, 41.2% in the cold-dry season, 45.7% in the cold-wet season and 41.5% in the warm-dry season. For the component scores, the first PC for the warm-wet season, as shown in Figure 5a (with 25% of the explained variance), has high negative scores in the north and southwest of the study region and medium-high positive signals in the central region. Stronger positive and negative signals exhibit an affinity of higher variation in the precipitation amount, while weaker scores indicate low variability. The second PC for the warm-wet season, as shown in Figure 5e (with 19.8% of the explained variance), has mostly positive component scores in the southeast, except for medium and strong negative signals in several southeast and central regions, respectively. In the cold-dry season, we observe higher negative variability in the southwest based on PC1, as shown in Figure 5b and in small areas of the southeast. Strong positive signals were detected in the northwest area of the study region. For PC2, high positive component scores were observed in the southeast, as shown in Figure 5f (with 18.8% of the explained variance), whereas strong negative scores occurred in the southwest area. In the cold-wet season, the variability in the precipitation amount was high in the southwest region for PC1 (with 27% of the explained variance in the data), as shown in Figure 5c and high positive signals in the southeast and central areas of the region. For PC2, cold-wet season (with 18.7% of the explained variance), shown in Figure 5g, had strong positive scores in the southeast area of the region but medium to high negative scores in the rest of the study area. In the warm-dry season, high variability was observed in the precipitation in the northern and southwest regions of the area for the PC1 component score (with 22.7% of the explained variance), as shown in Figure 5d. The PC2 had high positive variability in the southwest and medium negative signals in the rest of the region, as shown in Figure 5h. The first two PCs were noted to explain a significant variability for the precipitation data in all seasons.

Agglomerative Hierarchical Clustering (AHC)
We then identified all the potentially homogeneous regions inside the study area. The extent of the homogeneity in these homogeneous regions was identified using different cluster validity indices. The first 20 PCs were employed in the AHC analysis. These PCs were obtained in the previous step using the Python scikit module. The validity of different site/station clusters was evaluated by performing the silhouette score, CH and S Dbw tests. Figure 6 summarizes the cluster validation results. The machine learning algorithms for the cluster validity were run repetitively using different input values for the maximum Euclidean distance, followed by a comparison of the resulting clusters to obtain their validity. According to the results, the optimized partitioning of the data based on the CH Test was investigated while the lowest scores were estimated corresponding to clusters 13, 15, 14 and 14 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively, as shown in Figure 6a. The partitioning based on the S Dbw test suggests that optimized clustering can be performed when the data is partitioned into clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively, as shown in Figure 6b. These clusters were selected based on the lowest S Dbw scores. The partitioning of data into clusters 11, 15, 13 and 13 corresponds to the highest silhouette test scores for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively, as shown in Figure 6c. Based on Figure 6, all validity measures agree with the same number of clusters. Therefore, based on the recommendations of previous studies [70], we partitioned the stations for optimized clustering as clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively. The members/stations in different clusters were obtained via the truncation of the dendrogram at the Euclidean distances that correspond to the optimized clusters. The dendrograms were obtained from the AHC of the selected PCs. The maximum Euclidean distances were identified as corresponding to the proposed cluster numbers, where Figure 7 shows the cutoff/truncation bars for each season. Based on this truncation, we identified the groups of Based on Figure 6, all validity measures agree with the same number of clusters. Therefore, based on the recommendations of previous studies [70], we partitioned the stations for optimized clustering as clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively. The members/stations in different clusters were obtained via the truncation of the dendrogram at the Euclidean distances that correspond to the optimized clusters. The dendrograms were obtained from the AHC of the selected PCs. The maximum Euclidean distances were identified as corresponding to the proposed cluster numbers, where Figure 7 shows the cut-off/truncation bars for each season. Based on this truncation, we identified the groups of stations/sites in each cluster. Figure 7 shows the dendrograms that present the cluster trees for all four seasons. Figure 6, all validity measures agree with the same number of clusters. Therefore, based on the recommendations of previous studies [70], we partitioned the stations for optimized clustering as clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively. The members/stations in different clusters were obtained via the truncation of the dendrogram at the Euclidean distances that correspond to the optimized clusters. The dendrograms were obtained from the AHC of the selected PCs. The maximum Euclidean distances were identified as corresponding to the proposed cluster numbers, where Figure 7 shows the cutoff/truncation bars for each season. Based on this truncation, we identified the groups of stations/sites in each cluster. Figure 7 shows the dendrograms that present the cluster trees for all four seasons.

Climate Zones and Reference site
This study reveals useful details on the efficiency of machine learning algorithms when formulating large, homogeneous regions of precipitation for different seasons. Homogeneous precipitation sites were identified based on the validated number of cluster-values. These station clusters were then plotted on a study area map, as shown in Figure 6. All homogeneous regions are differentiated by different colors to show the general spatial patterns associated with precipitation, the approximate boundaries are drawn on the maps to allow for a clearer representation. The clusters of stations transformation into different climate zones for all seasons have been shown in Figure 8.
The Jhelum and Chenab river basins are partitioned into clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively, based on the analysis. There were three outliers in the clusters of the warm-wet season, one outlier in the clusters of the cold-dry season, one outlier in the cold-wet and one outlier in the warm-dry season. After merging the outliers with the neighboring clusters, the final demarcations of the basins were mapped for every season. The basins are finally classified into 10, 12, 13 and 13 climate zones for warm-wet, cold-dry, cold-wet and warm-dry seasons. The reference station was identified in each climate zone for the selection of the best-correlated GCM, using the past performance analysis of the GCM at the reference station in each climate zone.
formulating large, homogeneous regions of precipitation for different seasons. Homogeneous precipitation sites were identified based on the validated number of cluster-values. These station clusters were then plotted on a study area map, as shown in Figure 6. All homogeneous regions are differentiated by different colors to show the general spatial patterns associated with precipitation, the approximate boundaries are drawn on the maps to allow for a clearer representation. The clusters of stations transformation into different climate zones for all seasons have been shown in Figure 8. The Jhelum and Chenab river basins are partitioned into clusters 13, 13, 14 and 13 for the warm-wet, cold-dry, cold-wet and warm-dry seasons, respectively, based on the analysis. There were three outliers in the clusters of the warm-wet season, one outlier in the clusters of the cold-dry season, one outlier in the cold-wet and one outlier in the warm-dry season. After merging the outliers with the neighboring clusters, the final demarcations of the basins were mapped for every season. The basins are finally classified into 10, 12, 13 and 13 climate zones for warm-wet, cold-dry, cold-wet and warm-dry seasons. The reference station was identified in each climate zone for the selection of the best-correlated GCM, using the past performance analysis of the GCM at the reference station in each climate zone.

GCM Selection
We investigated the capabilities of 21 CMIP5 GCMs to reconstruct the highly resolved observation data for precipitation and, using Pearson correlation analysis, identified the most correlated GCMs for a dependable climate projection in a climate zone. The GCMs were selected in every climate zone of the study area using the GCM and APHRODITE seasonal precipitation data at the reference station. The reference stations in every climate zone for every season have been shown in Figure 8. These selected GCMs were then used for the projections of changes in the precipitation trends at every grid point throughout the respective zone in the river basins. Figure 9 shows the spatial distribution of the highly correlated GCMs in each climate zone. For the warmwet season, we selected the BNU-ESM, MIROC5, CESM, IPSL-CM5A-MR and GFDL-ESM2G GCMs. For the cold-dry season, we selected the CCSM4, MIROC-ESM, IPSL-CM5A-LR, MIROC5

GCM Selection
We investigated the capabilities of 21 CMIP5 GCMs to reconstruct the highly resolved observation data for precipitation and, using Pearson correlation analysis, identified the most correlated GCMs for a dependable climate projection in a climate zone. The GCMs were selected in every climate zone of the study area using the GCM and APHRODITE seasonal precipitation data at the reference station. The reference stations in every climate zone for every season have been shown in Figure 8. These selected GCMs were then used for the projections of changes in the precipitation trends at every grid point throughout the respective zone in the river basins. Figure 9 shows the spatial distribution of the highly correlated GCMs in each climate zone.  Table S1 in supplementary documents.

Validation of Selected GCMs
The selected GCMs and simple mean of MME (21 GCMs) daily precipitation data at every grid station were resampled for the seasonal cycles. The KS test was applied to every grid station and the data distribution was compared with the APHRODITE data distribution. p-values < 0.05 demarcate between the null and alternative hypotheses. When the p-value < 0.05, the similar distribution hypothesis is rejected and vice versa. Figure 10 shows the spatial distribution of the p-value for the KS test. The result shows that the selected GCMs present a higher degree of similarity in the daily seasonal distribution to the APHRODITE distribution as compared with the simple mean of MME distribution. The results show that 45% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), on the other hand, when using the datasets of simple mean MME, 28% of the datasets did not fall in the rejection zone for the warm-wet season. 76% of the time-series distribution of the datasets of the selected GCMs did not fall in the rejection zone (when α = 5%) for the cold-dry season, whereas only 55% of the datasets were not in the rejection zone when using the datasets of simple mean MME. For the cold-wet season, 76% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), whereas 0% of the datasets were not in the rejection zone when the simple mean MME was tested. For the warm-dry season, 73% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), whereas only traces of 10% of the datasets were not in the rejection zone when the simple mean MME was tested. The results validated the GCM selection procedure when assessed against past performance.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 27 and MPI-ESM-LR GCMs. For the cold-wet season, we selected the MRI-CGCM3, MIROC-ESM, NorESMI, CESM, IPSL-CM5A-MR and GFDL-ESM2G GCMs. For the warm-dry season, we selected the IPSL-CM5A-MR, BNU-ESM, GFDL-CM3, bcc-CSM1-1, inmcm4 and GFDL-CM3 GCMs. The description of the aforementioned GCMs can be checked from the Table S1 in supplementary documents. Figure 9. The spatial distribution of the highly correlated general circulation models (GCMs) with the observed data in each respective climate zone for the (a) warm-wet, (b) cold-dry, (c) cold-wet and (d) warm-dry seasons.

Validation of Selected GCMs
The selected GCMs and simple mean of MME (21 GCMs) daily precipitation data at every grid station were resampled for the seasonal cycles. The KS test was applied to every grid station and the data distribution was compared with the APHRODITE data distribution. p-values < 0.05 demarcate between the null and alternative hypotheses. When the p-value < 0.05, the similar distribution hypothesis is rejected and vice versa. Figure 10 shows the spatial distribution of the pvalue for the KS test. The result shows that the selected GCMs present a higher degree of similarity in the daily seasonal distribution to the APHRODITE distribution as compared with the simple mean of MME distribution. The results show that 45% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), on the other hand, when using the datasets of simple mean MME, 28% of the datasets did not fall in the rejection zone for the warmwet season. 76% of the time-series distribution of the datasets of the selected GCMs did not fall in the rejection zone (when α = 5%) for the cold-dry season, whereas only 55% of the datasets were not in the rejection zone when using the datasets of simple mean MME. For the cold-wet season, 76% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), whereas 0% of the datasets were not in the rejection zone when the simple mean MME was tested. For the warm-dry season, 73% of the time-series distribution of the selected GCMs' datasets did not fall in the rejection zone (when α = 5%), whereas only traces of 10% of the datasets were not in the rejection zone when the simple mean MME was tested. The results validated the GCM selection procedure when assessed against past performance.

Seasonal Precipitation Trend Projection
The significance of seasonal precipitation trends for the period of 2006-2099, as detected by the MK trend detection test, for forcing scenarios RCP 4.5 and 8.5, are presented in Figure 11. The distribution of p-value spatially for the MK test has been presented. Spatial patterns of the estimated Sen's slope exhibit negative and positive slopes for the trends in different seasons, as shown in Figure 12. The significance of the trends was identified via the MK test of trend detection, Figure 10. The p-values for the K-S Test comparing the APHRODITE data with selected GCMs data and conventional composite mean data of 21 GCMs for warm-wet, cold-dry, cold-wet and warm-dry season (a) Using Selected GCMs (b) Simple mean Multimodel Ensemble. The green color bands are presenting the zones where the null hypothesis of "similar distribution" is not rejected.

Seasonal Precipitation Trend Projection
The significance of seasonal precipitation trends for the period of 2006-2099, as detected by the MK trend detection test, for forcing scenarios RCP 4.5 and 8.5, are presented in Figure 11. The distribution of p-value spatially for the MK test has been presented. Spatial patterns of the estimated Sen's slope exhibit negative and positive slopes for the trends in different seasons, as shown in Figure 12. The significance of the trends was identified via the MK test of trend detection, as shown in Figure 11. For the warm-wet season, with the RCP 4.5 forcing scenario (refer to Figure 12a), 28% of the total study region yielded very strong and 5% of the total study region yielded strong evidence for decreasing precipitation at values of 2.2-3.29 mm yr −1 over the central and northwest of the area. For RCP 4.5, 67% of the region showed no or weak decreasing trends. However, for the RCP 8.5 (Figure 12e), weak evidence for increasing trends was detected in the central and southeast region at 0.5-2 mm yr −1 , except for certain areas in the central to western regions, where we detected a non-significant decreasing trend, with a slope of 1-3.29 mm yr −1 . Figure 11 depicts the p-values for the MK trend detection, which shows significant and non-significant trends.  For the cold-dry season and RCP 4.5 (Figure 12b), nearly all the regions of the study area had weak, increasing trends with Sen's slopes of 0-0.5 mm yr −1 , which are not significant, as depicted in Figure 11. On the other hand, for RCP 8.5 (Figure 12f), weak decreasing trends were detected in the central and southwest regions and certain areas of the southeast, with Sen's slopes of 0.2 mm yr −1 . Here, 80% of the total area had weak evidence of increasing trends at 0.5-1.5 mm/season/year.
For the cold-wet season for the RCP 4.5 (Figure 12c), the weak decreasing trend, with Sen's slopes of 0.2-0.5 mm yr −1 , was detected in the central and southwest parts of the area. High and weak increasing trends were detected in the northwest and southeast parts of the area, respectively. The positive slope varies between 0.5 and 2 mm yr −1 .
For the cold-wet season in the RCP 8.5 (Figure 12g), a non-significant decreasing trend was detected with Sen's slopes varying between 0.5 and 0.2 mm yr −1 in the central area, whereas in other regions, a weak, increasing trend was detected with slopes varying between 0.5 and 2 mm yr −1 . Figure 11. The p-values of the MK test for the best-correlated models for the period (2006-2099) for warm-wet, cold-dry, cold-wet and warm-dry for (a) Forcing Scenario of RCP 4.5 (b) Forcing Scenario of RCP 8.5. The blue and green shades are depicting that the Null Hypothesis of "no trend" is rejected. For the cold-dry season and RCP 4.5 (Figure 12b), nearly all the regions of the study area had weak, increasing trends with Sen's slopes of 0-0.5 mm yr −1 , which are not significant, as depicted in Figure 11. On the other hand, for RCP 8.5 (Figure 12f), weak decreasing trends were detected in the central and southwest regions and certain areas of the southeast, with Sen's slopes of 0.2 mm yr −1 .
Here, 80% of the total area had weak evidence of increasing trends at 0.5-1.5 mm/season/year. Unlike for the warm-dry season and RCP 4.5 (Figure 12d), there is strong evidence for an increase in the precipitation for the whole study area with p-values from the MK trend detection at less than 0.05 for nearly the entire study area, with the Sen's slope indicating an increase in the precipitation rate with values varying between 1.5 and 3.5 mm yr −1 . The same is the case for RCP 8.5 in the warm-dry season based on strong evidence of an increase in the precipitation throughout the entire area, with Sen's slope values varying between 1.5 and 3.5 mm yr −1 (Figure 12h). Figure 11 shows the p-values for the MK test for the best-correlated models (2006-2099) for RCPs 4.5 and 8.5.

Conclusions
Previous studies have not currently arrived at a consensus on a universally accepted method and criteria used for GCM selection [14,78]. Studies continue to investigate the methods, whether dynamic or statistical, to reduce the uncertainty in predicting climate change impacts. Intending to reduce the uncertainty, we presented a novel and flexible framework for the selection of GCMs in homogeneous climatic zones, which are based on daily seasonal precipitation data statistics spanning the baseline/historical period  for the Jhelum and Chenab river basin study areas. The GCM selected in each of the climate zones of every season can reproduce the climate distribution patterns in the study area. This has been proved by comparatively analyzing the precipitation data sets of GCMs selected and the simple composite mean of the ensemble of the 21 CMIP5 GCMs. Using these homogeneous precipitation regions, we suggest a different approach for selecting the GCMs. The GCMs were selected based on prior performance and compared with the highly resolved and gridded APHRODITE data. The GCM selection was validated for its performance to emulate the spatial precipitation patterns during the baseline period  for the study region. We sampled the daily precipitation data for the projected period using the selected GCM. The trends in the precipitation variability, based on the MK trend detection statistics, show that the regions in the Jhelum and Chenab river basins have negative, positive and no trends from 2006 to 2099.
The Jhelum and Chenab river basins are scarcely gauged regions; hence, the APHRODITE gridded datasets were used in this study because this type of climate zoning requires a dense network of climate observation stations. Numerous previous studies have verified the reliability of the APHRODITE data in this region [22,40]. However, at high altitudes, the deviation between the gridded data and in situ observations is more than that of flat regions due to a reduced number of gauging stations in mountainous areas, which poses a challenge to precise interpolation [79]. Uncertainties in the forcing datasets persist due to the scarcity of gauging stations in high elevation zones. Furthermore, gauging stations in valleys cannot represent peripheral higher elevation zones due to the high vertical lapse rate of precipitation [41]. Currently, there is no consensus among the climate research community regarding the methods of bias correction for high elevation areas. The comparative analysis shows that the APHRODITE correlates well with the observed dataset of some gauging stations at high altitudes, leading to the assumption of having sufficient accuracy for the regionalization process. Nevertheless, this regionalization framework is equally applicable to other datasets having comparatively higher resolution and accuracy.
All the analyses were conducted using the Python programming language, which is known to be powerful in machine learning. Several machine learning modules from the scikit-learn library, such as Principal Component Analysis (PCA), Agglomerative Hierarchical Clustering (AHC) and clustering validity indices, correlation coefficient and the KS test, were used in the analysis. All these machine learning algorithms were compiled to develop a program for regionalization and climate change trend detection. The program can be provided to the climate research community to augment decision-and policymaking for water resource planning and management.
Potential forecasts from the GCMs, however, are fundamentally uncertain and decision-makers still find it difficult to understand or use such predictions of climate change. The major cause of this uncertainty in the assessments derives from the inherent uncertainties in GCM outputs [80][81][82][83]. Based on the results and interpretations presented in this study, we can draw the following conclusions:

1)
Due to its highly volatile hydroclimatic conditions, this area of Pakistan poses major scientific challenges; hence, investigations require complex methods. Multivariate techniques, such as PCA and Agglomerative Hierarchical Clustering (AHC) algorithms, are used to develop decisive precipitation statistics and delineate homogeneous precipitation regions, respectively. The entire study area was divided into 13 homogeneous precipitation regions for the warm-wet season, 13 homogeneous precipitation regions for the cold-dry season, 14 homogeneous precipitation regions for the cold-wet season and 13 homogeneous precipitation regions for the warm-dry season. The reference station, which is representative of the respective homogeneous precipitation regions, was obtained in each homogeneous precipitation region for GCM selection. Seasonal rainfall characteristics were incorporated to define the homogeneous climate regions to design the climate zones in the river basins of Pakistan. Representative rainfall station/grid points were selected in each climate zone for the climate change impact assessment. This study provides an objective demarcation structure for homogeneous precipitation regions based on the statistical characteristics of seasonal precipitation. 2) The best-correlated GCM was identified in each climate zone in the Jhelum and Chenab river basins, followed by data extraction. This selection was validated and compared spatially with the KS test. Two schemes of the daily precipitation data were derived: (1) the precipitation data sampled from the selected GCMs based on this framework and (2) the precipitation data from the MME mean (21 models of the CMIP5 mean). The KS test was used to measure the compatibility of the observed data with data from scheme 1 and then the data from scheme 2. The results clearly show the improved performance of the present framework for the selected GCMs compared with the conventional method involving the equally weighted means of all the available GCMs.
3) The precipitation pattern can represent the climatic dynamics because the atmospheric circulation covariates have a strong correlation with the precipitation patterns in the region [2]. The daily projected precipitation data , synthesized for the regions of the Jhelum and Chenab river basins, are publicly available. The creation of homogeneous climatic zones accommodates a deeper understanding of the dynamic spatio-temporal precipitation variability over a given region. Infrastructure development and climate forecasting are based on effective regionalization using reliable estimates of climate variables. The region should be conditionally considered as climatically homogeneous to proceed with the selection of GCMs for the climate change impact assessments. 4) The projections of the changes in the precipitation trends were presented based on two scenarios, that is, RCP 4.5 and RCP 8.5. For the RCP 4.5, a significant positive precipitation trend was observed for the warm-dry season, while insignificant positive trends were observed for the cold dry and cold wet seasons. A negative precipitation trend was observed throughout the entire study area, in which 40% of the area had a significant negative trend for the warm-wet season.
For the RCP 8.5, the warm-dry season again exhibited a significant positive precipitation trend, whereas insignificant positive trends were observed for the warm-wet, cold-dry and cold-wet seasons. The effect of the precipitation change trends may have greater implications on water resources and its management. This unique study provides spatial changes in patterns, as well as temporal changes. Strong increasing and decreasing signals occurred during the warm-dry season for both scenarios in the river basins. Precipitation trends were detected within this framework, which were mapped spatially across the entire study region. 5) Future studies should redefine the homogeneous precipitation region and change detection using the projected data to obtain the hydrological response in a basin using this projected input data for the rainfall. The same framework, with a set of machine learning modules, can be employed to select and detect precipitation trends using the next generation of GCMs. More than 125 years of precipitation data of daily scale, were used for the analyses in this study, which depicts the strength of machine learning.
There is extensive literature available, addressing the uncertainties in climate projections through climate models [5,17]. The uncertainty in predicting the near future day to day weather of dynamic chaotic weather systems has been discussed by meteorologist Edward N Lorenz [16], in which he emphasized the theory of chaos, that is, the importance of the sensitivity of initial state simulation and parametrization of a model for plausible weather forecast. The knowledge of the precise initial state of the weather conditions is essential for the robust whether forecast, furthermore, it is also useful to predict the climate having initial conditions of weather. However, the predictability of the weather is limited to few days due to the reason of atmospheric dynamics. Weather patterns are affected by atmospheric oscillations such as El Nino and La Nina. Climate projections are done using the heat energy changes on Earth's surface and the changes in greenhouse gas emissions in the atmosphere. The projections of these changes to eventually project the climate is significantly easier than the weather forecasting of a week. To summaries, long term changes in atmospheric composition are significantly more predictable than the short term weather forecast [84].
The sources of confidence in future projections by the models are their physical basis and their performance in simulating the past/historical climate. According to IPCC, the models have been proven as important tools to project the future climate [84]. Their performance in representing large scale climate features are extensively evaluated through the comparative analysis of its simulations with the observations and they have been used to predict the ancient climate such as the warm Holocene period, the last glacial maximum and the last ice age [84].
The paper is envisioned to provide a framework that can easily be replicated to project the climate change trends for an area by utilizing the stochastic analysis strengths of machine learning and the available observed and GCMs' data. It tends to be an important reference to those who are working on impact modeling. However, it does not aim to discuss the dynamic processes and associated uncertainties present in all GCMs. These are presented in the literature and their simulation outputs have been used in the analysis. To gain a higher confidence level, the end-users of the present framework are suggested to use large ensemble of GCMs while selecting a suitable model in each zone of the study area. The codes prepared in this study are provided upon request.
Further studies are required to analyze the seasonal climate shift and how accurately the GCM outputs can project the present climatic zones, which are delineated based on the observed gridded historical data.