A robust multi-criteria and rank-stability framework for station-scale evaluation and ranking of CMIP6 GCM precipitation over the Jhelum River Basin

1 Introduction

Climate change has become a critical global concern, particularly since the latter part of the twentieth century (Abbass et al., 2022). Among hydro-climatic variables, precipitation is especially important, as its magnitude, frequency, and spatiotemporal distribution are projected to undergo substantial alterations under changing climatic conditions (Suo et al., 2025). Such changes in precipitation regimes can trigger a range of hydro-meteorological hazards, notably droughts and floods, resulting from rainfall deficits and excesses, respectively (Almheiri et al., 2024). Variations in precipitation further exert profound influences on streamflow, ground water recharge, soil vegetation interactions, agricultural and prevalent hydrological extremes (Yu et al., 2025).

In mountainous environments, this process is controlled by rugged elevation gradients, snow-rain phase and orographically induced precipitation enhancement (Cui et al., 2022). These factors lead to a strong spatial heterogeneity and temporal variability. So, making accurate representation of precipitation a fundamental prerequisite for reliable hydrological modeling, water-resources management, and climate-risk assessment. Errors in precipitation simulation in such regions can propagate through hydrological models, resulting in large uncertainties in impact assessments.

The Himalayan region is one of the most complicated regions of the world in terms of hydro-climate. They are being home to more than one-fifth of the global population. Further this region demonstrates pronounced hydroclimatic sensitivity, as cryosphere-dominated and hydrological processes are closely governed by precipitation and temperature variability (Wester et al., 2019; Sharma et al., 2022). Precipitation across Himalayas is shaped by the interaction of large-scale synoptic system, particularly western disturbances and monsoon circulations, with a strong local topographic forcing (Dahri et al., 2021; Krishnan et al., 2020).

Within this region, the Jhelum River Basin (JRB) is a highly hydrologically and socio-economically significant site, which sustains irrigation, hydropower, domestic water supply etc. The spatial and temporal trend in the precipitation of the basin is strongly variable with a combination of Western Disturbances, monsoon intrusion, valley-scale circle and localized convective activity (Ansari and Grossi, 2022; Nadeem et al., 2025; Ansari et al., 2025). These interacting processes create high heterogeneity in rainfall amount, frequency and intensity over relatively compact spatial scales and therefore posing a major challenge for climate-model simulations.

General Circulation Models (GCMs) especially from the Coupled Model Intercomparison Project Phase 6 (CMIP6) are generally utilized to study climate variability in the past and to project future hydro-climatic projections (Song et al., 2021; Ma et al., 2024). Despite advances, the problem of how to simulate precipitation in a complicated mountainous background has not been resolved completely (Lalande et al., 2021). The systematic biases that effect computer-simulated precipitation are usually due to using coarser grid spacing, simplified parameterizations of cloud microphysics and convection, and incomplete representation of fine-scale orographic lifting (Salahi et al., 2024; Prasanna et al., 2024). These deficiencies are particularly severe with regards to precipitation over temperature when considering that rainfalls are very sensitive to the minor atmospheric processes and focal convergence of moisture (Zhang et al., 2021). Consequently, the precipitation as determined by GCM often has significant deviations as compared to that of ground-based observations in high-relief areas like the Himalayas (Assamnew et al., 2025).

Precipitation in the Jhelum Basin also exhibits pronounced elevation-dependent behavior (Ahmad et al., 2022). Large scale synoptic systems dominate the conditions in low lying valley stations (Ansari and Grossi, 2022). Mid-elevation zones are greatly orographically enhanced, and high elevations zones contribute a large portion of precipitation in the form of snowfall (Dahri et al., 2016). Such elevation-dependent behavior cannot be adequately captured through basin-averaged statistics or coarse grid-scale evaluations (Malla and Arya, 2025). Therefore, station-level assessment is important to capture microclimatic variability of precipitation and to diagnosing model performance under distinct micro-climatic regimes. However, besides being a hydrological feature, station scale evaluation of climate-model precipitation in the JRB has been limited. Existing studies predominantly focus on basin-averaged or gridded analyses, often emphasizing temperature trends or mean precipitation totals (Islam and Chakma, 2024; Mahmood et al., 2018; Akhter and Manzoor Ahmad, 2017; Munawar et al., 2022b; Ahsan et al., 2023). Although these methods are informative at macroscales, the local variations caused by topography are obscured, and these methods do not reveal the dynamics of precipitation that are essential in hydrological and impact-based studies.

One more significant limitation of existing literature is the dependence on a restricted set of evaluation metrics, including correlation coefficients (CC) or root-mean-square error (Ahsan et al., 2023; Munawar et al., 2022b; Islam and Chakma, 2024). Although these metrics offer valuable information, they do not comprehensively describe the distributional structure or event-based performance of precipitation simulations (Cao et al., 2023; Deb et al., 2026). A robust assessment requires complementary indicators that jointly capture accuracy, variability, efficiency, and extreme behavior; without this, model evaluation remains fragmented and rankings can be incomplete or misleading. No clear, repeatable framework currently exists for prioritizing precipitation models at the station level, so model selection is often based on subjective judgment or limited indicators, producing inconsistent conclusions across regions and applications. This is especially consequential in mountain basins, where spatial heterogeneity means a model can perform well at one station and poorly at another.

In addition, the majority of the existing evaluations fail to specifically assess the consistency of model rankings over seasons and ways of evaluation (Ansari et al., 2023; Mahmood et al., 2018; Munawar et al., 2022a; Islam and Chakma, 2024). Rank instability introduces additional uncertainty into model selection, while it is rarely used to quantify or incorporated into final model recommendations (Rahman and Pekkat, 2024; Raju and Kumar, 2014). It is on this background, there is a clear need for a precipitation-focused, station-scale evaluation framework that integrates multiple performance metrics, employs objective multi-criteria decision-making (MCDM) techniques, and explicitly examines ranking stability. Addressing this need is important for improving the reliability of precipitation inputs used in hydrological simulations, flood forecasting, agricultural planning, and climate-change impact assessments in the Jhelum Basin. In response to these gaps, the present study develops and applies a robust, multi-criteria, and rank-stability-based framework for the station-scale evaluation and ranking of CMIP6 GCM precipitation simulations over the Indian part of the Jhelum River Basin. The specific objectives are to;

• Evaluate the performance of CMIP6 GCMs at individual meteorological stations using a comprehensive set of precipitation-relevant statistical indicators that capture accuracy, variability, efficiency, and bias.
• Examine both seasonal and annual precipitation behavior, thereby accounting for elevation-dependent climatic regimes and synoptic variability across the basin.
• Rank GCMs using multiple multi-criteria decision-making techniques [Technique for Order Preference by Similarity to Ideal Solution (TOPSIS), Preference Ranking Organization Method for Enrichment Evaluation (PROMETHEE), and VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR)] to reduce methodological bias in model selection.
• Quantify rank stability across stations, seasons, and ranking methods, enabling an objective assessment of model robustness.
• Derive an integrated consensus ranking by combining annual and seasonal rankings, thereby identifying the most reliable precipitation models for basin-scale and station-specific applications.

By moving beyond basin-averaged analyses and single-metric evaluations, this work provides a more physically meaningful and operationally relevant basis for climate-model selection. The results offer actionable guidance for hydrological modeling and climate-risk assessment in one of the most hydro-climatically sensitive regions of the Indian Himalayas.

2 Study area and data

2.1 Study area: Jhelum River Basin

Jhelum River Basin is a significant trans-boundary sub-catchment within the Indus River system, with a total drainage area of approximately 33,867 km² (Figure 1).

Geographically the basin lies at 32°20′N to 35°20′N and longitudes from 73°40′E to 75°35′E. On the physiographic side, the basin is characterized by high-gradient topography with mountainous sections at the northern part of the basin characterized by high geography to a relatively low one at the southern part of the basin. This elevation-induced heterogeneity is very critical in the development of the hydrological response of the basin in seasons.

Climatologically, Jhelum River Basin is predominantly controlled by the Indian summer monsoon and mid-latitude Western Disturbances. In terms of precipitation the annual averages in the lower catchment reaching approximately 700 mm at Srinagar, and up to 1,499 mm at higher catchments like Garhi Dopatta, then decreasing back to around 853 mm at Mangla. Approximately 38% of the yearly rainfall falls in the winter season with the balance of 62% being received in the summer seasons. It is also notable that the Kashmir Valley has a characteristic precipitation pattern, the proportion of the annual precipitation falls higher during the winter and is lower during the summer, downstream areas receive relatively more precipitation. The basin is extremely reliant on the aspect of snow and glacier melts, especially the higher altitudes. The climatic warming has increased the rate of streamflow runoff unevenly and its effects have a direct bearing on water availability, and risk of floods. These changes in such fragile mountainous system are alarming for the need of assessing the future threats.

2.2 Data acquisition and pre-processing

The observed data for conducting this study was obtained from Indian Meteorological Department (IMD) during the years 1985–2014 (Table 1). Climatic data in form of GCMs were obtained from an earlier work conducted by Mishra et al. (2020). This dataset is bias-corrected climate projections at a spatial resolution 0.25°. The dataset contains 13 GCMs, where each model covers five scenarios (Historical, SSP126, SSP245, SSP370, and SSP585), with daily data available for precipitation, maximum temperature (Tmax), and minimum temperature (Tmin; Table 2). Bias correction of these fields was performed by Mishra et al. (2020) using Empirical Quantile Mapping (EQM) against IMD gridded observations over the historical reference period (1951–2014) at 0.25° spatial resolution, and applied to both the historical and projected (2015–2100) simulations for all 13 GCMs. EQM itself was not re-performed or re-parameterized in this study; the present contribution is the station-scale evaluation and multi-criteria ranking framework applied to this pre-corrected product. Also, the digital elevation model (DEM) was acquired from Shuttle Radar Topography Mission (SRTM) at a resolution of 30 m. The information employed in this study was collected from different sources as indicated in Table 1. In order to be consistent, the resampling of all datasets was performed to a shared spatial resolution of 0.25°, equivalent to this of data on climate.

Table 1. Data sources, spatial resolution, meteorological stations and corresponding CMIP6 GCM grid points used for model evaluation.

Table 2. List of GCMs, organisations and their spatial resolutions.

To equate the projections of the gridded model to the observed climate data, the bias-corrected GCM output was down-extracted to specific meteorological stations of the Jhelum Basin. To accomplish this, the station coordinates were mapped to the nearest locations within the grid in 0.25° × 0.25° dataset (Table 1). In case of stations between grid nodes, bilinear interpolation was used, so that the local variability would be better reproduced (Kim et al., 2019). As part of the preprocessing, the completeness of the observed precipitation records was assessed for all meteorological stations. All stations record 100% record completeness over the study period, and no gap-filling or imputation was required.

3 Methodology

3.1 Performance evaluation of GCMs

The performance evaluation was used to check the annual and seasonal data separately. To do seasonal analysis, the dataset on the was stratified by the seasons and the performance of each seasonal subset was considered separately, which also saved intra-annual climatic variability. In this study four major seasons were considered i.e., Winter (12, 1, 2), Spring (3, 4, 5), Summer–Monsoon (6, 7, 8, 9) and Autumn (10, 11). The degree of linear correlation between simulated and observed precipitation was measured by correlation coefficient (CC) and it is the capability of the model to reproduce the temporal variability and phase fidelity (Moriasi et al., 2007). Although CC does not measure magnitude errors, it gives fundamental information regarding whether models are appropriate in terms of whether they represent right dynamics of precipitation over time. To calculate the CC Equation 1 was implied.

CC =∑i=1n(Oi-O¯)(Si-S¯)∑i=1n(Oi-O¯)2∑i=1n(Si-S¯)2(1)

Further, Root mean square error (RMSE) and normalized version (NRMSE) were utilized to measure the total size of the deviation between simulations and observations (Khoshvaght et al., 2025; Willmott and Matsuura, 2005). RMSE is more focused on bigger error and will be sensitive to extreme inconsistencies and thus it is more applicable in the precipitation datasets that are of high variability (Equation 2). NRMSE also allows scale-free comparison among the stations and climatic regimes, which is essential in multi-station evaluations (Equation 3).

RMSE=1n∑i=1n(Si-Oi)2(2)

NRMSE =RMSEO¯(3)

To quantify sheer systematic overestimation or underestimation in simulated volumes of precipitation, the absolute normalized mean bias error (ANMBE, percent) was added (Equation 4). The bias measures play a critical role in precipitation researches, because climate models tend to have wet or dry systematically biases that has a direct effect upon hydrological and climate-impact measurements.

ANMBE(\%)=|∑i=1n(Si−Oi)∑i=1nOi|×100(4)

In order to measure predictive skill, a hydrological approach was embraced and the Nash Sutcliffe efficiency (NSE) was adopted (Equation 5). NSE compares observed values to the simulated values of reported precipitation against observed mean and it is confidently said to be a bench mark measure in hydrology because it is sensitive to both bias and variance (Gupta et al., 2009).

NSE=1−∑i=1n(Si−Oi)2∑i=1n(Oi−O¯)2(5)

To eliminate some of the limitations of NSE, the Kling-Gupta efficiency (KGE) was added, which breaks down the model performance into correlation, variability and biasing components (Equation 6). KGE offers a more diagnostic and balanced evaluation of precipitation simulations, so KGE is specifically well-suited to ranking climate models once they are biased (Gupta et al., 2009).

KGE=1-((r-1)2+(α-1)2+(β-1)2)(6)

where, r = CC, α=σSσO, and β=S¯O¯.

Finally, one of the measures was the probability of skill (PSS) applied in determining the probabilistic skill of the models, particularly the possibility of distinguishing various precipitation states. PSS can be used in the process of determining the reliability and utility of precipitation simulations in decision and risk-based processes (Jolliffe and Stephenson, 2012; Ferro and Stephenson, 2011). PSS is calculated using Equation 7 as shown below.

PSS = POD - POFD(7)

where,

POD=HH+M, POFD=FF+C

and H = hits, M = misses, F = false alarms and C = correct negatives.

3.2 GCMs ranking

In order to arrive at a holistic, unbiased, and well-constructed ranking of GCMs, this study will employ an integrated, multi-criteria decision-making (MCDM) framework agenda that incorporates Technique for Order Preference by Similarity to Ideal Solution (TOPSIS), Preference Ranking Organization Method for Enrichment Evaluation (PROMETHEE), and VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR). The combination of various MCDM methods is supported theoretically by the nature of the climate model assessment, when the performances of diverse models are different with respect to statistical measures, weather variables, seasons and stations, and no algorithm can be condensed to reflect model reliability in all cases (Palutikof et al., 2019). The use of several, conceptually independent MCDM procedures helps to decrease the methodological bias and to strengthen the confidence in the obtained rankings by the models which stand out as well in different, complementary decision philosophies (Baydaş et al., 2024).

3.2.1 Multicriteria decision making (TOPSIS, PROMPTEE and VIKOR)

The TOPSIS was used in determining the comparative proximity of each GCM toward an ideal state of performance (Hwang and Yoon, 1981). TOPSIS is based on the fact that the optimal alternative should simultaneously minimize the distance from a positive ideal solution. This represents the best performance that can be obtained under all criteria, and maximize the distance from a negative ideal solution, which is a measure of the poorest possible solution (Behzadian et al., 2012). First of all, a decision matrix was prepared where rows are GCMs (alternatives) and columns are standardized performance measures (criteria). The matrix was scaled using min–max normalization in order to eliminate the effects of scale among dissimilar metrics and map each criterion onto a common [0, 1] benefit scale. TOPSIS was calculated using Equation 8.

rij={xij−mini(xij)maxi(xij)−mini(xij), for benefit criteria, maxi(xij)−xijmaxi(xij)−mini(xij), for cost criteria. (8)

where (xij) refers to the value of criterion (j) of model (i). The criterion weights were then used to multiply the normalized values to get the weighted decision matrix (vij = wj rij). The ideal solutions of positive ((A+)) and negative ((A–)) were then established, using the criterion of benefit and cost. Euclidean distances from these reference solutions were calculated as Equation 9.

Di+=∑j(vij-vj+)2, Di-=∑j(vij-vj-)2(9)

Further, the Relative closeness coefficient (TOPSIS) is expressed as Equation 10.

Ci = Di-Di++Di-(10)

This served as the final TOPSIS ranking score, with higher values indicating superior overall performance.

To complement the distance-based assessment provided by TOPSIS, the PROMETHEE was utilized in order to explicitly analyze the pairwise dominance associations between GCMs. PROMETHEE compares the extent to which one model is better than another with regards to each criterion; hence giving understanding of competitive advantages and trade-offs that might be obscured in a summary score (Christian et al., 2016). For each pair of models a and b, preference functions were applied to convert performance differences into preference degrees. The aggregated preference index was computed as (Equation 11):

π(a,b)=∑j=1nwjPj(a,b)(11)

where Pj(a, b) is the preference function for criterion j. based on this index, the positive (outgoing) preference flow ϕ⁺(a) and negative (incoming) preference flow ϕ⁻(a) were computed as (Equation 12);

ϕ+(a)=1m−1∑b≠aπ(a,b), ϕ−(a)=1m−1∑b≠aπ(b,a)(12)

The Net Flow,

ϕ(a)=ϕ+(a)-ϕ-(a)

This indicates the dominance of a model as compared to all the others. The greater the net flow values, the greater the overall preference and performance. Also, PROMETHEE is especially useful to find models that reliably perform better than others when evaluated based on a number of criteria, even when overall aggregate rankings are tightly bundled around.

Finally, the VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR) approach was used to come up with a compromise solution which clearly trades the group utility and individual worst-case performance (regret) (Opricovic and Tzeng, 2004). In contrast to TOPSIS and PROMETHEE, VIKOR is proposed to solve the decision-making problems with conflicting criteria and identify solutions that provide the acceptable trade-offs. For each criterion, the best (fj) and worst (fj-) values were identified. The group utility measure (S_i) and individual regret measure (R_i*) were then computed as (Equation 13):

Si=Σjwjfj∗−fijfj∗−fj−, Ri=maxj[wjfj∗−fijfj∗−fj−](13)

These were collectively summed up into VIKOR index which is computed using Equation 14.

Qi=vSi-SS--S+(1-v)Ri-RR--R(14)

where v (set to 0.5) represents the decision-maker's preference for group utility vs. individual regret. Lower (Q_i) values indicate superior compromise solutions.

Combined use of TOPSIS, PROMETHEE, and VIKOR makes GCM ranking strong and credible both in terms of distance to ideal performance (TOPSIS), dominance and pairwise competitiveness (PROMETHEE), and compromise optimality over conflicting criteria (VIKOR). Those GCMs ranking high in all three processes can thus be deemed as statistically sound and physically sound and thus can be advised to be applied in station scale precipitation analysis, seasonal diagnostics and construction of ensembles of the models in the Jhelum Basin.

3.3 Integrated consensus ranking and seasonal and annual rank stability

After ranking the GCMs using three different methods, identification of an optimal GCM is an underdetermined decision problem. To reduce this uncertainty and prevent the methodological bias, the strategy of methodological pluralism was taken meaning that three conceptually different MCDM methods were combined: TOPSIS, PROMETHEE and VIKOR. A GCM model, which tracks uniformly on all three ranking paradigms is more stable to shifts in ranking logic such that it is a better choice to use in estimating precipitation at the station scale and building ensembles. This multi-method strategy ensures that the final selection is not sensitive to the choice of a specific MCDM algorithm but instead reflects comprehensive and defensible model skill.

In order to derive a robust overall ranking, the final ranks obtained independently from TOPSIS, PROMETHEE, and VIKOR were integrated using a consensus-based approach. GCM identifiers of each station were standardized and rankings were narrowed down to the results of station-specific interest; thus, methodological consistency of the methods was maintained. Let Ri(k) denote the rank assigned to GCM I by method k, where k = 1, 2, 3 correspond to TOPSIS, PROMETHEE, and VIKOR, respectively. The annual consensus mean rank for each GCM was computed as the arithmetic mean of ranks obtained from the three MCDM methods (Equation 15);

(R̄iannual=1KΣk=1KRi(k), K=3)(15)

To measure how much the ranking methods agreed on the ranking, the standard deviation of the ranks was used as a measure of the usual degree of stability of the rankings per year (Equation 16);

(σiannual=1K-1Σk=1K(Ri(k)-R̄iannual)2)(16)

GCMs were ranked in ascending order of R-iannual, with σiannual used as a secondary criterion to prioritize models exhibiting greater inter-method consistency.

In addition, to achieve intra-annual variability in model skill the consensus procedure was also employed at the seasonal level. Seasonal consensus mean rank of a specified season s was calculated as (Equation 17):

(R̄i,sseasonal=1KΣk=1KRi,s(k))(17)

The corresponding seasonal rank stability was calculated as (Equation 18):

(σi,sseasonal=1K-1Σk=1K(Ri,s(k)-R̄i,sseasonal)2)(18)

Within each season, GCMs were ordered in ascending order of R-i,sseasonal. To obtain a single representative seasonal indicator for each GCM, seasonal consensus ranks were averaged across all seasons (Equation 19):

(R̄iseason=1SΣs=1SR̄i,sseasonal)(19)

where S denotes the total number of seasons considered. Seasonal rank stability was similarly summarized by averaging seasonal standard deviations across all seasons.

Finally, an integrated annual-seasonal consensus ranking was developed to provide a holistic assessment of GCM performance by jointly considering long-term climatological accuracy and seasonal reliability. The integrated mean rank was computed by assigning equal importance to annual and seasonal consensus information (Equation 20):

(R̄iintegrated=Riannual¯+Riseason¯2)(20)

Ranking stability at the integrated scale was evaluated as the mean of annual and seasonal rank standard deviations (Equation 21):

(σiintegrated=σiannual+σiseason2)(21)

GCMs were then ordered in ascending order of R-iintegrated, with σiintegrated used to resolve ties and emphasize robustness. This integrated framework ensures that selected models exhibit strong interannual performance, consistent seasonal behavior, and stability across multiple ranking algorithms, thereby providing a scientifically robust basis for GCM selection and ensemble development in station-scale precipitation studies.

To interpret the resulting standard deviations consistently, GCMs were classified into three rank-stability classes: stable, where the standard deviation of a GCM's ranks did not exceed one rank unit; moderately stable, between one and two units; and unstable, above two units. The one-unit threshold corresponds to disagreement of at most one position among the three MCDM methods. The two-unit threshold represents approximately half the standard deviation expected under a null model of randomly permuted ranks over the 13 GCMs considered here (σ_null ≈ 3.74, computed from the variance of a discrete uniform permutation over 13 items). This classification was applied uniformly to inter-method (σ_method), inter-station (σ_station), inter-season (σ_season), and integrated stability values.

3.4 Inter-method agreement assessment

Before combining the three MCDM methods into a consensus ranking, the degree of agreement between them was quantified at each station-season combination using two complementary measures. Kendall's coefficient of concordance (W) was computed across the three rank vectors (TOPSIS, PROMETHEE, VIKOR) as (Equation 22):

W=12Sk2(n3-n)(22)

where S is the sum of squared deviations of the total ranks assigned to each GCM from their mean, k is the number of ranking methods (k = 3), and n is the number of GCMs (n = 13). W ranges from 0 (no agreement) to 1 (perfect agreement); values above 0.70 were treated as indicating strong concordance, following conventional interpretation thresholds for Kendall's W. Pairwise Spearman rank correlations (ρ) were additionally computed between each pair of methods (TOPSIS-PROMETHEE, TOPSIS-VIKOR, PROMETHEE-VIKOR) at every station-season combination, to characterize agreement at the level of individual method pairs rather than only across all three simultaneously.

3.5 Sensitivity and robustness analysis

To test the robustness of the integrated consensus ranking to methodological choices, sensitivity analyses were performed along four independent axes, following recommendations for reporting MCDM robustness in climate-model evaluation (Rahman and Pekkat, 2024).

• a Weighting scheme: Equal weights (w_*j* = 1/7 for each of the seven criteria: CC, RMSE, NRMSE, ANMBE, NSE, KGE, PSS) were used as the baseline, reflecting the absence of an accepted physical basis for privileging any one criterion in climate-model evaluation. Four alternative weighting schemes were tested: (i) an entropy-based scheme, in which weights were derived from the information entropy of each criterion's distribution across GCMs, assigning greater weight to criteria with higher discriminating variance; (ii) an efficiency-emphasis scheme, up-weighting the efficiency-family criteria (NSE, KGE, and CC) (iii) an error-emphasis scheme, up-weighting RMSE, NRMSE, and ANMBE; and (iv) a probabilistic-emphasis scheme, up-weighting PSS. Each alternative weighting was applied to the full TOPSIS-PROMETHEE-VIKOR–consensus pipeline, and the resulting GCM ranking was compared against the equal-weight baseline using Spearman rank correlation (ρ).
• b Normalization scheme. Because TOPSIS's Euclidean-distance calculation is sensitive to how the decision matrix is normalized, four normalization schemes were tested: min-max normalization (the baseline; Equation 8), vector normalization, sum normalization, and maximum normalization. PROMETHEE and VIKOR are unaffected by this choice, since their calculations are invariant to positive linear rescaling of the criteria.
• c Method parameters. For PROMETHEE, eight combinations of the indifference threshold (q, ranging 0–0.15) and preference threshold (p, ranging 0.20–0.50) were tested against the baseline thresholds. For VIKOR, the compromise parameter v which balances group utility against individual regret in Eq. 14 was varied across five values (v = 0, 0.25, 0.50, 0.75, 1.00) against the baseline of v = 0.5.
• d Monte-Carlo weight uncertainty. To characterize uncertainty in the ranking arising from the arbitrary choice of criterion weights, two Monte-Carlo ensembles of 2,000 weight vectors each were generated: (i) weights drawn uniformly from a Dirichlet distribution over the seven-dimensional weight simplex, and (ii) weights generated by perturbing the equal-weight baseline by up to ±50% per criterion, followed by renormalization to sum to unity. For each of the 4,000 total draws, the full ranking pipeline was recomputed, and the resulting rank distribution for each GCM was summarized as its median rank, 90% rank interval (5th−95th percentile), and probability of achieving rank 1 or a top-3 placement.

All sensitivity configurations, together with the Monte-Carlo ensembles, were compared against the equal-weight baseline ranking using Spearman rank correlation, with ρ ≥ 0.70 treated as indicating substantial agreement (Schober and Schwarte, 2018).

3.6 Post-correction distributional diagnostics

To evaluate the post-correction agreement between the bias-corrected CMIP6 output and station observations at the level of the full precipitation distribution beyond the summary metrics in Section 3.2, wet-day precipitation (defined as daily totals ≥ 1 mm) was compared between the consensus best-ranked and worst-ranked GCMs and IMD observations at three representative stations spanning the basin's elevation gradient (Gulmarg, high-altitude; Pahalgam, mid-elevation; Srinagar, valley). Two complementary diagnostics were produced for each station: (i) kernel density estimates (KDE, Gaussian kernel) of wet-day daily precipitation on a logarithmic intensity scale, and (ii) quantile–quantile (Q–Q) plots comparing the 1st through 99th percentile of the simulated wet-day distribution against the corresponding observed percentiles, with the 1:1 line shown as a reference for perfect agreement. Because raw, pre-correction CMIP6 output was not available to this study [bias correction having been performed independently by Mishra et al. (2020), these diagnostics characterize post-correction agreement only, and do not quantify the magnitude of the correction itself; this is stated as an explicit limitation in Section 7.5.

3.7 Seasonal wet-day frequency and intensity diagnostics

To provide empirical support for the physical mechanisms proposed to explain seasonal and elevation-dependent performance contrasts (Section 7.2), two station-level diagnostics were computed for each of the four seasons, at the same three representative stations used in Section 3.7: (i) seasonal wet-day frequency, defined as the count of days with precipitation ≥1 mm within each season across the full study period, computed separately for IMD observations and for the consensus best-ranked GCM; and (ii) seasonal mean wet-day intensity, defined as the mean daily precipitation on wet days within each season. Both diagnostics were computed identically for observed and simulated series to allow direct comparison of frequency bias and intensity bias by season and by station.

3.8 Spatial bias mapping

To visualize the spatial distribution of station-level bias across the basin, the annual ANMBE (Equation 4) for the consensus best-ranked GCM at each of the six stations was interpolated across the basin extent using inverse-distance weighting (IDW) (Equation 23):

Z(x,y)=∑iwiZi∑iwi, where wi=1dip(23)

where Z(x, y) is the interpolated ANMBE value at grid location (x, y), Z_**i** is the observed ANMBE at station i, di is the Euclidean distance from (x, y), to station i, and p is the power parameter (set to p = 2, the conventional default). Interpolation was performed on a regular grid (300 cells along the longer basin dimension) spanning the extent of the basin DEM (SRTM, 30 m resolution), and grid cells outside the basin outline were masked using the DEM's own valid-data footprint, identified via the DEM's nodata value and each grid cell's corresponding raster row/column computed from the DEM's affine transform. Because interpolation is based on only six station values, values between stations should be interpreted as indicative of spatial trend rather than as independently verified estimates. The overall methodological framework adopted in this study is illustrated in Figure 2.

4 Results

4.1 Seasonal station-wise performance evaluation of CMIP6 GCMs

Performance at this high-altitude station varies strongly by season. In Gulmarg during autumn, correlations are weak (CC ≈ −0.03 to 0.06) and NSE and KGE are negative for every model (NSE −2.34 to −0.42; KGE −0.50 to 0.05). ACCESS-ESM1-5, INM-CM5-0, and ACCESS-CM2 give the lowest RMSE (7.36–8.7 mm) and the least-negative NSE. MRI-ESM2-0, EC-Earth3-Veg, and CanESM5 perform worst; CanESM5 alone shows an extreme bias (ANMBE ≈ 105.8%). In spring, RMSE narrows to 13.27–15.42 mm. BCC-CSM2-MR, CanESM5, and NorESM2-LM give the lowest RMSE and highest NSE, while MPI-ESM1-2-HR and MPI-ESM1-2-LR are consistently weakest. Summer monsoon is the station's best season overall: CanESM5, MPI-ESM1-2-LR, and MPI-ESM1-2-HR post the lowest RMSE (8.47–9.7 mm) and least-negative NSE. BCC-CSM2-MR has the highest RMSE and the most negative NSE (−1.49) here, despite the best PSS (0.964). In winter, BCC-CSM2-MR, INM-CM4-8, and ACCESS-ESM1-5 give the best deterministic scores; ACCESS-CM2 gives the best PSS (0.932). MRI-ESM2-0, NorESM2-LM, and INM-CM5-0 are weakest. Efficiency metrics stay negative for every model in every season, consistent with structural limits in reproducing snow-regime precipitation. Seasonal multi-metric performance is shown in Figure 3.

The valley station Srinagar shows weak performance across all seasons. In autumn, correlations sit near zero and NSE/KGE are strongly negative (NSE down to −3.60). ACCESS-ESM1-5, INM-CM5-0, and ACCESS-CM2 perform best (RMSE 5.47–7.4 mm). CanESM5 shows the most extreme overestimation in the whole dataset (ANMBE ≈ 232%); MRI-ESM2-0 and EC-Earth3-Veg are also weak. Spring improves modestly (RMSE 8.63–11.72 mm), led by BCC-CSM2-MR, CanESM5, and NorESM2-MM; MPI-ESM1-2-HR, MPI-ESM1-2-LR, and INM-CM5-0 are weakest. In summer monsoon, CanESM5, MPI-ESM1-2-HR, and MPI-ESM1-2-LR give the lowest RMSE and best NSE; ACCESS-ESM1-5 shows the highest bias (ANMBE ≈ 97.5%). Winter is this station's best season: ACCESS-ESM1-5, BCC-CSM2-MR, and EC-Earth3 give the lowest RMSE and highest NSE. NSE and KGE stay negative for every model in every season here.

In Pahalgam performance sits between Srinagar and Gulmarg. In autumn, ACCESS-ESM1-5, INM-CM5-0, and ACCESS-CM2 give the lowest RMSE (6.57–8.0 mm) and best NSE; CanESM5, MRI-ESM2-0, and EC-Earth3-Veg are weakest, with CanESM5 reaching ANMBE ≈ 179%. Spring improves (RMSE 12.68–16.33 mm), led by BCC-CSM2-MR, CanESM5, and NorESM2-LM, though NSE stays negative throughout. Summer monsoon is the strongest season here: MPI-ESM1-2-HR, CanESM5, and MPI-ESM1-2-LR give the lowest RMSE and best NSE, and NorESM2-MM gives the best PSS (0.946). In winter, ACCESS-ESM1-5, BCC-CSM2-MR, and EC-Earth3 remain the top performers, with PSS reaching 0.978 (MPI-ESM1-2-HR). At Kupwara seasonal behavior here resembles the other mid-altitude stations. In autumn, ACCESS-ESM1-5, INM-CM5-0, and BCC-CSM2-MR give the lowest RMSE and best NSE; CanES5's ANMBE reaches ≈ 180%. Spring improves moderately (RMSE 11.80–14.91 mm), led by BCC-CSM2-MR, NorESM2-LM, and CanESM5. Summer monsoon is comparatively favorable: CanESM5, MPI-ESM1-2-LR, and INM-CM5-0 give the lowest RMSE and best NSE. In winter, ACCESS-ESM1-5, BCC-CSM2-MR, and EC-Earth3 perform best. Bias stays comparatively low here mostly under 10% for the best models even though efficiency metrics remain negative throughout.

At Qazigund performance is strongly season-dependent and broadly resembles Srinagar. Autumn is the poorest season, with NSE as low as −3.08 and CanES5's ANMBE reaching ≈ 213%. ACCESS-ESM1-5, INM-CM5-0, and ACCESS-CM2 perform best. Spring shows lower RMSE and bias for BCC-CSM2-MR, CanESM5, and NorESM2-MM. Summer monsoon improves for CanESM5, MPI-ESM1-2-HR, and MRI-ESM2-0. Winter is this station's best season: ACCESS-ESM1-5, BCC-CSM2-MR, and EC-Earth3 lead. Deterministic performance stays weak for every model in every season. At Kokernag performance sits between the high-altitude and valley stations. Autumn is by far the weakest season, with NSE down to −4.19 and CanESM5 reaching the most extreme bias in the entire dataset (ANMBE ≈ 282%). ACCESS-ESM1-5, INM-CM5-0, and ACCESS-CM2 perform best. Spring favors BCC-CSM2-MR, CanESM5, and NorESM2-MM. Summer monsoon is the most favorable season here: CanESM5, MPI-ESM1-2-HR, and MRI-ESM2-0 improve on RMSE and NSE. In winter, ACCESS-ESM1-5, BCC-CSM2-MR, and EC-Earth3 lead, with EC-Earth3 reaching the highest PSS in the seasonal dataset (0.986). NSE and KGE stay negative for every model.

4.2 Annual station-wise GCM performance evaluation

At Gulmarg the annual CC stays weak overall (−0.018 to 0.031). RMSE ranges from 12.23 to 13.23 mm; CanESM5, INM-CM5-0, and NorESM2-MM give the lowest errors. NSE stays negative for every model (−0.666 to −0.423); KGE ranges from −0.121 to −0.052. ANMBE ranges from 26.7% (MRI-ESM2-0) to 35.2% (ACCESS-CM2). Annual multi-metric performance is shown in Figure 4. At Srinagar Annual CC is weak (−0.018 to 0.031). RMSE ranges from 9.13 to 10.25 mm the lowest of any station, though this reflects Srinagar's shorter record and lower absolute precipitation rather than better skill. NSE is strongly negative (−1.129 to −0.688); KGE ranges from −0.053 to 0.004. ANMBE ranges from 13.3% (CanESM5) to 26.8% (MRI-ESM2-0). Srinagar records the weakest annual efficiency scores among all stations (Figure 4).

At Pahalgam the annual CC ranges from −0.017 to 0.036. RMSE ranges from 12.50 to 13.66 mm; CanESM5, NorESM2-LM, and BCC-CSM2-MR give the lowest errors. NSE stays negative (−0.829 to −0.530); KGE ranges from −0.053 to 0.003. ANMBE ranges from 14.1% to 22.6% among the narrowest bias spreads of the six stations. Pahalgam shows the best annual PSS values in the dataset (0.968–0.978). At Kupwara, Annual CC ranges from −0.002 to 0.037. RMSE spans 10.80–11.64 mm; EC-Earth3, CanESM5, and ACCESS-ESM1-5 give the lowest errors. NSE stays negative (−1.110 to −0.816); KGE ranges from −0.007 to 0.035, the least-negative range of any station. ANMBE ranges from just 1.5% to 8.2% the lowest annual bias of all six stations. No single GCM performs consistently well across all metrics here.

At Qazigund the Annual CC ranges from −0.018 to 0.043. RMSE is among the highest across stations (13.49–14.25 mm); CanESM5, INM-CM5-0, and ACCESS-CM2 give the lowest errors. NSE ranges from −0.663 to −0.492; KGE from −0.059 to 0.012. ANMBE ranges from 5.8% (NorESM2-MM) to 15.4% (CanESM5) far below the extreme autumn bias at this station, showing how annual averaging masks acute seasonal overestimation. At Kokernag, Annual performance sits between the high-altitude and valley stations. CC ranges from −0.019 to 0.027. RMSE ranges from 12.24 to 13.02 mm; CanESM5, ACCESS-CM2, and NorESM2-LM give the lowest errors. NSE spans −1.027 to −0.791; KGE ranges from −0.027 to 0.019. ANMBE ranges from 4.3% (CanESM5) to 16.0% (NorESM2-MM) again far below the autumn extreme, where CanESM5 reaches nearly 282% bias at this station.

4.3 Post-correction distributional diagnostics

To support the magnitude-agreement statistics, the bias-corrected GCM precipitation distribution was compared against IMD observations at three representative stations spanning the basin's elevation range. We compare wet-day (≥1 mm) precipitation using kernel density estimates on a log scale, and quantile–quantile plots against the 1st−99th observed percentiles, for the best- and worst-ranked models in the integrated consensus ranking. Raw, pre-correction CMIP6 output was not available to this study, since bias correction was performed independently by Mishra et al. (2020). These diagnostics therefore show post-correction agreement, not the size of the correction itself.

Both models track the observed distribution closely at low-to-moderate intensities, at all three stations. Above roughly 20 mm/day, both models, including MPI-ESM1-2-LR, the top-ranked GCM fall below the 1:1 line, showing an under-prediction of high-intensity wet-day precipitation. This under-prediction is largest at Gulmarg, the highest station, where the gap opens up beyond roughly the 90th percentile. MRI-ESM2-0 and MPI-ESM1-2-LR track almost identically in this upper tail at all three stations (Figure 5). This residual tail under-prediction is consistent with prior bias-correction work in the Jhelum basin: linear-scaling-corrected CMIP GCM precipitation has similarly under-represented high-intensity events despite good agreement in the bulk of the distribution (Mahmood et al., 2018), and a multi-method bias-correction evaluation for the basin likewise found that correction reduces but does not eliminate underestimation of precipitation extremes, particularly over its higher-elevation reaches (Ansari et al., 2023). This suggests the tail bias is a structural limitation of statistical bias correction in this basin rather than one specific to EQM.

5 Ranking of GCMs

5.1 Seasonal and annual station-wise GCM ranking using TOPSIS

Autumn shows the clearest and most consistent model hierarchy in the basin. ACCESS-ESM1-5 gives the highest closeness coefficient at five of six stations: Qazigund (0.803), Kupwara (0.821), Pahalgam (0.779), Kokernag (0.819), and Gulmarg (0.778). Only at Srinagar does it fall to second (0.750), where INM-CM5-0 leads (0.774). INM-CM5-0 and MPI-ESM1-2-LR are consistently the closest competitors. Spring shows a more station-specific pattern. CanESM5 leads at Kokernag (0.698), Kupwara (0.590, second to NorESM2-LM), Pahalgam (0.666, second to EC-Earth3-Veg), and Srinagar (0.733); NorESM2-MM leads at Gulmarg (0.670). Leading scores in spring range from 0.59 to 0.73 narrower than in autumn, showing weaker discrimination among models in this transitional season. In the summer monsoon, MPI-ESM1-2-LR dominates: it leads at Kokernag (0.744), Pahalgam (0.711, second to INM-CM48's 0.850), Qazigund (0.804), and Srinagar (0.853), and places second at Gulmarg and Kupwara. Winter shows the most station-dependent pattern of any season: NorESM2-MM leads at Gulmarg and Qazigund, INM-CM4-8 at Kokernag, Pahalgam, and Srinagar, and EC-Earth3 at Kupwara. No GCM leads in more than one season at any station. Only ACCESS-ESM1-5 in autumn and MPI-ESM1-2-LR in the summer monsoon show consistent basin-wide dominance within a single season. Seasonal and station-wise TOPSIS scores are shown in Figure 6.

At the annual scale, TOPSIS again shows strong station dependence, but with tighter score margins than the earlier, uncorrected analysis reflecting the narrower true separation between models once metrics are computed on correctly aligned daily series. NorESM2-MM leads at Gulmarg (0.877), ahead of INM-CM4-8 (0.853) and NorESM2-LM (0.827) the closest three-way margin of any station. ACCESS-CM2 leads at Kokernag (0.805). CanESM5 leads at both Kupwara (0.733) and Srinagar (0.781). NorESM2-LM leads clearly at Pahalgam (0.883). INM-CM5-0 leads at Qazigund (0.773). No GCM leads at more than one station annually confirming that the annual ranking, like the seasonal one, is dominated by local station character rather than by any single model's overall skill. Annual TOPSIS scores are shown in Figure 6.

5.2 Seasonal and annual station-wise GCM ranking using PROMETHEE

PROMETHEE-II broadly confirms the TOPSIS hierarchy. In autumn, ACCESS-ESM1-5 records the highest net outranking flow at Kokernag (Φ = 0.6892, the worked example in Section 6.1) and leads consistently across the basin, closely followed by INM-CM5-0 and ACCESS-CM2. MPI-ESM1-2-LR and NorESM2-MM form a stable second tier. CanESM5, EC-Earth3, and MRI-ESM2-0 record the weakest autumn flows basin-wide; CanESM5 records the single lowest net flow anywhere in the study (Φ = −0.8488 at Kokernag). The wide gap between the leading and trailing groups confirms autumn as the season with the strongest model discrimination. Seasonal PROMETHEE ranks are shown in Figure 7.

Annual PROMETHEE rankings track the TOPSIS results closely. NorESM2-MM leads at Gulmarg (Φ = 0.4554); CanESM5 leads at Kokernag (0.5205, narrowly ahead of ACCESS-C2's 0.5117), Kupwara (0.6027), and Srinagar (0.6089); NorESM2-LM leads at Pahalgam by a wide margin (0.7052, more than double the second-placed INM-CM5-0); INM-CM5-0 leads at Qazigund (0.5422). The close margin at Kokernag (CanESM5 and ACCESS-CM2 separated by only 0.0088) foreshadows the VIKOR compromise-set result for that station (Section 5.3). Annual PROMETHEE ranks are shown in Figure 7.

5.3 Seasonal and annual station-wise GCM ranking using VIKOR

VIKOR results for autumn agree closely with TOPSIS and PROMETHEE. Because Q is min–max normalized within each comparison set, the leading model at any station or season lands at or near Q = 0. The index expresses relative, not absolute, compromise quality, and should be read together with the acceptability conditions below. Following (Opricovic and Tzeng, 2004), we tested two acceptability conditions at all 30 station–season combinations. The acceptable-advantage condition requiring the Q-gap between the first- and second-ranked model to exceed DQ = 1/(m – 1) = 0.083 for 13 alternatives holds in 19 of 30 cases. In the other 11, a compromise set of two or three models applies instead of a single winner. This includes the Gulmarg (advantage 0.042), Srinagar (0.020), and Kupwara (0.042) annual comparisons, and the Kokernag autumn comparison (0.018) used as the worked example in Section 6.1. The acceptable-stability condition holds in 29 of 30 cases, failing only at Kupwara in spring. This roughly one-in-three failure rate of the acceptable-advantage condition is itself a finding it is VIKOR's own confirmation of the narrow inter-model separation discussed at the consensus level in Section 6.4. Seasonal VIKOR rankings are shown in Figure 8.

Annual VIKOR rankings closely track PROMETHEE and TOPSIS. NorESM2-MM is the best compromise solution at Gulmarg, but its advantage over INM-CM4-8 and NorESM2-LM (0.042) falls below the acceptability threshold, so all three form a compromise set there. ACCESS-CM2 leads clearly at Kokernag (Q = 0.006); CanESM5 leads at Kupwara (Q = 0.000, though it fails the advantage test against MPI-ESM1-2-LR); NorESM2-LM leads decisively at Pahalgam (Q = 0.000); INM-CM5-0 leads at Qazigund (Q = 0.000). At Srinagar, CanESM5, BCC-CSM2-MR, and ACCESS-CM2 form a three-way compromise set. Four of the six stations fail the acceptable-advantage condition at the annual scale. Annual VIKOR rankings are shown in Figure 8.

6 Rank stability and integrated consensus ranking

6.1 Inter-method agreement

Before combining the three methods, we checked how well they agree. Kendall's coefficient of concordance across the three rank vectors averages W = 0.922 over the 30 station–season combinations, ranging from 0.752 at Pahalgam (winter) to 0.988 at Kupwara (winter; Figure 9). Pairwise Spearman correlations average 0.941 between TOPSIS and PROMETHEE, 0.897 between TOPSIS and VIKOR, and 0.812 between PROMETHEE and VIKOR. All three exceed the conventional 0.70 threshold for strong concordance at every one of the 30 combinations, confirming the three methods identify substantially the same hierarchy despite different underlying logic.

Table 3 gives the full worked calculation for one station and season, chosen because it also shows a VIKOR compromise-set case (Section 5.3). At Kokernag in autumn, ACCESS-ESM1-5 ranks first by all three methods a mean rank of 1.00 and zero standard deviation, the stable class. INM-CM5-0 receives ranks 2, 2, and 3 (mean 2.33, SD 0.58). MPI-ESM1-2-LR shows the difference between rank magnitude and rank stability: it ranks 4th by TOPSIS and PROMETHEE but 2nd by VIKOR (mean 3.33, SD 1.15, moderately stable), reflecting VIKOR's 0.018 advantage over ACCESS-ESM1-5 falling below the acceptability threshold.

Table 3. Worked example of rank-stability calculation Kokernag, autumn.

ACCESS-ESM15's VIKOR advantage over the second-ranked model here is 0.018 below the acceptable-advantage threshold (DQ = 0.083 for 13 alternatives). ACCESS-ESM1-5 and MPI-ESM1-2-LR therefore form a VIKOR compromise set at this station–season, consistent with Section 5.3.

For Rank-stability classification we classify a GCM as stable where its rank standard deviation stays at or below one rank unit, moderately stable between one and two units, and unstable above two units. The one-unit threshold means disagreement of at most one position among methods. The two-unit threshold is roughly half the standard deviation expected from randomly permuted ranks over 13 alternatives (3.74).

6.2 Seasonal rank stability

We evaluated rank stability in three dimensions: across the three MCDM methods, across the six stations, and across the four seasons. Table 4 presents seasonal consensus ranks, averaged over all stations, for all four seasons together.

Table 4. Seasonal rank stability of CMIP6 GCMs across the Jhelum River Basin.

In autumn, ACCESS-ESM1-5 ranks first (mean rank 1.50) with the lowest inter-method variability among the leaders (σ method = 0.58), followed by INM-CM5-0 (2.44). Both are also spatially consistent (σ station 0.75 and 0.58). MPI-ESM1-2-LR (4.00) and ACCESS-CM2 (4.39) follow. CanESM5 is unambiguously the weakest performer, ranking last at every station (mean 12.83, σ method = 0.00); MRI-ESM2-0 is nearly as weak (12.06). In spring, CanESM5 leads (2.06), followed by EC-Earth3-Veg (3.33) and NorESM2-MM (4.89), the latter showing the highest spatial variability among the leaders (σ station = 3.26). BCC-CSM2-MR shows the highest inter-method disagreement of any model in any season (σ method = 2.87). MPI-ESM1-2-HR is the weakest spring performer (12.39) with very low variability, showing a consistent rather than occasional weakness. In the summer monsoon, MPI-ESM1-2-LR leads (1.50) with low inter-method variability (σ method = 0.36), followed by INM-CM5-0 (3.22) and ACCESS-CM2 (3.94). ACCESS-ESM1-5 is weakest here (12.28) a sharp reversal from its first-place standing in autumn, and the clearest single illustration in this study of how much model skill depends on season. In winter, INM-CM4-8 leads (2.22), followed by NorESM2-MM (2.89). EC-Earth3-Veg shows the season's highest spatial variability (σ station = 3.13), consistent with the elevation dependence of winter snowfall processes discussed throughout this study. MRI-ESM2-0 is again weakest and most consistently so (12.89).

No model ranks first in more than one season: ACCESS-ESM1-5 in autumn, CanESM5 in spring, MPI-ESM1-2-LR in the summer monsoon, INM-CM4-8 in winter. The autumn leader (ACCESS-ESM1-5) is the summer-monsoon's worst performer. The winter leader (INM-CM4-8) ranks only 9th in autumn. This is the strongest evidence in this study that model selection needs seasonal, not just annual, evaluation. Figure 10a shows the basin-averaged version of this reordering; Figure 10b shows the same pattern at the level of individual stations.

6.3 Annual rank stability

At the annual scale, MPI-ESM1-2-LR ranks first with a mean rank of 3.78 across the six stations and the lowest inter-method variability of any leading model (σ method = 0.44), followed by CanESM5 (3.94) and INM-CM5-0 (5.00; Table 5). CanESM5 shows much higher spatial variability than MPI-ESM1-2-LR (σ station 3.34 vs. 1.46), reflecting its strong showing at Kupwara and Srinagar against weaker standing elsewhere (Figure 11). ACCESS-CM2 shows the largest spatial spread of any model in the annual ranking (σ station = 4.41). MRI-ESM2-0 is unambiguously last (12.89) with minimal variability in either dimension a consistently, not occasionally, weak model.

Table 5. Annual rank stability of CMIP6 GCMs over the JRB.

Table 6 summarizes the three best- and three worst-performing GCMs at each station, with the metrics the reviewer requested. At Gulmarg, the top three (NorESM2-MM, INM-CM4-8, NorESM2-LM) are separated by consensus rank alone, with zero inter-method disagreement. At Qazigund, the third-place model (CanESM5) shows an inter-method SD of 1.73, meaning a genuine three-way split rather than a clear third place.

Table 6. Top three and bottom three GCMs at each station by integrated consensus rank.

6.4 Integrated consensus ranking

The integrated annual–seasonal consensus identifies MPI-ESM1-2-LR as the best-performing GCM overall, with an annual mean rank of 3.78, a seasonal mean rank of 4.53, and an integrated mean rank of 4.15 (Table 7). Its integrated rank stability (1.57) comes almost entirely from low inter-method disagreement (0.64); its seasonal variability (2.31) sits close to the basin median. This combination a low mean rank achieved without extreme variability in any single dimension is why we identify it as the most defensible single choice for station-scale precipitation analysis and ensemble building in this basin.

Table 7. Integrated annual/seasonal consensus ranking of CMIP6 GCMs over the JRB.

Integrated stability is the mean of the three variance components (method, station, season).

At the opposite end, MRI-ESM2-0 is unambiguously the weakest model, with an integrated mean rank of 11.81 and the lowest inter-method variability of any model in the study (0.41). Its poor standing is not an artifact of disagreement between methods all three agree. These two results MPI-ESM1-2-LR best, MRI-ESM2-0 worst are the only positions in the integrated ranking we treat as firm conclusions. Section 6.5 shows both are robust to every weighting scheme, normalization choice, and method parameter tested, and that MPI-ESM1-2-LR is the Monte-Carlo modal top-ranked model under 67.5% of randomly drawn weight vectors.

The 11 positions between rank 2 and rank 12 should not be read as a strict ordering. The 90% weight-uncertainty interval for INM-CM5-0 (rank 2) spans ranks 1-5; for CanESM5 (rank 6) it spans 1–12 effectively the whole table (Figure 12). Seven models INM-CM5-0, INM-CM4-8, ACCESS-CM2, NorESM2-MM, CanESM5, NorESM2-LM, and ACCESS-ESM1-5 have overlapping 90% intervals covering ranks 3 through 10. The data cannot tell us which of these is genuinely 3rd-best and which is genuinely 10th-best. This follows directly from the narrow inter-model separation confirmed in Section 4: once metrics are computed on correctly aligned daily data, the 13 GCMs differ from each other by about as much as the uncertainty introduced by the choice of criterion weights.

Breaking rank stability into its three components (Figure 13) shows that seasonal variability (0.96–4.62) generally exceeds both inter-method variability (0.41–1.24) and spatial variability (0.96–3.48). Rank instability in this basin is mostly seasonal, not methodological. CanESM5 and ACCESS-ESM1-5 show the most extreme seasonal variability (4.62 and 4.57); both swing from a top-two seasonal ranking in one season to bottom-two in another (Section 6.2) a swing the integrated mean rank alone does not show.

6.5 Sensitivity of the ranking to weighting and method parameters

We tested the robustness of the integrated consensus ranking along four independent axes weighting scheme, TOPSIS normalization, PROMETHEE preference thresholds, and the VIKOR compromise parameter together with two Monte-Carlo ensembles over the weight vector.

For Weighting equal weights were adopted as the baseline: 1/7 (≈0.143) for each of the seven criteria (CC, RMSE, NRMSE, ANMBE, NSE, KGE, PSS), in the absence of an accepted physical basis for favoring any one metric over another in climate-model evaluation (exact weights in Table 8). This avoids reintroducing the analyst subjectivity a multi-criteria framework is meant to remove. We tested four alternative schemes (Figure 14). Entropy-based weights (ρ = 0.9286) and an efficiency-emphasis scheme (ρ = 0.984) both keep MPI-ESM1-2-LR as the top-ranked model. An error-emphasis scheme, which up-weights RMSE, NRMSE, and ANMBE, departs more (ρ = 0.786) but still ranks MPI-ESM1-2-LR first. Only a probabilistic-emphasis scheme, which weights PSS heavily, changes the leader (ρ = 0.593; ACCESS-CM2 first, MPI-ESM1-2-LR second). This is the only weighting scheme tested that displaces MPI-ESM1-2-LR from first place.

Table 8. Criterion weights under the equal-weight baseline used in the integrated consensus ranking.

For Normalization, TOPSIS under vector, min–max, sum, and maximum normalization were computed. Because the vector, sum, and maximum schemes are simple linear rescalings of the raw criteria, PROMETHEE and VIKOR are unaffected by this choice; only TOPSIS, which uses absolute Euclidean distances, depends on normalisation. All four schemes give ρ ≥ 0.934 against the baseline and keep MPI-ESM1-2-LR as the leading model. In Method parameters, eight combinations of the PROMETHEE indifference and preference thresholds (q from 0 to 0.15, p from 0.20 to 0.50) all give ρ ≥ 0.984, with MPI-ESM1-2-LR first every time. The VIKOR compromise parameter varied from v = 0 to v = 1; all five settings give ρ ≥ 0.967 and keep MPI-ESM1-2-LR first. Neither parameter matters to the headline result.

For Monte-Carlo weight uncertainty, two ensembles of 2,000 draws each were generated: a uniform Dirichlet distribution over the weight simplex, and the equal weights perturbed by up to ±50%. Under uniform sampling, MPI-ESM1-2-LR ranks 1st in 67.5% of draws and reaches the top three in 86.6%, with a 90% rank interval of 1 to 6 (Figure 15). Under the constrained ±50% perturbation, it reaches the top three in 100.0% of draws. INM-CM5-0 is the only other model to exceed 70% top-3 probability under either scheme (73.3% uniform, 98.8% constrained). No model ranked 8th or below in the baseline ever reaches the top three under either scheme.

Across all configurations tested, 21 of 22 give ρ ≥ 0.786 against the equal-weight baseline, and 19 of 22 give ρ ≥ 0.929 (Figure 16 and Tables 9, 10). The consensus ranking's top position is robust to essentially every methodological choice we tested, with the single exception of an extreme probabilistic-emphasis weighting. This robustness at the extremes, set against the fragility of the middle ranks under the exact same tests (Section 6.4), is the central finding of the sensitivity analysis: a multi-criteria ranking can deliver a highly robust best/worst identification and a genuinely indeterminate middle order at the same time—and the two should be reported differently.

Table 9. Spearman rank correlation (ρ) against the equal-weight baseline across all 22 tested configurations.

Table 10. Monte-Carlo summary: probability of rank-1 and top-3 placement, with 90% rank-uncertainty intervals, under uniform (MC) and constrained ±50% (CP) weight sampling.

7 Discussion

7.1 Overall interpretation of CMIP6 precipitation skill in the Jhelum Basin

The station-scale evaluation shows that bias correction has cut, but not eliminated, systematic error in CMIP6 precipitation over the basin. Basin-wide, annual RMSE ranges from 9.13 to 14.25 mm and NRMSE from 3.16 to 5.05 across all 13 models and six stations. While ANMBE ranges from 1.5% to 35.2%, and consistent with EQM meaningfully reducing mean bias at the station scale. Time-dependent and efficiency-based indicators stay weak throughout: annual CC lies between −0.019 and 0.043, NSE between −1.129 and −0.423, and KGE between −0.121 and 0.035. No station or model combination in the dataset reaches a positive NSE.

This split real improvement in magnitude statistics alongside weak correlation and efficiency matches the established view that bias correction adjusts the statistical distribution of model output without improving its ability to reproduce day-to-day sequencing (Ehret et al., 2012; Maraun, 2016). CMIP6 historical simulations run freely rather than being nudged to observed weather, so a simulated day is not expected to line up with the matching observed day (Eyring et al., 2016). Near-zero daily correlation and negative daily NSE are therefore expected for any free-running GCM at daily resolution not a defect specific to this basin or this bias-correction method. This is consistent with similarly negative daily-scale NSE values reported for CMIP6 precipitation elsewhere in South Asia against the same IMD gridded observational network (e.g., NSE as low as −0.48 at daily and −0.78) at monthly resolution (Vinod and Agilan, 2024). This qualification applies to the daily-scale metrics (CC, NSE, KGE), not to the magnitude-based metrics (RMSE, NRMSE, ANMBE, PSS), which keep their ordinary meaning regardless of day-to-day alignment.

The elevation-dependent signal in the data is real, but it isn't a single clean gradient. NRMSE improves with elevation (r = −0.85 across the six stations), and so does NSE (r = +0.67). This fits with Gulmarg and Pahalgam sitting under the more organised, large-scale precipitation from Western Disturbances systems (Dimri et al., 2015) CMI6's native resolution captures better than the convective, orographically shielded processes that dominate lower elevations (Prein et al., 2015). But ANMBE and KGE move the opposite way with elevation (r = +0.74 and r = −0.83): Gulmarg, the highest station, has both the largest annual bias in the basin (30.3%) and the most negative KGE (−0.076), while Kupwara, one of the lowest, has the smallest bias (4.8%) and the least negative KGE (+0.011; Figure 17).

The spatial pattern of station-level bias (Figure 17) confirms this uneven elevation dependence directly. The highest annual ANMBE in the basin is at Gulmarg (30.3%), the highest station; the lowest is at Kupwara (4.8%). Critically, this relationship isn't monotonic: Srinagar, the lowest-elevation station in the basin, shows the second-highest bias (21.9%) higher than both Pahalgam (16.0%) and Kokernag (11.8%), both considerably higher stations.

Part of the apparent elevation effect on RMSE is a magnitude artifact, not a real skill difference. Srinagar's low absolute RMSE (9.13–10.25 mm) reflects a mean daily precipitation roughly half that of the higher stations, so its NRMSE (4.67, the basin's second-highest) is the fairer comparison. Taken together, elevation shifts which metrics improve and which get worse it doesn't give a uniform high-altitude advantage. Using RMSE alone to argue for better high-altitude performance in this basin would be misleading without the NRMSE and ANMBE context alongside it.

7.2 Seasonal contrasts and process-level insights

Seasonal analysis shows large, physically explainable contrasts in model performance. Autumn is the weakest season on four of five metrics: basin-mean NRMSE (7.36, roughly double any other season), ANMBE (82.7%), NSE (−2.00), and KGE (−0.45) are all worst in autumn. This fits with autumn being a transitional period between monsoon withdrawal and winter's onset, when precipitation is more episodic, more spatially patchy, and more driven by local-scale forcing that coarse-resolution GCMs are structurally poor at capturing (Prein et al., 2015). Two further mechanisms likely compound this transitional-season weakness. Orographic shadowing where a station sits in the lee of a ridge relative to the basin's dominant storm tracks can locally suppress or enhance precipitation in ways a 0.25° grid cell cannot resolve, and this effect is plausibly strongest precisely when the basin's dominant synoptic control is itself in transition, as in autumn (Roe, 2005). The basin's comparatively sparse six-station observation network further limits the ability to disentangle these mechanisms empirically; a denser network would be needed to separate the convective, orographic, and sampling contributions to the autumn signal with confidence. Correlation is the one metric where autumn isn't the weakest; basin-mean CC is instead lowest in the summer monsoon, though barely distinguishable from autumn and spring (all four seasonal means sit within ±0.005 of zero) underlining that daily correlation carries little useful information in any season, for the reasons given in Section 7.1. Winter shows the best probabilistic skill of any season (basin-mean PSS = 0.922) and comparatively good NRMSE and KGE, consistent with Western Disturbance-driven and snowfall precipitation being better represented at CMI6's resolution (Dimri et al., 2015) than the convective processes that dominate other seasons (Prein et al., 2015). Even so, efficiency metrics stay negative in winter for every station and model a season's relative improvement over autumn doesn't mean real skill in reproducing observed variability. The scale of the seasonal reordering is large: the leading model changes completely from one season to the next (Section 6.2), and models that lead one season can fall to the bottom of the basin-wide ranking in another. ACCESS-ESM1-5 is the clearest case first basin-wide in autumn, last in the summer monsoon a complete reversal that annual-only evaluation would hide entirely.

To back up the physical mechanisms behind these regional and seasonal contrasts, we evaluated seasonal wet-day frequency and mean-intensity diagnostics across the elevation gradient (Figure 18). High-altitude stations, like Gulmarg, show a small wet-day frequency bias (+7%), with their largest relative intensity under-prediction in autumn. This fits with winter and spring precipitation there being driven by widespread, stratiform Western Disturbance systems, which scale predictably across GCM grid cells, while autumn's post-monsoon transition brings localised convective variability that's harder to resolve. Srinagar, the valley station, shows a much larger wet-day frequency bias (+66%, vs. +7–7.5% at the two higher-elevation stations), and this holds across all four seasons. This fits with valley precipitation being more locally convective and topographically confined something a 0.25° grid cell can't resolve at the sub-grid scale, so the model marks days as wet when the station itself stayed largely dry (Prein et al., 2015). This frequency bias, not a seasonal one, is the main elevation-linked gap in the ranking framework's underlying error metrics.

We also mapped station-level annual bias directly (Figure 17). It shows the same non-monotonic elevation pattern: Srinagar, despite its low elevation, carries the second-highest bias in the basin, ahead of two higher stations. Together, the frequency/intensity diagnostics and the bias map give direct empirical support for the physical mechanisms proposed here, rather than relying on the seasonal statistics alone.

7.3 Implications of divergent seasonal and annual rankings

A central finding here is the split between seasonal and annual model rankings once we apply station-level evaluation consistently across all six stations. ACCESS-ESM1-5 and INM-CM5-0 dominate the autumn seasonal ranking, but neither holds that position annually, where MPI-ESM1-2-LR and CanESM5 lead instead (Sections 6.3 and 6.4). Similar splits between annual and seasonal model skill have been reported in CMIP6 evaluations over India and High Mountain Asia (Vinod and Agilan, 2024; Lalande et al., 2021) where annual metrics tend to favour models that get the cumulative total right even when they miss season-specific dynamics.

This has direct practical consequences. A model chosen on annual criteria alone may fit poorly for seasonal or event-based hydrological work. Examples include flood forecasting during the monsoon, snowmelt timing after winter accumulation, and drought characterisation during autumn's weak-signal transition. This mismatch can occur even if the model's annual performance looks fine. The sensitivity analysis (Section 6.5) adds a further caution. Even the annual and integrated rankings are only trustworthy at their extremes. Identifying a single best or worst model is robust. The order of models between those extremes is not, and shouldn't be treated as more precise than the data support.

7.4 Rank stability as a measure of model robustness

Quantifying rank stability explicitly across methods, stations, and seasons gives a fuller picture of model robustness than a mean rank alone. Two results anchor this. MPI-ESM1-2-LR combines the lowest integrated mean rank in the basin (4.15) with low inter-method disagreement (0.64), and stays top-ranked under 21 of 22 alternative weighting and parameter setups, including 67.5% of 2,000 randomly drawn weight vectors in the Monte-Carlo analysis. MRI-ESM2-0 shows the mirror result: the highest integrated mean rank (11.81) with the lowest inter-method variability of any model in the study (0.41) its poor performance isn't an artifact of disagreement between methods; it's a finding all three agree on. A stable best model and a stable worst model together is consistent with rank stability being reported as a meaningful robustness check in other CMIP6 evaluation studies.

Between these two extremes, the picture looks quite different from what a single consensus table suggests. Breaking rank stability into its three parts (Section 6.4) shows seasonal variability (0.96–4.62) is consistently larger than either inter-method variability (0.41–1.24) or spatial variability (0.96–3.48), for every model. Instability in this basin is mostly seasonal the three MCDM methods agree with each other far more than any model agrees with its own performance across seasons. CanESM5 and ACCESS-ESM1-5 show this most sharply, each swinging from a top-two seasonal rank to bottom-two in another season (Section 6.2). Low rank variability doesn't by itself mean high skill: MRI-ESM2-0 is both the most stable and the worst-performing model here, a combination that reporting either stability or mean rank alone would miss (Sreekumar and Valayamkunnath, 2025).

The integrated consensus framework, and its sensitivity analysis, makes this distinction between a robust extremum and a genuinely indeterminate middle order an explicit part of the evaluation rather than something buried in supporting tables. This speaks directly to recent calls for MCDM frameworks in CMIP6 evaluation to be paired with their own sensitivity account (Rahman and Pekkat, 2024), and to longstanding concerns about subjectivity in climate-model evaluation more broadly (Knutti et al., 2013). The framework here measures not just which model performs best, but how much confidence that identification deserves relative to the rest of the table.

7.5 Implications, limitations, and future directions

These findings matter directly for hydrological modeling, climate-risk assessment, and water-resources planning in mountain regions. Weak efficiency metrics persisting after bias correction mean bias-corrected CMIP6 precipitation shouldn't be used as a direct, unadjusted station-scale input to hydrological models that need day-to-day accuracy (Teutschbein and Seibert, 2012). Magnitude-based improvements cut volumetric error, but the models' inability to reproduce temporal coherence and event sequencing can still push substantial uncertainty into simulated streamflow, flood peaks, and snowmelt timing. This matters most at valley stations, where local-scale processes CMI6's native grid can't resolve have an outsized effect on hydrological response.

The uneven elevation dependence in Section 7.1, together with the strong seasonal reordering in Sections 6.2 and 7.2–7.3, argue against picking one model uniformly across the basin. Hydrological work relying on a single annually-ranked model risks getting seasonal extremes wrong most acutely in autumn, where model skill is weakest and least stable across the basin. The station- and season-specific consensus rankings developed here (Tables 5–8) give a practical basis for building multi-model ensembles tailored to a station, season, or application, rather than defaulting to one basin-wide choice.

This study has several limitations. First, it relies on a single bias-correction method Empirical Quantile Mapping. A different correction or stochastic post-processing method might change the magnitude-based results here, though it wouldn't be expected to fix the temporal-skill limits discussed in Section 7.1, since those come from free-running GCM simulation itself, not the correction method. Second, the evaluation covers six stations in the Indian portion of the Jhelum Basin. These span the basin's main elevation regimes, but may not capture its full micro-climate variability, and the elevation relationships in Section 7.1 rest on only six data points they're indicative, not definitive. Third, this evaluation covers precipitation only, not temperature, snowfall fraction, or snow-water equivalent, all important for coupled cryosphere-hydrology work in this basin. Fourth, and directly relevant to the ranking itself, the sensitivity analysis (Section 6.5) shows the consensus ranking's middle order roughly ranks 3 through 10 isn't robust to reasonable changes in criterion weighting. Only the single best and single worst model should be treated as firm conclusions; this qualification should carry into any later use of the ranking, including ensemble construction.

A further limitation is relying on a single, pre-computed bias-corrected CMIP6 product (Mishra et al., 2020), which rules out independent cross-validation split-sample or leave-1-year-out testing of the EQM step itself within this study. Any overfitting risk in the correction step comes from the source dataset, not from the station-scale evaluation framework built here. The post-correction diagnostics (Figure 5) also show both the best- and worst-ranked models under-predict high-intensity wet-day precipitation above roughly the 90th percentile, most strongly at the highest station (Gulmarg). This tail bias should be kept in mind when reading how the ranking framework treats extreme events, since it isn't weighted separately in the current MCDM criteria.

Future work should pursue dynamically downscaled or convection-permitting regional simulations, to better resolve the mesoscale and orographic processes shown here to be poorly captured at CMI6's native resolution. It should also explore process-informed or machine-learning post-processing alongside conventional bias correction, aimed specifically at improving temporal coherence rather than just magnitude. Extending the evaluation to snowfall and snow-water-equivalent variables would address this basin's cryospheric sensitivity directly. And applying the rank-stability and sensitivity framework built here to future-climate projections would test whether the same split between robust extremes and an indeterminate middle order holds up, and matters, for adaptation planning.

8 Conclusions

This study presents a station-scale, multi-method evaluation of bias-corrected CMIP6 precipitation over the Indian Jhelum River Basin. It integrates deterministic, efficiency-based, and probabilistic metrics with three independent MCDM methods (TOPSIS, PROMETHEE-II, VIKOR). An explicit rank-stability analysis and a systematic sensitivity analysis are also included. This moves beyond the basin-averaged, single-metric evaluations that typically mask local hydro-climatic variability. Bias correction substantially reduces systematic error, but does not eliminate it. Basin-wide annual RMSE ranges 9.13–14.25 mm, NRMSE 3.16–5.05, and ANMBE 1.5–35.2%. Correlation and efficiency metrics remain weak throughout: CC ranges −0.019 to 0.043, NSE ranges −1.129 to −0.423, and no station-model combination achieves positive NSE. This is consistent with bias correction improving distributional statistics without restoring day-to-day sequencing. It reflects an expected limitation of free-running historical simulations at daily resolution, not a basin-specific deficiency. Elevation modulates skill non-monotonically. NRMSE and NSE improve with elevation, while ANMBE and KGE do not. Srinagar's anomalously high bias despite its low elevation rules out a simple altitude gradient. Autumn is consistently the weakest season across four of five metrics. This reflects monsoon-withdrawal transition dynamics and orographic effects poorly resolved at CMI6's native resolution. No single model leads more than one season, underscoring the risk of annual-only model selection. The integrated consensus ranking identifies MPI-ESM1-2-LR as the most robust model and MRI-ESM2-0 as the least robust. Both results are stable under 21 of 22 alternative weighting and parameter configurations. They also hold under the large majority of Monte-Carlo weight draws: a 67.5% modal top rank and over 86% top-three placement. The remaining 11 models' relative order is not robust to reasonable weighting variation. This instability traces predominantly to seasonal disagreement rather than methodological disagreement. We present this as a substantive finding, not a limitation. It distinguishes a robust extremum from a genuinely indeterminate middle order. These results indicate that bias-corrected CMIP6 precipitation cannot be applied directly as a station-scale hydrological input in this basin. Reliable model selection requires station- and season-specific evaluation. It also requires an honest accounting of ranking uncertainty. The framework developed here is reproducible. It is extensible to other data-sparse, topographically complex basins and to future-climate projections. Future work should pursue convection-permitting regional simulations. It should also explore process-informed or machine-learning post-processing to improve temporal coherence beyond what conventional bias correction achieves.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Guotao Cui, Sun Yat-sen University, China

Reviewed by: Houssam Ayt Ougougdal, Cadi Ayyad University, Morocco

Sumana Sarkar, Fourth Paradigm Institute (CSIR), India