Hydrological model performance and uncertainty in a snow-dominated Himalayan Basin: a systematic quantitative synthesis for the Jhelum River Basin (2000–2024)

1 Introduction

Globally, hydrological systems are facing increasing stress from urbanisation, land cover change, changing precipitation patterns and increasing water demands. This has resulted in destabilising watershed responses in many regions (Yang et al., 2021; United Nations, Department of Economic and Social Affairs, Population Division, 2024). In particular, hydro-meteorological extremes and changing precipitation patterns are modifying river streamflow dynamics, which causes greater hydrological uncertainty (Mendoza et al., 2023). Mountain systems with strong cryospheric processes and altitude dependent climatic gradients are especially sensitive to these changes and are likely to exhibit strong hydrologic response (Chang et al., 2024). These changes lead to more frequent and severe floods, droughts, and water scarcity events (United Nations World Water Development Report, 2024). As a result, hydrological modelling has become essential for water resource management, adaptation, and disaster risk reduction.

Hydrological models simplify the complex water cycle by representing processes such as rainfall, snowmelt, and infiltration, and are used for water resource planning, flood forecasting, and drought management (Mosavi et al., 2018). However, these models are subject to considerable limitations in mountain basins, where complex terrain, steep elevation gradients, varied land covers, and limited observation networks compromise accuracy and reliability (Rasheed et al., 2024; Imamura, 2022). This limitation is well established in general terms, but how it plays out in basins subject to two distinct atmospheric forcing systems acting on the same catchment as in the Jhelum River Basin (JRB), has received less attention.

A central source of this limitation is precipitation forcing uncertainty. Sparse gauge networks across high-altitude catchments, combined with a scarcity of hydro-meteorological stations, mean many catchments lack direct snowfall measurements (Noor et al., 2023; Rasheed et al., 2024). To compensate, researchers commonly interpolate temperature and precipitation using lapse rates, or to use reanalysis and satellite data (TRMM data, MODIS snow cover). However, these datasets introduce systematic forcing uncertainty due to poor phase discrimination over rugged terrain (Dogra et al., 2023), while anthropogenic land-use change adds further nonlinear complexity that is difficult to simulate reliably (Ojelabi et al., 2025; Wester et al., 2019). That such constraints degrade model calibration is well established; how they specifically affect cross-model performance comparison within a single, dual-forcing basin such as the Jhelum River Basin (JRB) is not.

These forcing constraints are compounded by cryospheric change. Across the Himalaya, glacier melting has accelerated with rising temperatures, driving declining snow-cover fractions and shifts in the timing and volume of seasonal runoff (Krishnan et al., 2020; Uprety et al., 2019). Several Himalayan basins are experiencing a shift to a rainfall dominated river regime. This shift increases flood risk, reduces seasonal water storage, and introduces strong non-stationarity into hydrological systems (Pokhrel et al., 2026). In snow-fed systems specifically, cryospheric processes and temperature-driven melt behaviour largely control hydrology (Abbass et al., 2022; Chakrabortty et al., 2021; Seneviratne et al., 2021; Malik et al., 2020; Pal et al., 2021; Saha et al., 2021). Yet most of the studies reviewed here rely on temperature-index snowmelt formulations, which introduces systematic uncertainty in runoff timing and magnitude. This becomes particularly consequential under warming conditions. This structural weakness is well documented within individual models. What remains far less clear is how it manifests differently across modelling frameworks applied to the same basin, and whether the regional-scale shift toward rainfall dominance is already detectable in reported model performance, rather than only in observed streamflow trends.

Comparable challenges have been documented across other snow-dominated Himalayan basins, offering important reference points for interpreting the JRB literature. In the Koshi River Basin predictive performance has been shown to depend more on forcing data quality than on hydrological model structural sophistication (Khatakho et al., 2024). In the Upper Indus Basin, energy-balance snowmelt representations have outperformed temperature-index approaches under warming scenarios (Liu et al., 2023). Studies in the Brahmaputra headwaters confirm that structural limitations particularly the inability to represent dynamic glacier retreat are the primary source of predictive divergence under projected climate change (Mo et al., 2025; Sun et al., 2024). In Narayani-Gandaki Basin, machine learning approaches combining Artificial Neural Networks (ANN), Random Forest, and XGBoost with satellite-derived climate data have demonstrated the potential of data-driven methods as complements to physically based modelling in cryosphere-influenced systems (Suwal et al., 2025). Collectively, these cross-basin studies reveal a consistent pattern: performance divergence in Himalayan hydrological modelling is governed not by inherent model superiority but by the alignment between model structure, data availability, and simulation objective. This insight, well established in neighbouring basins, has never been systematically tested within the JRB.

The JRB represents a particularly complex and policy-relevant case. As a major tributary of the transboundary Indus system shared between India and Pakistan, the JRB is defined by steep topography and diverse microclimates (Azmat et al., 2016). Snowmelt constitutes a dominant fraction of annual runoff, making temperature the overriding climatic influence on runoff generation. Increasing temperatures, declining snowfall fractions, and earlier snowmelt onset are already shifting seasonal flow regimes (Sheikh et al., 2022), with projected warming of 1.5–3.8 °C relative to the 1986–2005 baseline by 2041–2070 under RCP 4.5 to RCP 8.5 expected to further intensify these changes (Sheikh et al., 2023; Munawar et al., 2021). Streamflow projections indicate increasing trends of 3% to over 30% across sub-basins (Saddique et al., 2019; Azmat et al., 2018), heightening both flood risk and threats to water availability.

The JRB has been modelled using a range of hydrological modelling frameworks, like Soil and Water Assessment Tool (SWAT) (Saddique et al., 2019; Yaseen et al., 2020; Saddique et al., 2020; Gull and Shah, 2021; Malik et al., 2022a) etc., the Hydrologic Engineering Center-Hydrologic Modeling System (HEC-HMS) (Mahmood and Jia, 2016; Shah and Lone, 2022a; Naqash et al., 2023; Shah and Lone, 2022b; Altaf and Romshoo, 2022; Azmat et al., 2018; Meraj et al., 2018), etc., the Snowmelt Runoff Model (SRM) (Munawar et al., 2021) etc., the Variable Infiltration Capacity (VIC) model, MIKE-based hydrologic-hydraulic systems (Altaf et al., 2020; Parvaze et al., 2022; Parvaze et al., 2023), and data-driven approaches such as ANN (Sofi et al., 2024) etc., These models vary significantly in process representation, spatial discretisation, temporal resolution, and calibration approach. SWAT is used for long term continuous simulation while, HEC-HMS is mainly applied to flood modelling (Azmat et al., 2017; Wijayarathne and Coulibaly, 2020). What is not established is whether the performance differences reported between these frameworks reflect genuine differences in predictive capability, or simply differences in the type of application each was designed for. Each of these studies reports model performance individually, typically through NSE values, but without any cross-study statistical framework to compare or contextualise results. NSE values from SWAT, HEC-HMS, and SRM applications have therefore accumulated in the literature in isolation, interpreted as indicators of model capability rather than evaluated against differences in simulation objective, temporal scale, or data availability. No prior study has performed a systematic quantitative synthesis of reported NSE values across the full range of modelling frameworks applied in the JRB.

Uncertainty treatment across this literature has received comparatively little attention. Most reported NSE values are presented as single deterministic estimates, with formal probabilistic or ensemble-based uncertainty analysis rarely implemented, despite streamflow projections for the basin varying from 3% to over 30% across sub-basins (Saddique et al., 2019; Azmat et al., 2018). This near-absence of uncertainty quantification is consistent with patterns reported in other Himalayan basins (Khatakho et al., 2024; Liu et al., 2023; Mo et al., 2025; Sun et al., 2024) but has not itself been systematically documented for the JRB, where non-stationary hydro-climatic forcing arguably makes deterministic predictions least defensible.

This body of evidence points to three critical gaps in the JRB literature. First, reported NSE values have never been statistically synthesised across the full range of modelling frameworks applied in the basin, leaving cross-study performance variability unexplained. Second, model performance has not been interpreted relative to model structure or simulation objective; NSE values are routinely attributed to model superiority rather than to differences in application context, temporal scale, or data availability. Third, formal probabilistic uncertainty analysis is largely absent from existing studies a critical deficiency given that streamflow projections for the basin vary from 3% to over 30%, rendering deterministic predictions inadequate for climate adaptation and transboundary water resource planning.

This study directly addresses all three gaps through a Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) based synthesis of 41 peer-reviewed studies (2000–2024). Four specific objectives are pursued: (i) quantitative synthesis of model performance across dominant modelling frameworks; (ii) comparison of structural characteristics and process representation in snow-dominated environments; (iii) evaluation of uncertainty treatment and propagation across the reviewed literature; and (iv) identification of future modelling pathways including hybrid and cryosphere-informed approaches for improving predictive reliability under climate change. Objectives (i) and (ii) directly address Gaps 1 and 2, respectively, by providing the first cross-study statistical performance framework and a structure-informed interpretive basis for NSE variability; objectives (iii) and (iv) together address Gap 3 by synthesising uncertainty treatment across all reviewed frameworks and identifying the methodological transition required for credible prediction under hydro-climatic non-stationarity. The paper proceeds as follows: Section 2 describes the JRB hydro-climatic context; Sections 3–4 present methodology and structural model analysis; Sections 5–6 provide quantitative performance and uncertainty synthesis; and Sections 7–9 present discussion, limitations, and conclusions.

2 Hydro-climatic and cryospheric context of the Jhelum Basin

The hydro-climatic and physiographic features of the Basin provide the background of the hydrological model choice, process modeling, and uncertainty analysis. The Basin (33,867 km²), shared between India (56%) and Pakistan (44%) and extends over 32°20′N - 35°20′N and 73°40′E - 75°35′E (Figure 1). This basin is a transboundary system with steep elevation gradients (200–6,200 m a.s.l.), with bidirectional atmospheric forcing, and strong seasonality. The basin originates from the Verinag spring in Jammu and Kashmir, flowing through the Kashmir Valley before entering Pakistan and joining the Chenab River. Major river tributaries of this basin include Kishenganga (Neelum) River, which joins near Muzaffarabad, Kunhar River, Sind, Veshaw, Sandran, Bringi, Arapath, etc. Most of these are driven by snowmelt streamflow and finally they join the Indus River.

The Basin experiences a complex hydroclimatic regime governed by the interaction of Western Disturbances and the Indian Summer Monsoon. Western Disturbances provide a substantial proportion of winter and spring precipitation, particularly at higher elevations. While monsoonal systems dominate summer rainfall and contribute to peak-flow generation. Mean annual precipitation varies spatially from approximately 700 mm at Srinagar to about 1,499 mm at Garhi Dupatta, with nearly 38% occurring during winter–spring and 62% during summer. The basin-wide mean annual temperature and precipitation are approximately 13.7 °C and 1,200 mm, respectively. Streamflow exhibits strong seasonality, with discharge beginning to rise during March–April due to snowmelt and reaching maximum levels during May–August when snowmelt and monsoon rainfall occur simultaneously. Mean annual discharge at Azad Pattan is approximately 829 m³ s⁻¹ (equivalent to 1,002 mm runoff depth), representing runoff contributions from nearly 80% of the basin area. Although glacierized areas occupy less than 1% of the Basin, seasonal snow storage and melt play a critical role in regulating spring and early-summer flows. Recent warming trends, particularly at higher elevations, have increased minimum temperatures, reduced snowfall fractions, advanced snowmelt timing, and heightened flood susceptibility while altering seasonal water availability (Sheikh et al., 2022; Azmat et al., 2016).

2.1 Seasonal runoff partitioning

The Jhelum Basin runoff partitioning is the interaction between snowmelt and precipitation seasonally. Seasonal hydrology has four unique stages where winter snowfall (December–February) maintains low baseflow levels at basin outlet. The spring snowmelt (March–May) is the main annual discharge peak in snow dominated sub-basins. This is followed by summer flows that is a combination of residual snowmelt and monsoon precipitation; and baseflow recession during autumn. In such regions, temperature governs the melt onset, hydrograph shape, and interannual discharge. Simulations studies of climate impact indicate that increasing temperature by +1.5 to +2 °C can advance peak discharge in July to May–June as the melting would accelerate (Sheikh et al., 2023; Saddique et al., 2019). The flows in winter and early spring are expected to rise and late-summer discharge might decrease because of decreased residual snowpack. Studies further reveal that increase in temperature enhances the evapotranspiration, seasonal runoff distribution, and baseflow contribution across sub-basins (Saddique et al., 2019; Hassan et al., 2024).

2.2 Observational data gaps

The limited number of hydro-meteorological observations is a structural constraint that directly limits the quality of model calibration and the quantification of uncertainty across the basin. There is also a critical scarcity of gauge stations above 3,000 m a.s.l., and most studies have used a single meteorological station to represent large elevation ranges by extrapolating the lapse rate (Mahmood and Jia, 2016; Saddique et al., 2019). To compensate for spatial observational gaps, satellite and reanalysis products such as TRMM, WFDEI, NASA POWER, and MODIS snow cover are often incorporated into modelling systems. These datasets improve the spatial coverage, but they cause uncertainties in the discrimination of precipitation phases, snowfall undercatch, and high-altitude representation (Dogra et al., 2023). More importantly, the total lack of snow water equivalent data across the basin necessitates the use of calibrated degree-day aspects as proxy measures of on-the-ground observations of the snowpack.

2.3 Climate projection spread

Climate projections for the Jhelum Basin are quite varied with respect to the selection of GCM, emission scenario and downscaling method. RCP 4.5–8.5 temperature increases of indicate temperature increases of +1.5–3.8 °C at the end of the century (Sheikh et al., 2023; Munawar et al., 2021). Further certain VIC-based models project warming greater than +5 °C under high-emission scenarios at the end of the century (Singh Jasrotia et al., 2021). The variation in the projected precipitation is between 2 and 19%, showing greater enhancement in the autumn months (Yaseen et al., 2020). The corresponding streamflow projections show an increase of about 3% to over 30%, depending on the nature of sub-basins and assumptions of the scenario (Azmat et al., 2018; Saddique et al., 2019). Figure 2 summarizes the projected hydro-climatic evolution of the Jhelum River Basin during 2020–2,100, indicating increasing trends in temperature, precipitation, and discharge, while end-of-century assessments suggest a 31% increase in winter–spring flow, a 10% decrease in summer flow, and a sediment yield of 8.54 Mg ha⁻¹ yr.⁻¹ relative to baseline conditions. This figure illustrates projected changes in temperature, precipitation, discharge, seasonal variability, and sediment yield. Moreover, assessments are based on single GCMs such as HadCM3, CanESM2, CCSM4 and NorESM1-M that do not involve ensemble averaging enhancing scenario uncertainty in predicted hydrological reactions (Mahmood and Jia, 2016). Statistical downscaling methods add an extra variability, especially in complicated landscapes where the orographic influence is challenging to address at GCM spatial scales. The hierarchical uncertainty structure scenario, model and observational that determines the predictive reliability in this system affected by the cryosphere are reflected by the wide projection spread.

3 Methodology

3.1 Search strategy and study identification (PRISMA framework)

To ensure transparency and reproducibility, this review followed the PRISMA framework. The objective was to systematically identify peer-reviewed studies applying SWAT, HEC-HMS, or comparable hydrological modelling frameworks in the Jhelum River Basin. A systematic literature search was conducted across four databases: Scopus, Web of Science, ScienceDirect, and Google Scholar. The search period was limited to 2000–2024 to capture recent advances in hydrological modelling under changing climatic conditions. The following Boolean search string was used:

(“hydrological modeling” OR “hydrological modelling”) AND (“Jhelum Basin” OR “Jhelum River”) AND (“SWAT” OR “HEC-HMS” OR “SRM” OR “VIC” OR “SWIM” OR “MIKE” OR “ANN”)

The search identified 207 records, including 60 from Scopus, 55 from Web of Science, 42 from ScienceDirect, and 50 from Google Scholar. After removing 27 duplicates, 180 studies were screened by title and abstract. Of these, 120 were excluded because they did not apply a relevant hydrological modelling framework. They reported only remote sensing or geospatial analyses without a simulation component or were purely conceptual or review-based. Sixty full-text articles were then assessed for eligibility. Nineteen were excluded due to insufficient methodological detail or non-extractable performance information. Ultimately, 41 studies were included in the qualitative synthesis, of which 34 provided sufficiently comparable methodological and/or performance data for quantitative comparative synthesis, while 7 were retained for qualitative interpretation only (Figure 3). Of the 41 included studies, 34 provided sufficiently comparable and extractable NSE-based performance data for quantitative synthesis, while the remaining 7 contributed relevant methodological or contextual insight but lacked standardised or extractable performance metrics for cross-study statistical comparison. A numerical summary of the study selection process is provided in Table 1, and the database-specific search details are summarized in Table 2. Screening and eligibility assessment were conducted by three reviewers, with disagreements resolved through discussion and consensus. Although a structured quality assessment framework was applied, no standardized risk-of-bias tool currently exists for hydrological modelling studies. Formal protocol pre-registration was not undertaken, which represents a limitation of the review.

Table 1. PRISMA-based numerical summary of study selection process.

Table 2. Database-specific search strings and number of records retrieved.

3.2 Inclusion and exclusion criteria

Clear and operational eligibility criteria were established prior to full-text assessment to ensure analytical consistency and minimize subjective selection bias. Studies were considered eligible if they reported the application of hydrological or hydro-climatic modelling frameworks within the JRB or its associated sub-basins, and addressed at least one major hydrological application domain such as streamflow simulation, flood modelling, climate change impact assessment, LULC change analysis, or snowmelt-related hydrological processes. To be retained in the review, studies were required to present an applied modelling framework and provide sufficient methodological and/or performance-related information to permit structured evidence extraction and interpretation. This included, where available, reporting on model setup, calibration or validation procedures, input datasets, and model performance indicators. However, studies were not excluded solely on the basis of incomplete reporting of standardized performance metrics if they remained relevant to the review objectives and contributed useful methodological or contextual insight.

For the purpose of quantitative comparative synthesis, a stricter set of criteria was applied. Only studies reporting extractable and sufficiently comparable model evaluation information primarily performance metrics such as Nash-Sutcliffe Efficiency (NSE) and/or coefficient of determination (R²) were included in the cross-study quantitative comparison. Where standard validation metrics were not explicitly reported but a single representative performance value was available, this value was retained with appropriate notation in the extracted dataset. For example, single-period simulations without a formal calibration-validation split such as studies reporting one representative NSE or R² value for a single hydrological year were recorded as ‘single-period’ or ‘overall’ performance rather than forced into either category, and are flagged accordingly in the extracted dataset. Where both calibration and validation performance metrics were reported, validation values were prioritised for retention because they more directly reflect out-of-sample predictive skill, whereas calibration performance can overstate model reliability due to overfitting to the parameter-tuning period. In cases where validation metrics were not available, calibration or overall reported performance values were used with appropriate caution. Studies lacking sufficiently standardized or extractable performance information were retained in the qualitative synthesis only. Studies were excluded if they were purely conceptual or review-based, lacked a primary modelling application, fell outside the geographical scope of the JRB, or did not provide sufficient methodological or performance-related information to support inclusion in either qualitative interpretation or quantitative synthesis. Given the heterogeneity in modelling objectives, temporal scales, and reporting practices across studies, these criteria were applied flexibly to preserve both methodological rigor and representativeness of the evidence base. This criterion-driven screening framework ensured that the final dataset captured both the breadth of modelling applications in the basin and the subset of studies suitable for structured comparative evaluation.

3.3 Study quality assessment

Besides eligibility screening, a systematic methodological quality review was done to assess the transparency, strength, and comparability of the included studies. A scoring system with four criteria was created on the basis of: calibration transparency, uncertainty analysis, data sufficiency, and temporal scale clarity. Each of the criteria was rated between 0 and 2 where 0 representing absent or very limited reporting, 1 representing partial reporting, and 2 representing clear and interpretable reporting. The transparency of calibration was evaluated according to the reports of the calibration and validation periods, choice of parameters, and evaluation of models. Uncertainty analysis considered whether studies included sensitivity analysis, parameter uncertainty methods such as SUFI-2 or GLUE, or predictive uncertainty bounds. Data sufficiency was evaluated based on the adequacy, completeness, and spatial representativeness of hydro-meteorological and geospatial inputs. Temporal scale clarity was assessed whether the simulation design clearly defined as event-based, seasonal, annual, or continuous long-term. The scoring criteria are summarized in Table 3.

Table 3. Methodological quality assessment framework used for included studies.

The overall quality score was between 0 and 8 and was utilized to categorize the studies according to the level of methodological strength as high (78), moderate (46), and limited (03). Overall, there was evident variance in reporting quality in the reviewed studies. Calibration design and simulation set up were more frequently reported than uncertainty treatment. Other studies used formal methods like SUFI-2, GLUE, Monte Carlo simulation or sensitivity-based uncertainty analysis. However, a substantial number of them were based on deterministic calibration without quantifying the uncertainty. Event-based flood and rainfall-runoff studies tended to exhibit high calibration transparency at the expense of uncertainty reporting. Conversely, long-term hydro-climatic and climate impact research tended to focus more on scenario design and the extent of simulation, though often failed to provide formal uncertainty treatment. This heterogeneity of methodology directly translates into cross-study comparison of reported model performance. Quality scores were used to support structured interpretation and comparison, not to exclude otherwise relevant studies. Scoring was applied only where sufficient methodological detail could be consistently interpreted from the full text.

3.4 Data extraction framework

The reviewed studies applied seven modelling frameworks: SWAT, a semi-distributed physically based model using Hydrologic Response Units; HEC-HMS, a semi-lumped conceptual model supporting event and continuous simulation; SRM, a temperature-index model estimating snowmelt across elevation zones; VIC, a physically based land-surface model with energy-balance snow representation; MIKE-based models, coupled hydrologic-hydraulic systems used for flood routing and forecasting; SWIM, a physically based continuous model structurally similar to SWAT; and ANN, a data-driven approach with no explicit process representation.

A standardized data extraction framework was developed to ensure consistency across studies and to support both qualitative synthesis and quantitative comparison. For each included study, information was extracted on model type, study region, application objective, model configuration, calibration and validation design, reported performance metrics, and input datasets. Particular attention was given to the structural and process-level characteristics of the applied models. For SWAT-based studies, extracted variables included Hydrological Response Unit (HRU) definition, land-use and soil discretization, and representation of snowmelt and elevation-dependent processes. For HEC-HMS studies, extracted information included loss methods, transform functions, baseflow representation, and routing approaches, alongside spatial configuration variables such as sub-basin delineation thresholds and elevation-band discretization where reported (Section 7.5). For other frameworks (e.g., SRM, VIC, SWIM, MIKE-based, and ANN models), equivalent model-specific information was compiled wherever reported to support cross-model comparison.

Because the Jhelum Basin is strongly influenced by snowmelt and cryospheric processes, special attention was given to the representation of snow cover, snowmelt factors, elevation zoning, and temperature-index approaches. Additional variables included the type and source of hydro-meteorological and geospatial datasets, such as meteorological forcing, snow and ice information, topographic data, streamflow observations, and climate projection inputs where applicable. Calibration strategy, validation period, and performance metrics (e.g., NSE, R², PBIAS, and RMSE) were also recorded where available. For calibration strategy specifically, each study was classified according to whether parameter estimation followed a formal probabilistic approach (e.g., SUFI-2, GLUE) or a deterministic/manual procedure, and the specific calibrated parameters reported for each modelling framework (e.g., SMFMX, SFTMP, and CN2 for SWAT; the initial abstraction ratio, storage coefficient, and recession constant for HEC-HMS; the degree-day factor for SRM) were recorded where available; the full calibration and uncertainty-treatment synthesis is presented in Section 4.4. The hydro-meteorological and geospatial datasets reported across the reviewed studies are summarized in Table 4.

Table 4. Hydro-meteorological and geospatial datasets reported across the reviewed studies.

Among the reviewed studies, SRTM (30 m) and ASTER (30 m) DEMs were the most frequently used elevation sources, together accounting for the majority of studies, followed by ALOS-PALSAR (12.5 m) in three studies and coarser-resolution USGS or MERIT DEMs in a small number of cases. For snow and cryosphere representation, MODIS-derived snow cover products (principally MOD10A1/MOD10A2) were used as an explicit validation or forcing input in roughly half of the reviewed studies, while the remaining studies relied solely on implicit temperature-index or elevation-band snowmelt representation within the hydrological model itself, without independent satellite validation. Meteorological forcing was drawn primarily from national agency records the India Meteorological Department (IMD) on the Indian side of the basin and the Pakistan Meteorological Department (PMD) or Water and Power Development Authority (WAPDA) on the Pakistani side supplemented in approximately a third of studies by global reanalysis or satellite-based products (e.g., TRMM, WFDEI, NASA POWER, GLDAS, CFSR) to fill spatial gaps in station coverage.

3.5 Quantitative synthesis approach

To complement the structural and methodological review, a comparative quantitative synthesis of reported model performance was conducted for studies with sufficiently extractable and comparable evaluation metrics. To preserve statistical independence, only one primary NSE value per study per modelling framework was retained. Where both calibration and validation values were available, validation NSE was prioritized. Where validation metrics were not explicitly reported but a single representative performance value was available, that value was retained with appropriate notation.

In cases where daily and monthly simulations were reported, daily-scale was selected to ensure there was a consistent basis of comparison. If multiple sub-basin values were available, basin-outlet performance value was applied. NSE was used as the primary performance indicator since it measures predictive performance in terms of discharge variability that is observed. As a complementary measure to determine the strength of association between the simulated and observed streamflow, the R² was retained. Since these metrics represent various variables of model performance their combined interpretation makes a strong foundation for evaluation. The number of studies available for each modelling framework varied, and comparative interpretation should therefore consider differences in sample size.

Summary statistics, including mean NSE, standard deviation, and observed performance range, were computed primarily for SWAT and HEC-HMS studies. This is because these frameworks provided the most consistently comparable data. PBIAS and RMSE were taken as supporting indicators where possible. Of the 41 reviewed studies, 35 (85%) reported R² values, while formal PBIAS, RMSE, log-NSE, and Kling-Gupta Efficiency (KGE) values were reported far less consistently and, in most cases, not at all. Of the 34 studies identified as providing extractable performance data, three (two SWAT and one HEC-HMS) were excluded from the NSE-based synthesis due to non-standardised reporting. From the remaining studies, 29 independent NSE values were retained following the prioritisation protocol described above, reflecting that some studies contributed retained values to more than one modelling framework. This confirms that R² and NSE together dominate current reporting practice in the JRB literature, while alternative diagnostic metrics remain largely absent from the reviewed record. Where possible, the interpretation of results was also done according to simulation type (event based or continuous), aim of application (historical or climate impact) and the uncertainty analysis. Since model structures, availability of data and reporting practices are heterogeneous, the findings must be viewed as an organized comparative synthesis, as opposed to a formal statistical meta-analysis.

3.6 Statistical analysis and performance variability assessment

The performance variability of extracted NSE values were analysed descriptively to facilitate performance cross-framework comparison. As the most consistently reported metric across the reviewed studies, NSE was the main foundation of comparative statistical summarization. Mean, median, standard deviation (SD), minimum, maximum and interquartile range (IQR) were calculated in each modelling framework where adequate data were available. To give normalized measure of variability, coefficient of variation (CV) was determined as (Equation 1):

C V = S D M e a n (1)

In cases where the sample size allowed, the 95% confidence intervals (CI) of the mean NSE were also calculated using (Equation 2):

C I = M e a n ± ( t × S D n ) (2)

Where, mean is used to indicate the average NSE, SD is the standard deviation, n is the number of studies and tis the t-value of the student at the 95% confidence level and n-1 degrees of freedom. These confidence intervals must be viewed not as a strict inferential limit but simply as a descriptive range of uncertainty due to the relatively small sample size and the heterogeneous reporting structure across researches. Where sufficiently reported, R² values were viewed together with NSE to promote a more holistic interpretation of model behaviour, especially in the differentiation of the strength of correlation and predictive ability. This statistical descriptive framework was applied to establish performance dispersion, consistency and relative robustness among modelling frameworks.

4 Hydrological models and structural characteristics

4.1 Evolution and distribution of hydrological modelling approaches

The 41 reviewed studies indicate a clear evolution in hydrological modelling approaches within the Jhelum River Basin. The distribution of model types derived from the extracted dataset is summarized in Table 5 and Figure 4 illustrated the temporal distribution of hydrological model applications in the Jhelum River Basin. Early studies were primarily dominated by the Snowmelt Runoff Model (SRM), reflecting the initial emphasis on temperature-index snowmelt simulation in high-altitude Himalayan catchments. Between 2016 and 2018, applications of HEC-HMS increased, particularly in flood simulation and rainfall-runoff assessment, owing to its suitability for event-based hydrological response in steep terrain. From 2019 onward, SWAT became the most frequently applied modelling framework, especially for climate change assessment, land-use land-cover (LULC) impact analysis, sediment simulation, and basin-scale hydrological evaluation. The extracted evidence also indicates increasing diversification in recent years. Studies published between 2022 and 2024 include not only SWAT, SRM, and HEC-HMS, but also VIC, MIKE-based hydrologic-hydraulic models, and data-driven approaches such as ANN. This trend reflects a gradual transition from single-process snowmelt simulation toward more integrated, multi-process modelling approaches capable of representing hydrological complexity under increasing hydro-climatic variability.

Table 5. Distribution of model types across reviewed studies.

Hydrological models range from empirical data-driven approaches, which rely on statistical pattern recognition without explicit process representation, to physically based frameworks such as SWAT, HEC-HMS, and VIC, which use mass-balance principles to explicitly parameterise land surface, terrain, and climatic variability. Model selection in data-scarce mountain basins requires balancing process realism against data and computational demands. The two most widely applied frameworks in the JRB (SWAT and HEC-HMS) are compared structurally in the following section, as they account for the majority of retained performance observations in this synthesis.

4.2 Structural comparison of SWAT and HEC-HMS

SWAT and HEC-HMS are the most widely applied frameworks in the reviewed JRB literature and form the primary basis for structural comparison in this section. SWAT functions as a physically based, semi-distributed model, whereas HEC-HMS operates primarily as a conceptual, semi-distributed framework. Their inherent differences in spatial representation, temporal simulation, and process integration predetermine their particular relevance to snow-dominated catchments, including the JRB (Azmat et al., 2017; Wijayarathne and Coulibaly, 2020).

The primary structural difference between these two models is the spatial discretisation. SWAT employs a highly granular Hydrologic Response Unit (HRU) structure which divides sub-basins by distinctive intersections of land use, soil typology and topographical slope. This explicit heterogeneous spatial parameterization greatly improves the process of localized hydrological processes, snow accumulation, infiltration, evapotranspiration, and groundwater recharge (Arnold et al., 1995; Nesru, 2023). On the other hand, HEC-HMS uses a semi-lumped sub-basin design, aggregating hydrological processes of a sub-basin, and then directing downstream discharges through a conceptual network. Although this method is effective in reducing the data requirements and computational burdens, the representation of the spatial variability which is important in topographically rugged mountainous landscapes is automatically compromised (Azmat et al., 2017).

Moreover, the models have different abilities to simulate temporal and degrees of process integration. SWAT includes an advanced water balance framework simulating long-term and continuous processes. This makes it extremely effective for analyses seasonal changes in hydrologic conditions, climate change impacts, and the interaction of land-use changes (Nguyen et al., 2020). On contrary, HEC-HMS is structurally adapted to event-based precipitation-runoff modeling, flood prediction, and peak flow definition (Wijayarathne and Coulibaly, 2020; El-Bagoury and Gad, 2024). However, it is limited in its application to large-scale, long-term assessments of the water balance because of its dependence on conceptual loss and transport mechanisms.

Regarding cryospheric dynamics, the two models traditionally use temperature-index snowmelt algorithms. However, SWAT addresses these processes at the HRU level and couple’s snowmelt with soil moisture and groundwater stores, thus defining the quantification of seasonal melt contributions. HEC-HMS only considers snowmelt as a direct macro-scale input to sub-basin runoff, but this simplifies the computation at the cost of the spatial resolution (Azmat et al., 2017).

4.3 Structural behaviour under snow-dominated hydrology

Snow and glacier melt contributes approximately 50–60% of annual runoff in upper Jhelum catchments (Marazi and Romshoo, 2018; (Dar and Ramshoo, 2013). This makes cryospheric process representation a dominant structural control on model performance across all frameworks reviewed. Explicit elevation zoning models, especially SRM and SWAT with HRU-based discretization, have been shown to be more seasonally effective. They are able to represent altitude-dependent dynamics of melt onset and snow depletion (Alam et al., 2011; Saddique et al., 2019). By comparison, semi-lumped models like HEC-HMS dilute spatial heterogeneity of snow storage, which may decrease fidelity to elevation-dependent melt variability over broad sub-basins.

One of the most fundamental weaknesses of all of the reviewed models is that they all are based on formulations of temperature-index snowmelts. Although computationally efficient and typically sufficient to calibrate historical melting these formulations do not explicitly model radiative energy balance, albedo dynamics, and turbulent heat exchange, which grow increasingly important in driving melt processes with suggested warming of 1.5–3.8 °C (Uprety et al., 2019; Sheikh et al., 2022). Such systematic uncertainty in the timing of melting and the magnitude of seasonal runoff, which cannot be resolved by calibration alone, is introduced with this structural simplification. This limitation is especially important in a non-stationary climatic environment. Simulations of climate change generally indicate earlier snow melting, higher spring discharge, and lower late-summer flows (Etter et al., 2017; Sheikh et al., 2023; Munawar et al., 2021). Divergence between model projections is, however, more of a representation of structural differences in snowmelt and runoff representation than of climate forcing differences. Models that are more explicitly coupled to the cryosphere, including VIC, show relatively better performance when applied to long-term climate forecast settings (Singh Jasrotia et al., 2021), which supports the fact that structural sufficiency of cryosphere representation is a leading predictive reliability determinant in changing conditions.

4.4 Calibration practice and uncertainty treatment

Out of the 41 reviewed studies, 11 unique studies (26.8%) implemented explicit probabilistic parameter uncertainty analysis. For the purpose of reporting the overall proportion, each study was counted only once, even where multiple model frameworks were evaluated within a single study. The remaining 30 studies (73.2%) relied on deterministic calibration, sensitivity testing, or qualitative discussion without explicit quantification of predictive uncertainty. Formal uncertainty treatment was concentrated primarily in SWAT-based applications, although isolated implementations were also observed in HEC-HMS and SRM studies. Nine of the eleven studies employed SUFI-2 within SWAT-CUP to determine parameter uncertainty, sensitivity ranking as well as predictive uncertainty quantification by using parameters like the P-factor and R-factor. These comprised nine SWAT-based studies employing SUFI-2 within SWAT-CUP (Saddique et al., 2019; Saddique et al., 2020; Gull and Shah, 2021; Malik et al., 2022b; Mattoo et al., 2023; Hamid et al., 2024; Sofi et al., 2024; Naqash et al., 2023; Hassan et al., 2024), complemented by (Azmat et al., 2018) who applied GLUE combined with Monte Carlo simulations. The distribution of calibration and uncertainty treatment across modelling frameworks is summarized in Table 6.

Table 6. Numerical summary of calibration and uncertainty treatment across reviewed hydrological models.

Counts under “Studies with Formal Probabilistic Uncertainty (n)” represent model-wise occurrences rather than unique studies. Azmat et al. (2018) evaluated both HEC-HMS and SRM and is therefore counted in both model categories, while the unique-study total remains 11. Some studies applied SUFI-2 but did not explicitly report final P-factor/R-factor values in the manuscript. *SWIM included forcing/attribution-based uncertainty assessment but not explicit probabilistic parameter uncertainty.

4.5 Model suitability and structural specialization

The Model suitability across the JRB varied strongly with simulation objective, temporal scale, and dominant hydrological process rather than with any universal model superiority. The reviewed evidence indicates clear structural specialization among modelling frameworks, with each model performing most effectively within the domain for which its architecture was designed.

HEC-HMS was most suitable for flood-event simulation and short-term hydrological assessment, particularly in steep Himalayan sub-basins. Its modular event-oriented structure, including flexible loss, transform, and routing methods, supported strong hydrograph reproduction and peak-flow simulation in flood-focused studies such as (Azmat et al., 2017; Ranjan et al., 2022; Shah and Lone, 2022a; Altaf and Romshoo, 2022). However, its simplified spatial representation and weaker treatment of long-term water balance limited its suitability for integrated basin-scale assessment. SRM was most appropriate for snowmelt-season discharge simulation, particularly in elevation-controlled snow-fed catchments where runoff is governed primarily by temperature-index melt processes and snow cover dynamics. This pattern was evident in studies such as (Romshoo et al., 2015; Azmat et al., 2016; Azmat et al., 2018; Kumar et al., 2022; Sheikh et al., 2023; Munawar et al., 2021) which showed strong suitability for seasonal snowmelt runoff estimation but more limited applicability for broader watershed-scale process integration.

SWAT showed the broadest suitability for long-term water balance analysis, climate impact assessment, and land-use change studies. Its semi-distributed HRU structure and continuous simulation framework enabled stronger representation of seasonal hydrology, groundwater response, and land-surface interactions across studies such as (Saddique et al., 2019; Saddique et al., 2020; Mattoo et al., 2023; Hassan et al., 2024; Naqash et al., 2023). VIC demonstrated strong suitability for climate projection and basin-scale hydro-climatic assessment, particularly in (Singh Jasrotia et al., 2021) owing to its physically realistic land–atmosphere coupling and energy-balance-based representation, although its application remained limited by higher data and computational demand. MIKE NAM-HD showed clear advantages in real-time forecasting and integrated hydro-hydraulic simulation, especially in (Parvaze et al., 2022; Parvaze et al., 2023; Altaf et al., 2020) where operational flow prediction and routing dynamics were central objectives. Overall, model selection in snow-dominated Himalayan basins should be guided by simulation objective, temporal scale, data availability, and process representation, not by aggregate performance statistics alone. The comparative suitability of reviewed frameworks is summarized in Table 7.

Table 7. Evidence-based model suitability and structural specialization across reviewed studies.

5 Quantitative performance synthesis

Building on the structural analysis presented in Section 4, this section provides a quantitative synthesis of hydrological model performance across the reviewed studies. Following the prioritization protocol described in Section 3, 29 independent NSE values were retained for the main synthesis, with one representative streamflow-based NSE per study per modelling framework. Where studies reported multiple variables, temporal resolutions, or calibration–validation combinations, values were selected according to the predefined hierarchy to ensure comparability and avoid pseudo-replication. Studies comparing multiple models contributed one retained value to each relevant model category. The retained dataset spans SWAT, HEC-HMS, SRM, MIKE-based models, SWIM, VIC, and ANN frameworks. The sample sizes used in the quantitative performance synthesis differ from the total model frequencies reported in Table 5 because only studies providing extractable and sufficiently comparable retained streamflow performance metrics were included. As a result, two SWAT studies (Mir and Ahmed, 2021; Akhtar et al., 2023) and one HEC-HMS study (Meraj et al., 2018) were excluded from the NSE-based synthesis. Statistical summaries are presented in Table 8.

Table 8. Summary statistics of retained NSE performance across hydrological models in the Jhelum River Basin.

5.1 NSE-based model performance comparison

The retained NSE synthesis reveals performance differences across modelling frameworks that are strongly conditioned by simulation objectives and application contexts rather than intrinsic model superiority (Figure 5). SWAT, representing the largest retained sample (N = 13), exhibited the widest performance range and greatest dispersion, with a mean NSE of 0.732 and values spanning from approximately 0.35 to 0.98. This broad spread reflects its use across diverse hydrological contexts including long-term basin-scale simulation, climate change assessment, and land-use analysis where continuous multi-process representation increases sensitivity to forcing uncertainty and parameter equifinality. The interquartile range shown in the boxplot further indicates that, despite an acceptable central tendency, a substantial proportion of SWAT applications fall within the lower performance band, consistent with the data limitations and cryosphere-related complexity of Himalayan sub-basins.

HEC-HMS (N = 8) showed a slightly higher mean NSE of 0.763 and lower dispersion than SWAT, with a comparatively compact interquartile range and fewer low-performing values. This tighter distribution is consistent with its predominant use in event-based and flood-oriented studies, where its modular rainfall-runoff structure is well aligned with short-duration hydrological response and calibration over shorter time windows. SRM achieved the highest mean NSE among frameworks with more than one retained study (mean = 0.802, N = 3), with a relatively narrow spread reflecting its structural specialization for temperature-index snowmelt simulation in elevation-controlled catchments. However, its small sample size and narrow application domain mean this result should be interpreted as indicative rather than broadly representative.

MIKE-based frameworks (N = 2, mean = 0.753), VIC (N = 1, mean = 0.790), SWIM (N = 1, mean = 0.760), and ANN (N = 1, mean = 0.700) are represented by limited observations and should therefore be interpreted descriptively only. Their retained mean or single-study NSE values all exceed the NSE = 0.65 satisfactory threshold shown in Figure 5, and three of the four exceed the NSE = 0.75 good threshold, indicating generally acceptable to good historical simulation skill within their reported application settings. Overall, the substantial overlap in NSE ranges across frameworks indicates that no single model consistently outperformed others across all use cases. Instead, observed variability particularly within SWAT appears to reflect differences in simulation objective, temporal scale, data availability, and basin complexity more than inherent model deficiency.

5.2 Supplementary synthesis of R² performance

R² values were retained using the same prioritisation protocol to assess whether correlation patterns supported the NSE-based interpretation. Only studies with explicitly reported R² values were included. Overall, the retained R² dataset indicates moderate to strong agreement between observed and simulated discharge across the reviewed studies. SWAT contributed the largest number of retained R² observations and exhibited the broadest performance range, reflecting its application across heterogeneous simulation objectives and basin conditions. HEC-HMS showed relatively high and more consistent R² values, largely corresponding to event-based and flood-oriented applications. MIKE-based, VIC, SWIM, and ANN frameworks were represented by fewer retained observations and therefore should be interpreted descriptively rather than comparatively. Summary statistics of retained R² values are presented in Table 9. Overall, the retained R² synthesis broadly supports the NSE-based findings, while also indicating that strong linear agreement does not necessarily imply uniformly high predictive efficiency across all modelling contexts.

Table 9. Summary statistics of retained R² performance across hydrological models in the Jhelum River Basin.

5.3 Structural controls, temporal scaling, and simulation objectives

The observed differences in model performance are primarily explained by structural characteristics, temporal resolution, and simulation objectives, rather than inherent superiority of any single modelling framework. Models applied at monthly or aggregated temporal scales generally exhibited higher performance due to reduced sensitivity to timing errors and short-term variability, whereas daily simulations showed greater dispersion, particularly in snow-influenced catchments. A clear distinction is observed between event-based and continuous modelling approaches. Models such as HEC-HMS, designed for event-scale simulation, tend to perform better in flood-focused applications where hydrograph shape and peak discharge are critical. In contrast, continuous models such as SWAT exhibit broader variability because they simulate a wider range of processes, including baseflow, evapotranspiration, and soil moisture dynamics over longer periods.

Process specialization also plays a key role. SRM performs well in snowmelt-dominated systems due to its temperature-index formulation, whereas SWAT provides more comprehensive basin-scale representation, integrating land surface, groundwater, and climatic processes. VIC demonstrates strength in climate-driven hydrological simulation, while MIKE-based models are better suited for operational forecasting and hydraulic routing applications. The only direct comparison between a physically based and data-driven model, as reported in the literature, indicates that ANN models can effectively capture peak flows and nonlinear patterns, while process-based models such as SWAT provide more consistent overall hydrological representation. This suggests that data-driven and physically based approaches may be complementary rather than competing frameworks. Overall, these findings confirm that model performance is strongly dependent on alignment between model structure and simulation objective, and that no single model is universally optimal across all hydrological applications.

5.4 Overall distribution and basin-wide performance pattern

The overall distribution of retained NSE values provides insight into basin-wide model performance characteristics. The histogram of retained NSE values shows that most observations are between 0.70 and 0.75. This range has the count with seven studies. Another six studies are between 0.75 and 0.80. These two ranges have 13 of the 29 retained values. This confirms that the best performance range for modelling in the JRB is between NSE = 0.70 and 0.80 (Figure 6). A secondary cluster of three studies appears in the 0.95–1.00 bin, corresponding to well-calibrated event-based simulations or studies with strong observational constraints and SUFI-2 configuration. The lower tail of the distribution, with two studies below NSE = 0.65, is associated with data-scarce high-altitude applications relying on single low-elevation meteorological stations and short calibration records. Only three studies fall below the satisfactory threshold of NSE = 0.65, while 23 of the 29 retained values exceed the good threshold of NSE = 0.75, indicating that hydrological modelling in the JRB is generally reliable across frameworks. Performance nonetheless remains sensitive to data availability, model structure, and application scale, and the spread from NSE = 0.35 to NSE = 0.98 confirms that configuration choices and calibration design drive more variability than model selection itself.

6 Uncertainty propagation and predictive reliability in the Jhelum River Basin

The JRB is influenced by interacting uncertainties in hydrological prediction because of the forcing data constraint, model structure, equifinality of parameters, and future climate scenarios. These uncertainties propagate through simulations and influence runoff timing, peak discharge, seasonal water balance, and climate sensitivity. These effects are particularly high in Himalayan basins that are dominated by snow, and whose hydrological behaviour is highly affected by elevation, snow storage, and snow melt (Mahmood and Jia, 2016; Uprety et al., 2019; Sheikh et al., 2022).

6.1 Forcing and structural uncertainty

The largest source of uncertainty is the unavailability of high-altitude hydro-meteorological observations. This forces various studies to estimate precipitation and temperature over large elevation ranges using lapse-rate extrapolation (Mahmood and Jia, 2016; Saddique et al., 2019; Marazi and Romshoo, 2018). Ultimately this affects precipitation phase, snow accumulation, and melt timing, particularly in upper catchments where snowmelt contributes substantially to annual runoff (Parvaze et al., 2022).

To reduce data gaps, a number of studies employed TRMM, WFDEI, NASA POWER, and MODIS data. These enhance spatial coverage but they also cause uncertainty because of undercatching snowfall, coarse resolution, and inability to differentiate rainfall and snowfall in mountainous areas (Dogra et al., 2023; Noor et al., 2023). These forcing errors have a direct impact on hydrograph timing and magnitude of runoff, particularly in continuous SWAT models and event-based HEC-HMS models (Hassan et al., 2024; Azmat et al., 2017). Structural uncertainty further increases predictive error. As discussed in Section 4.4, temperature-index formulations omit the radiation balance and turbulent heat exchange processes that will increasingly govern melt under projected warming, introducing structural uncertainty that no calibration procedure can eliminate (Uprety et al., 2019; Sheikh et al., 2022). Furthermore, the majority of long-term simulations assume static glacier extent and snowline position, which may overestimate future runoff under warming conditions (Azmat et al., 2018; Singh Jasrotia et al., 2021; Sheikh et al., 2023).

6.2 Parameter and scenario uncertainty

Parameter equifinality is where multiple parameters set yield statistically equivalent calibration performance but fundamentally different internal process behaviour. This represents the second major source of predictive uncertainty in the JRB (Whelan et al., 2019). In snow-fed mountain basins, runoff generation involves intensive interactions between precipitation, snowmelt, evapotranspiration, soil moisture, and groundwater recharge, making parameter compensation both structurally inevitable and difficult to detect. Across the reviewed studies, parameter selection followed three distinct approaches. The most common was auto-sensitivity ranking via SUFI-2 within SWAT-CUP, applied in nine studies including (Saddique et al., 2019; Mattoo et al., 2023; Sofi et al., 2024), where parameters were ranked by global sensitivity index before calibration sampling began. A second group relied on literature-based transfer, applying parameter ranges from prior Himalayan studies without basin-specific screening evident in (Akhter et al., 2022; Gull and Dar, 2020). A third group, predominantly HEC-HMS and SRM applications such as (Azmat et al., 2017; Mahmood and Jia, 2016), used expert judgement to define physically plausible parameter ranges, with sensitivity testing applied post-calibration rather than as a pre-calibration screening step. While SUFI-2 implements global sensitivity analysis by simultaneously sampling all parameters using Latin Hypercube Sampling, none of the reviewed studies applied a formal pre-screening step prior to SUFI-2 or GLUE sampling. Pre-screening using Morris screening (Morris, 1991) identifies and removes insensitive parameters before the main calibration sampling begins, reducing the dimensionality of the parameter space. Without this step, insensitive parameters remain in the sampling space and trade off against sensitive ones without physical justification, inflating equifinality even within a global calibration framework. The recommendation here is therefore not to replace SUFI-2 but to precede it with a Morris screening step to constrain the parameter space before global sampling begins.

Snowmelt parameters dominate JRB sensitivity analyses for clear physical reasons. SMFMX, SMFMN, SFTMP, and TIMP directly control the timing and magnitude of the largest annual discharge pulse; in a basin where snowmelt contributes 50–60% of annual runoff (Marazi and Romshoo, 2018; Dar and Ramshoo, 2013), errors in these parameters propagate through the full seasonal hydrograph. Their calibration must therefore be treated as the first priority. CN2 and SOL_AWC further compound equifinality through intercorrelation. This increased surface runoff from a higher CN2 can be offset by deeper subsurface storage from a higher SOL_AWC, producing near-identical hydrographs but fundamentally different water balance representations. One-at-a-time sensitivity testing cannot adequately identify these parameter compensations because it evaluates parameters individually and ignores interactions among them. Global sensitivity methods, such as those implemented in SUFI-2, are better suited for identifying these interactions and reducing equifinality.

In HEC-HMS, the most sensitive parameters are the initial abstraction ratio (Ia/S), the Clark unit hydrograph storage coefficient (R), and the recession constant (k). One-at-a-time testing applied in several reviewed studies (Azmat et al., 2017; Altaf and Romshoo, 2022) is insufficient to expose the Ia/S-R interaction, which jointly determines peak magnitude and timing. In SRM, the degree-day factor (DDF) is the single most sensitive parameter, exhibiting strong seasonal and elevation-dependent variability that a single basin-wide value cannot capture (Romshoo et al., 2015) demonstrated that DDF variability alone introduces ±15–20% uncertainty in peak melt volume estimates. Future HEC-HMS and SRM studies should adopt global sensitivity frameworks, consistent with the standard established for SWAT through SUFI-2.

Global sensitivity analysis via Latin Hypercube Sampling in SUFI-2 simultaneously samples all parameters and captures interactions invisible to one-at-a-time approaches. Khatakho et al. (2024) demonstrated this in the Koshi Basin using Regional Sensitivity Analysis across 1,500 Monte Carlo simulations, identifying the melt factor as the most influential parameter at 14.28% of model response, directly paralleling JRB snowmelt parameter dominance (Khatakho et al., 2024). These findings reinforce the importance of global sensitivity approaches for reducing equifinality and improving parameter identification in snow-dominated Himalayan basins. Equifinality can be reduced through several complementary approaches. One approach is multi-criteria calibration, where performance metrics such as Nash-Sutcliffe Efficiency (NSE), Percent Bias (PBIAS), logarithmic NSE (log-NSE), and Kling-Gupta Efficiency (KGE) are used simultaneously to constrain the parameter space from multiple perspectives. A second approach is physical regularisation, which eliminates parameter sets that produce unrealistic water-balance components or physically implausible hydrological behaviour. A third approach involves integrating independent observational datasets, such as MODIS-derived snow cover area (SCA) and GRACE-based terrestrial water storage anomalies. These additional datasets provide constraints beyond streamflow observations and help improve parameter identifiability. Scenario uncertainty further amplifies these parameter effects. Projected warming of +1.5 to +3.8 °C under RCP 4.5–8.5 and streamflow increases of 3% to over 30% reflect GCM and downscaling variability compounded by the forward propagation of parameter equifinality into future simulations. Model ensembles and probabilistic parameter estimation are therefore essential for generating credible future projections (Yaseen et al., 2020).

6.3 Synthesis and implications for predictive reliability

Uncertainty in the JRB propagates hierarchically from input data through parameterisation to predictive outcomes (Figure 7). The dominance of forcing and parameter uncertainties, combined with limited treatment of structural and scenario uncertainties, reveals a critical gap one confirmed by the finding that only 26.8% of reviewed studies applied formal probabilistic analysis. Among the reviewed models, SWAT showed the most consistent uncertainty treatment, particularly where SWAT-CUP/SUFI-2 was applied (Gull and Shah, 2021; Hassan et al., 2024).

Generally, HEC-HMS demonstrated weaker explicit uncertainty treatment although it was highly event-based (Azmat et al., 2017). SRM and VIC/MIKE-based applications accounted for forcing and scenario variability to some extent, but parameter and structural uncertainty remained only partially addressed (Singh Jasrotia et al., 2021; Parvaze et al., 2023). Overall, predictive reliability of the JRB cannot be evaluated based on past NSE or R² values only. Strong retrospective performance does not necessarily indicate robust future prediction in a cryosphere-influenced and climatically non-stationary basin. Model suitability is therefore objective-dependent, and future applications should increasingly adopt ensemble-based, uncertainty-aware, and cryosphere-informed modelling frameworks.

7 Discussion

7.1 Structural specialization rather than universal model superiority

This review provides the first basin-specific synthesis of hydrological modelling applications in the JRB through an integrated evaluation of model structure, reported performance, and uncertainty treatment. Rather than identifying a universally superior modelling framework, the evidence indicates that model behaviour in the basin is governed primarily by structural specialization, simulation objective, and uncertainty propagation under snow-dominated and data-constrained conditions. This interpretation builds on earlier basin-specific studies, but extends them through a comparative synthesis framework (Mahmood and Jia, 2016; Saddique et al., 2019; Shah and Lone, 2022a).

Both SWAT and HEC-HMS perform acceptably under historical conditions. HEC-HMS shows marginally higher central NSE values and lower dispersion than SWAT, reflecting its event-oriented application to flood hydrograph reconstruction and peak-flow simulation (Mahmood and Jia, 2016; Azmat et al., 2017; Altaf and Romshoo, 2022; Ranjan et al., 2022). Its modular design simulates short-duration rainfall-runoff accurately in steep terrain (Wijayarathne and Coulibaly, 2020), and event-scale calibration is by definition biased toward peak-flow accuracy. SWAT, applied to long-term climate change impact assessment, land-use change, and sediment modelling (Saddique et al., 2019; Saddique et al., 2020; Hassan et al., 2024), must simultaneously represent evapotranspiration, groundwater, and seasonal storage redistribution, increasing sensitivity to forcing uncertainty and parameter equifinality (Refsgaard and Abbott, 1996). The broader dispersion observed for SWAT is therefore structurally expected, manifested numerically in CV = 0.219 versus HEC-HMS CV = 0.159.

Within SWAT applications alone, the NSE range spans 0.35 to 0.98 a spread of 0.63 that substantially exceeds the largest mean NSE difference between any two frameworks in the retained dataset (0.07, between SWAT and SRM). This demonstrates that within-framework variability driven by configuration, data quality, and calibration design is larger than any between-framework performance gap. The productive question is therefore not which model is superior, but which configuration, objective function, and uncertainty treatment best fit a given simulation purpose. This reframes the modelling literature from a model-comparison exercise toward a configuration-and-uncertainty problem, constituting the principal scientific contribution of this synthesis.

This creates an apparent tension SRM, structurally the most criticised framework for future reliability due to its temperature-index formulation, achieves the highest historical mean NSE among reviewed frameworks (0.802). This is not a genuine contradiction strong historical calibration and structural adequacy under future non-stationary warming are different properties but the distinction is made explicit here, since historical performance metrics alone could otherwise be misread as endorsing continued reliance on temperature-index approaches.

7.2 Performance discrepancies and how this review advances prior work

Prior JRB studies have reported performance metrics without cross-study context. Mahmood and Jia (2016) present HEC-HMS results for the transboundary basin without contextualising calibration design or temporal scale. Saddique et al. (2019) report upper-band SWAT NSE values but attribute these to model capability rather than to the SUFI-2 configuration and multi-station forcing that produced them. Shah and Lone (2022a) report HEC-HMS performance for the Sindh sub-basin without elevation band configuration or uncertainty bounds. This synthesis places these findings within a broader comparative framework. The low variability in HEC-HMS performance (CV = 0.159) reflects its use in event-based simulations. It should not be interpreted as evidence of model superiority. Similarly, the highest SWAT performances were reported in studies that used SUFI-2 calibration and elevation-band discretisation. This indicates that model configuration strongly influenced the results. In addition, sub-basin studies that did not include uncertainty analysis cannot be directly compared with basin-outlet studies. Differences in spatial scale and uncertainty can introduce systematic bias. These observations do not question the validity of the individual studies. Rather, they show that cross-study comparisons can be misleading when descriptive statistical context is not considered.

Within SWAT applications (N = 13), NSE values span 0.35 to 0.98 a range driven not by model inconsistency but by application context. The lowest values (NSE < 0.55) occur in data-scarce high-altitude catchments forced by single low-elevation stations, short calibration records, and no elevation band configuration. The highest (NSE > 0.90) correspond to studies using SUFI-2 over longer records with MODIS snow cover as an auxiliary validation constraint. This directly parallels (Khatakho et al., 2024) were forcing data quality alone shifted NSE from −0.20 to 0.91 within the same model structure confirming that reported NSE values are co-determined by data quality, calibration design, and application context, not model capability alone. Prior JRB studies have presented NSE values in isolation, creating a misleading impression of cross-study comparability. The present synthesis exposes this: with CV = 0.219, the performance range within SWAT alone exceeds the mean difference between frameworks. The productive question is therefore not which model is best, but which structure is most appropriate for a given objective, data context, and temporal scale. For HEC-HMS (mean NSE = 0.763; SD = 0.122), the tighter distribution reflects its narrower event-based application domain not inherent superiority over SWAT.

As a systematic synthesis rather than a primary modelling study, this review does not propose a new modelling approach. The contribution is methodological and interpretive: by establishing the first cross-study descriptive statistical framework for the JRB, this synthesis advances the basis on which model outputs are evaluated and compared shifting the modelling literature from isolated performance reporting toward a configuration- and uncertainty-aware standard. Improvements in accuracy, applicability, and process representation will follow when future JRB studies adopt the configuration practices, uncertainty quantification standards, and ensemble-based frameworks identified through this synthesis.

7.3 Cryosphere influence as a first-order control on model behaviour

In the Jhelum Basin, hydrological response is highly controlled by cryospheric processes and steep elevation gradients. Thus, snow and glacier melt represent a significant fraction of seasonal runoff in upper sub-basins (Dar and Ramshoo, 2013; Parvaze et al., 2022) which increases the structural implications of model. Discretization models that maintain the spatial heterogeneity via elevation zoning or HRU are more appropriate to seasonal snowmelt dynamics. That is why SWAT and SRM are still used in the studies devoted to the sensitivity of climate and the variability of seasonal runoffs (Saddique et al., 2019; Azmat et al., 2017; Sheikh et al., 2023). Their advantage lies in structural alignment with basin processes but not necessarily higher performance metrics. Conversely, the semi-lumped models like HEC-HMS make spatial variability simpler and lump runoff at the sub-basin level (Azmat et al., 2017). Although this simplification works well in event-based flood simulation but there are limitations. It is not sufficient to model snow accumulation and melt over long periods of time at the elevation. This constraint is significant in Jhelum River Basin. Seasonal snow storage and melt in the basin have a significant influence on streamflow. So, the capacity to simulate snow pack dynamics and hydrological memory is a key consideration in the model selection.

7.4 Model setup, calibration, and uncertainty communication

In most JRB studies, NSE is used both for calibration and performance evaluation. NSE can be used to evaluate the overall performance of a model. However, it is strongly influenced by peak flows. Consequently, models calibrated from NSE only may not be able to accurately simulate low-flow conditions. These are critical times for water supply and groundwater recharge. Given the need to use multiple performance metrics, future research should be conducted using multiple performance metrics like NSE, PBIAS, log-NSE and KGE. Together these metrics are will provide more comprehensive evaluation of model performance.

It is also important to consider the length of the calibration record. Short records might not fully reflect the variability of snowpack and monsoon. This is an especially difficult problem with records that are under five years in length. Khatakho et al. (2024) observed that longer calibration periods resulted in better simulations for Koshi basin. The study suggested differential split-sample validation and a minimum of 10 years of data. Similar practices should be used in future JRB studies. Further the default station assignment method could cause biases in mountainous areas. A station at low elevation can be used to represent the whole elevation band. This may result in temperature and precipitation inaccuracies at higher elevations. Khatakho et al. (2021) suggested a point-pixel approach using satellite estimates. Future studies should be conducted to compare both approaches within the JRB. Another source of uncertainty is the choice of the baseline period. Climate-change signals may vary depending on the baseline period. This is because the warming rate is different at different elevation bands (Sheikh et al., 2022).

There should be more consistency in reporting uncertainty. The p-factor and the r-factor should always be given. The forcing data set used must also be specified. Khatakho et al. (2024) demonstrated that the uncertainty range can significantly differ depending on the precipitation product used. Evin et al. (2026) further emphasised the importance of clearly communicating uncertainty to decision-makers. Structural uncertainty can be practically handled by using the Bayesian Model Averaging (BMA). BMA is a combination of various models’ predictions. These may include SWAT, HEC-HMS, and SRM. Each model is weighted according to its performance. Khatakho et al. (2026) suggested this for the Himalayan basins. In addition to process-based models, machine-learning techniques can be used. Suwal et al. (2025) got the NSE value of 0.92–0.93 with ANN, RF and XGBoost models. Post-processing of the SWAT snow simulations with LSTM in the JRB could be of special interest.

Parameter regionalisation has not received much attention in the JRB. This is particularly unexpected as many catchments at high elevations are still ungauged. There are three basic methods. These include regression-based regionalisation, spatial proximity transfer and physical-similarity methods. Khatakho et al. (2026) concluded that the transferability of the parameters was not reliable without an appropriate representation of the elevation bands. It is therefore important to focus on physical-similarity methods in future studies. Matching criteria can be snow-cover from Moderate Resolution Imaging Spectroradiometer (MODIS) and hypsometric characteristics.

Finally, a comprehensive uncertainty framework should address four levels of uncertainty. First, forcing uncertainty can be reduced through data assimilation. One possibility is to use Ensemble Kalman filtering with MODIS and Gravity Recovery and Climate Experiment (GRACE) observations. Second, the uncertainty of the parameters can be measured using the Bayesian calibration approach. This method produces distributions of parameters rather than best estimates. Third, structural uncertainty can be dealt with via BMA. Fourth, residual timing errors can be minimized by post-processing with LSTM. This is particularly relevant for snowmelt timing. Adoption of this framework would strengthen uncertainty assessment in the JRB. It would also provide a stronger basis for climate adaptation planning.

7.5 Guidance for model configuration in snow-dominated Himalayan Basins

The performance variability documented in Section 5, including the NSE range of 0.35 to 0.98 within SWAT applications, is partly due to differences in model configuration. Studies that used elevation-band discretisation, multi-station forcing, and auxiliary snow-cover validation generally achieved higher performance. In contrast, studies that relied on single-station lapse-rate forcing and short calibration periods were concentrated in the lower performance range.

For basin-scale simulations in steep Himalayan catchments a drainage area threshold of 100–300 km² is typically appropriate, consistent with published SWAT configuration guidance for mountain systems (Arnold et al., 2012) and with the sub-basin configurations reported across the reviewed JRB studies, where effective applications used thresholds in this range (Saddique et al., 2019; Malik et al., 2022a). Thresholds below 100 km² produce excessive sub-basin fragmentation that increases parameter dimensionality without proportionate improvement in process representation, while thresholds above 300 km² lump heterogeneous elevation zones, diluting the spatial variability that drives differential snowmelt behaviour across the JRB’s 200–6200 m elevation range.

Similarly, HRU area thresholds of 2–5% for land use and soil type, and 5–10% for slope class, are recommended for Himalayan catchments with high land-use and land-cover heterogeneity. These recommendations are consistent with the sensitivity analyses of Arnold et al. (2012) and findings from high-altitude studies in the Koshi Basin. Thresholds greater than 10% can remove hydrologically important land-cover classes from the model. These include glaciers, permanent snowfields, and sparse alpine vegetation. Although these classes occupy relatively small areas, they contribute significantly to snowmelt runoff and groundwater recharge in the upper Jhelum River Basin.

Further the JRB spans 200 m to over 6,200 m a.s.l., a minimum of five to ten elevation bands at intervals of 200–500 m within each sub-basin is recommended for catchments spanning more than 2000 m of relief. In SWAT, elevation band-specific lapse rates should be calibrated separately rather than applied as single basin-wide constants. Model evaluation should also extend beyond streamflow-based metrics. MODIS-derived snow cover area products (e.g., MOD10A1 at 500 m resolution) provide a spatially distributed observational constraint on simulated snowpack independent of the streamflow record. Comparing simulated SCA depletion curves against MODIS SCA during the spring melt period provides critical process-level validation of snowmelt timing. Khatakho et al. (2024) demonstrated this approach in the Koshi Basin, recommending multi-response validation using modelled snow cover against remote sensing estimates as a standard step for uncertainty reduction in Himalayan basins.

7.6 Structural limits of temperature-index snowmelt modelling

The structural limitations of temperature-index snowmelt formulations are discussed in detail in Section 4.4. One additional concern specific to long-term projection studies is that the majority of reviewed models assume static glacier extent, which is inconsistent with observed cryospheric retreat across the JRB and may systematically overestimate future runoff under sustained warming. This distinction between calibration adequacy and structural adequacy is particularly consequential in cryosphere-dominated basins where process representation governs predictive reliability under non-stationary conditions.

7.7 Uncertainty as a hierarchical system property

In Jhelum Basin, hydrological modelling has uncertainty due to many interacting layers such as forcing uncertainty, structural simplification, parameter equifinality and variability in scenarios. These uncertainties propagate hierarchically rather than independently. Forcing uncertainty is especially significant due to the scarcity of high-altitude meteorological observations. Various studies rely on lapse-rate extrapolation of temperature and precipitation (Mahmood and Jia, 2016; Saddique et al., 2019). To compensate this Satellite-derived products such as TRMM and MODIS are used. This improves spatial coverage but introduce uncertainty in precipitation phase discrimination and snowfall estimation (Dogra et al., 2023).

These limitations are also complicated by structural uncertainty. Structural simplification, including semi-lumped spatial discretisation, reduces process fidelity under varying climatic conditions (Uprety et al., 2019). Parameter equifinality, especially in continuous models, allows multiple parameter sets to produce similar calibration performance but divergent predictive outcomes (Whelan et al., 2019). However formal probabilistic uncertainty analysis was applied in only 26.8% of the reviewed studies, and SUFI-2 and GLUE were primarily used in SWAT-based applications. SWAT-based studies are mainly determined by the application of such techniques as SUFI-2 and GLUE (Saddique et al., 2019; Naqash et al., 2023), and HEC-HMS and SRM applications mostly use deterministic calibration (Mahmood and Jia, 2016; Shah and Lone, 2022a). The uncertainty on scenarios due to GCM choice and downscaling further increases predictive variability (Yaseen et al., 2020; Munawar et al., 2021). Such results suggest that uncertainty is not just a methodological weakness but a major feature of hydrological forecasting in cryosphere-impacted basins.

7.8 Broader implications: comparison with other Himalayan Basins

The structural patterns, performance characteristics, and uncertainty challenges documented in the JRB are broadly consistent with findings from other snow- and glacier-dominated Himalayan basins. This reinforces the conclusion that the JRB represents a wider class of hydrologically complex, data-scarce mountain systems. This consistency extends beyond the Himalayan region: studies in the Alps and Andes have similarly demonstrated that continuous distributed models outperform event-based models for long-term water balance and climate estimation, while event-based models retain advantages for flood simulation confirming that the principle of structural specialisation governing model performance in the JRB is basin-independent (Portuguez-Maurtua et al., 2023; Calizaya et al., 2021).

In the Koshi River Basin of Nepal (Khatakho et al., 2024) applied the distributed J2000 model. Despite the model’s more physically explicit cryospheric representation relative to temperature-index frameworks used across JRB studies, calibration performance varied dramatically with forcing data: the combined gauge-APHRODITE dataset achieved NSE = 0.91, while MSWEP under ground-based calibration yielded NSE = −0.20, a spread larger than the entire cross-framework variability observed in the JRB synthesis. This reinforces the conclusion that forcing data quality is the primary determinant of predictive performance in data-scarce Himalayan basins. Further, GLUE-based uncertainty analysis confirmed that the proportion of observed data within the 95th percentile uncertainty band ranged from 41 to 72% solely depending on which precipitation product was used directly paralleling the uncertainty challenges documented across the JRB literature. In the Karnali River Basin, Khatakho et al. (2021) similarly evaluated satellite precipitation products using J2000, demonstrating that TRMM performed best under gauge-based calibration while MSWEP showed the greatest improvement under product-specific recalibration. A subsequent study evaluated SWAT parameter transferability between the Koshi and Karnali basins under 16 modelling scenarios (Khatakho et al., 2026). The results showed that incorporating elevation bands and combining gauge observations with CHIRPS precipitation data significantly improved model performance. These modifications also enhanced cross-basin parameter transferability. The authors further concluded that spatial parameter transferability is not viable in large, heterogeneous Himalayan basins unless elevation band configuration is included.

In the Narayani-Gandaki Basin of Nepal, Suwal et al. (2025) demonstrated that ensemble machine learning models like ANN, RF, and XGBoost achieved NSE = 0.92–0.93. ANN performed best for high flows while RF outperformed for low flows. This confirms that ensemble ML approaches are necessary to reproduce the full seasonal discharge range. The successful transfer of the ANN framework to the Trishuli tributary further confirmed the cross-basin applicability of data-driven approaches where physically based models are constrained by data scarcity. In the Tamor River Basin of eastern Nepal, HBV and VIC models applied demonstrated that temperature-index snowmelt formulations produce systematic melt-timing errors of 2 to 4 weeks under warming scenarios, with early-season discharge overestimated and late-season discharge underestimated a bias pattern structurally identical to that documented in JRB SWAT applications (Normand et al., 2010; Latif et al., 2020). Garee et al. (2017) evaluated the VIC model in the Hunza and Gilgit sub-basins of the Upper Indus Basin. Their results showed that energy-balance-based snowmelt representation reduced seasonal runoff timing errors compared with temperature-index approaches. The advantage became more pronounced under high-emission warming scenarios, where turbulent heat flux contributed substantially to snowmelt processes. Collectively, these findings confirm that the structural limitation of temperature-index formulations is basin-independent across the Hindu Kush-Himalaya system not a JRB-specific modelling artefact and that the transition toward energy-balance representation is a regional priority, not merely a JRB-specific recommendation.

Taken together, the cross-basin evidence confirms that structural adequacy, uncertainty transparency, and ensemble-based prediction must replace single-model deterministic assessment as the standard of scientific practice in Himalayan basin modelling. Comparable syntheses for Afghan Indus tributaries including the Kabul River Basin remain absent from the literature, representing a further regional gap that limits generalisability across the full Hindu Kush-Himalaya system. It is worth noting that cross-basin evidence does not point to a single dominant driver of predictive performance: while energy-balance representation improved timing accuracy in the Upper Indus Basin (22), forcing-data quality was shown to be the primary determinant of performance in the Koshi Basin regardless of model structure (21). Both factors are therefore relevant, and their relative importance may vary by basin and data context.

The structure-informed interpretive framework developed here links reported model performance to model design, data limitations, and the goal of each simulation. It does not treat performance as a fixed property of the model itself. This framework can be directly applied to other data-scarce, snow- and glacier-dominated basins. These include basins across the Hindu Kush-Himalaya and beyond, such as the Indus, Ganges, and Brahmaputra headwaters, as well as similar mountain systems like the Andes. Future studies in these basins can avoid the same recurring problems found in the JRB literature. This means adopting better practices from the start. Uncertainty methods should be planned before calibration begins, not added afterward. P-factor and R-factor values should be reported alongside NSE and R² in every study. Elevation-band or HRU-based spatial breakdowns should be used instead of relying on a single low-elevation weather station. Snowmelt timing should be checked against independent satellite data, such as MODIS snow cover. Making these practices standard, rather than optional, would make it much easier to compare results across studies. It would also give policymakers better evidence when managing water systems that cross borders and depend on melting snow and ice. Other cross-border river basins face the same kind of shifting, less predictable climate. For these basins too, using probabilistic, ensemble-based forecasts instead of single-number predictions would strengthen the evidence behind water-sharing agreements and infrastructure planning. This mirrors the role such approaches should play for the Indus Waters Treaty in the JRB itself.

8 Study limitations, future research and policy implications

8.1 Study limitations

Although a comprehensive synthesis presented in this review, several methodological and evidence-related limitations should be acknowledged when interpreting the findings.

• Small per-framework sample sizes: SWAT (N = 13) and HEC-HMS (N = 8) support descriptive comparison, but SRM (N = 3), VIC, MIKE-based, SWIM, and ANN (N = 1–2) are insufficient for statistical inference. SRM has a mean NSE of 0.802 which is based on three studies and is likely to change considerably with a single new observation.
• Reporting heterogeneity: Studies differ in whether calibration or validation values are reported, temporal resolution (daily vs. monthly), and spatial scope (outlet vs. sub-basin). The prioritisation protocol reduces but does not eliminate this heterogeneity, meaning the performance spread within frameworks likely understates true cross-study variability.
• Absence of a standardised risk-of-bias tool: No validated quality assessment framework exists for hydrological modelling systematic reviews. The four-criterion scoring system applied here calibration transparency, uncertainty analysis, data sufficiency, and temporal scale clarity provides a structured approximation but introduces subjectivity that a community-agreed standard would reduce.
• CMIP5/RCP reliance: The reviewed literature predominantly uses CMIP5-era GCMs and RCP scenarios. CMIP6 models generally project stronger warming over the Hindu Kush-Himalaya under comparable emissions pathways. The reported warming spread of +1.5–3.8 °C may therefore underestimate the upper bound, and corresponding streamflow projections may require upward revision as CMIP6-based JRB studies become available.
• Streamflow-centric performance synthesis: The synthesis evaluates NSE and R² for streamflow only. Model performance for sediment yield, water quality, groundwater recharge, and evapotranspiration all directly relevant to JRB water resource management remains unassessed. Models performing acceptably for streamflow may perform poorly for these outputs, and comprehensive evaluation would require multi-variable validation.

These limitations directly inform the future directions in Section 8.2: Limitation 1 motivates the model selection priority; Limitation 2 motivates standardised uncertainty reporting under Uncertainty Quantification; Limitation 4 motivates updating the synthesis as CMIP6-based JRB studies emerge; and Limitation 5 motivates expanding future assessments to sediment yield, water quality, groundwater recharge, and evapotranspiration.

8.2 Future research directions

While grounded in JRB-specific evidence, the recommendations below generalise directly to other data-scarce, cryosphere-dominated basins facing comparable structural and observational constraints, and are intended to help future studies in those basins avoid the same reporting and uncertainty gaps identified here. Hydrological modelling in the JRB has reached a methodological inflection point. Current frameworks provide useful historical simulation but fall short under accelerating cryospheric change, constrained by structural simplifications, observational gaps, and near-absent uncertainty treatment.

• Model Selection (Near-term): Model choice must follow the simulation objective, not convention. Event-based frameworks (HEC-HMS) suit flood hydrograph reconstruction and peak-flow estimation; continuous physically based models (SWAT, VIC) are better suited to climate sensitivity, land-use change, and long-term water balance studies. A formalised model-suitability framework linking model structure to management objective would support more principled, reproducible decisions.
• Cryosphere Representation (Medium-term): Virtually all reviewed models rely on temperature-index snowmelt schemes whose structural limitations are detailed in Section 4.4. Under projected warming of 1.5–3.8 °C, transition to energy-balance snowmelt methods, dynamic snowline tracking, and glacier mass-balance coupling is not optional it is a prerequisite for credible future projection in a basin where melt onset is already advancing.
• Uncertainty Quantification (Near-term): Only 26.8% of reviewed studies applied formal probabilistic uncertainty analysis. Addressing this requires a hierarchical approach:

• Forcing: Assimilate MODIS snow cover and GRACE water storage anomalies via ensemble Kalman filtering to continuously correct model states rather than relying on static calibration
• Parameter: Apply MCMC-based Bayesian calibration for true posterior distributions; use Bayesian Model Averaging (BMA) to combine SWAT, HEC-HMS, and SRM outputs as recommended in the discussion section.
• Structural: Deploy Physically Informed Machine Learning (PIML), feeding SWAT snowpack and soil moisture states into LSTM layers to learn runoff timing corrections
• Always report P-factor and R-factor alongside NSE and R², specifying the forcing dataset used

• Data Integration (Medium- to Long-term): Predictive reliability is ultimately bounded by observational scarcity. High-altitude measurements above 3,000 m a.s.l. remain critically absent the dominant source of performance variability across reviewed literature. Expanding gauge networks at elevation would directly reduce forcing uncertainty. Remote sensing should move beyond gap-filling into systematic data assimilation and multi-source validation workflows.
• Hybrid Modelling (Medium-term): Coupling physically based models with machine learning offers the most promising path for reducing both structural and parameter uncertainty. Suwal et al. (2025) achieved NSE = 0.92–0.93 in the Narayani-Gandaki Basin using ensemble ML (ANN, RF, XGBoost), with RF excelling for low flows and ANN for high flows. In the JRB, LSTM post-processing of SWAT snowpack states during spring melt is the highest-priority near-term application. These methods must complement physical reasoning, not replace it. Beyond LSTM, convolutional LSTM (ConvLSTM) architectures which integrate spatial convolutional operations with temporal recurrence offer a promising avenue for spatially distributed snowmelt correction across the JRB’s heterogeneous elevation zones, enabling simultaneous learning of spatial snow cover patterns and melt timing dynamics.
• Based on the structural gaps identified in this synthesis, three model directions are specifically recommended for future JRB and analogous Himalayan basin studies: energy-balance or hybrid energy-balance/temperature-index snowmelt models to address the structural snowmelt-timing limitation documented across all reviewed frameworks; LSTM or ConvLSTM architectures coupled with SWAT or HEC-HMS outputs to correct residual timing and low-flow errors; and Bayesian Model Averaging (BMA) across SWAT, HEC-HMS, and SRM outputs to produce ensemble, uncertainty-quantified streamflow projections rather than single deterministic estimates.

8.3 Policy and governance

• Indus Waters Treaty: Projected streamflow increases of 3–30% across sub-basins reflect genuine structural and scenario uncertainty that single deterministic projections cannot capture. Treaty assessments must report probabilistic bounds derived from multi-model ensembles, not point estimates. Critically, water allocation calculations under the Treaty should adopt a probabilistic allocation framework one that defines permissible abstraction volumes as a function of exceedance probability rather than mean annual flow, so that allocations remain valid across the full projection range. Ensemble-based probabilistic forecasting frameworks have demonstrated improved allocation reliability under non-stationary climatic conditions relative to static deterministic approaches (Yaseen et al., 2020), a finding directly applicable to the Jhelum-Chenab system.
• Reservoir Operations and Water Management: Future reservoir operation in the Jhelum River Basin should incorporate ensemble-based inflow forecasting and adaptive rule curves. Existing operating rules were developed under the assumption of stationary seasonal inflow patterns. However, projected warming of 1.5–3.8 °C is expected to advance snowmelt onset and shift peak inflows from July toward May–June. As a result, historical rule curves may no longer be optimal under future hydro-climatic conditions. Reservoir operators should therefore adopt adaptive management approaches that account for changing inflow timing and uncertainty. Robust decision-making frameworks and probabilistic forecasting can improve reservoir performance under non-stationary climate conditions.
• Flood Warning: Probabilistic peak discharge distributions enable more actionable, risk-proportionate civil protection than single-point estimates. Ensemble flood forecasting systems already operational in the European Flood Awareness System (EFAS) and the Global Flood Awareness System (GloFAS) provide templates directly transferable to the JRB operational context (Alfieri et al., 2013).
• Climate Resilience and Sustainable Planning: Building climate resilience in the JRB requires moving from deterministic infrastructure design toward scenario-neutral approaches that evaluate system performance across a wide range of possible futures rather than a single design scenario. For water supply planning, this means dimensioning storage and distribution infrastructure to remain functional across the full streamflow projection envelope (3–30% increase) rather than for a central estimate. Seasonal climate forecasts with quantified uncertainty, when integrated into reservoir operation models, have demonstrated 10–15% improvements in water supply reliability in analogous mountain systems.
• Transboundary Governance: High-altitude data gaps partly reflect limited hydro-meteorological data sharing between India and Pakistan. Joint monitoring protocols, open data access agreements, and shared ensemble modelling platforms are inseparable prerequisites for reliable prediction and responsible water governance across the transboundary system.

9 Conclusion

This study presents the first systematic quantitative synthesis of hydrological modelling in the Jhelum River Basin, drawing on 41 peer-reviewed studies (2000–2024) through a PRISMA-based protocol. The central finding is that performance divergence across frameworks is governed by structural specialisation and simulation objective, not universal model superiority. With within-framework NSE variability in SWAT (CV = 0.219) exceeding the mean performance difference between any two frameworks in the retained dataset, the productive question is not which model is best but which configuration, objective function, and uncertainty treatment are most appropriate for a given simulation purpose and data context. SWAT (mean NSE = 0.732; N = 13) exhibits broad dispersion consistent with long-term multi-process simulation; HEC-HMS (mean NSE = 0.763; N = 8) shows tighter distribution reflecting its event-based flood domain; and SRM achieves the highest mean NSE (0.802; N = 3) within its narrow snowmelt-specific scope. Complementary R² synthesis confirms moderate to strong agreement across all frameworks while reinforcing that correlation strength does not guarantee predictive efficiency. Three systemic constraints limit predictive reliability and cannot be resolved through calibration alone. First, universal reliance on temperature-index snowmelt formulations whose structural limitations are documented in Section 4.4 means acceptable historical calibration does not guarantee future reliability under projected warming of 1.5–3.8 °C. Second, near-total absence of high-altitude observations forces lapse-rate extrapolation that propagates precipitation phase and timing errors through every simulation stage a structural forcing problem no parameter optimisation can fully compensate. Third, only 26.8% of reviewed studies implemented formal probabilistic uncertainty analysis. In a basin where streamflow projections range from 3% to over 30% across sub-basins, this is not a methodological shortcoming but a policy risk with direct consequences for Indus Waters Treaty governance and Mangla Dam operations.

This review makes three principal contributions. It establishes the first structured basin-specific performance database for the JRB, enabling comparative statistical analysis unavailable from individual studies or prior narrative reviews. It develops a structure-informed interpretive framework linking NSE variability to model architecture, simulation objective, and data constraints advancing beyond context-free performance reporting toward a process-oriented basis for model selection. And it constructs a hierarchical uncertainty synthesis that exposes systemic reliance on deterministic calibration as a basin-wide methodological gap with direct consequences for climate adaptation planning and transboundary water governance. The overarching directive is clear: the JRB must move from isolated deterministic modelling toward ensemble-based, uncertainty-aware, and purpose-specific frameworks. Physically informed machine learning, energy-balance snowmelt representation, and Bayesian uncertainty quantification are not aspirational additions they are prerequisites for credible hydrological prediction in a transboundary, cryosphere-dominated, and climatically non-stationary basin. This synthesis provides the first structured foundation for that transition.

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: Amar Deep Tiwari, Michigan State University, United States

Reviewed by: Rajesh Khatakho, University of Georgia, United States

Romana Jamshed, COMSATS University Islamabad, Pakistan

Yonas Abebe Balcha, University of Tartu, Estonia